Seth,
this problem is most likely related to LAPACK's SLAMCH and DLAMCH
routines. If these routines get translated with f2c the compiler can
optimize things away and create infinite loops. This functions get
called when you start using singular and eigenvalue routines. It's
because Fortran 77 didn't have a way to specify "volatile" storage.
In any case, you can avoid the problem by linking different
SLAMCH and DLAMCH routines. I'm attaching two possible
impementations that use C99 and Fortran 95. Probably you'll have to
ensure that name mangling occurs correctly (rename slamch to slamch_
and dlamch to dlamch_).
Piotr
On Tuesday 19 September 2006 03:26, Dr. Seth Olsen wrote:
> Hi Travis,
>
> I would very happily accept the Scientific and MMTK patches. Thank
> you very much for the offer.
>
> I hadn't thought much about it until very recently, when the advent
> of a new IT structure in our laboratory forced me to upgrade my
> existing software. It was then that it became obvious that the LAPACK
> routines called by Numeric and MMTK were refusing to work. I've
> spent the day trying to wrestle with the problem (I am by no means an
> expert). I finally decided to get around the problems by altering
> the problem-solving routines in MMTK by using routines from numpy and
> then immediately applying a=Numeric.array( a.tolist()), which has
> stopped my script from gagging (although I have still to demonstrate
> to myself that everything is working). The uses to which I apply
> MMTK usually mean that the matrices in question are pretty small, so
> I don't have to worry too much about speed.
>
> I was suprised to note, however, that a straightforward application
> of F2PY to some fresh downloaded LAPACK F77 code did not work, and
> sent my machine into a similar endless whirr as I had seen at the
> beginning of the day. Apparently there's something sinister going
> on...
>
> Cheers,
>
> Seth
>
> On 9/19/06, Travis Oliphant <[EMAIL PROTECTED]> wrote:
> > Dr. Seth Olsen wrote:
> > > Hi Bill,
> > >
> > > MMTK has not made the conversion over to the new numpy module.
> > > It is built against the old Numeric code, and the word from its
> > > developers is that changing to numpy cannot be a priority now.
> >
> > My suggestion is to *kindly* put pressure on them.
> >
> > I've spent at least a hundred hours making it easy for people to
> > port Numeric and Numarray-built code to NumPy. Because of this,
> > I'm a little bit frustrated by this kind of response. I understand
> > it will take time for people to migrate, but it really does not
> > take that long to port code to use NumPy.
> >
> > I've offered to do it for any open source code. In fact, I just
> > spent 30 minutes and ported both Scientific Python and MMTK to use
> > numpy. I'll send you a patch if you want. It is true, that the
> > result needs to be better tested, but I'm confident that any errors
> > which might remain in the compatibility layer will be easily
> > fixable (and we need people who are willing to do the tests to fix
> > them).
> >
> > I'd rather not do this, but if necessary we can easily create an
> > SVN tree of third-party packages ported to use NumPy if the
> > package-owners are not willing to do it. Keeping Numeric packages
> > around except for legacy systems will only make things harder.
> >
> > I'll repeat the same offer I've made before: I will gladly give my
> > book and my help to any open source library author who will make
> > porting to NumPy a priority for their package. Note, however, my
> > (free) ports to use NumPy do not use any "numerix-style" layer.
> > The library is converted to work with NumPy alone. In other words,
> > I won't spend any more "spare" time supporting 3 array packages.
> >
> > Best regards,
> >
> > -Travis
> >
> >
> > -------------------------------------------------------------------
> >------ Take Surveys. Earn Cash. Influence the Future of IT
> > Join SourceForge.net's Techsay panel and you'll get the chance to
> > share your
> > opinions on IT & business topics through brief surveys -- and earn
> > cash
> > http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=
> >DEVDEV _______________________________________________
> > Numpy-discussion mailing list
> > Numpy-discussion@lists.sourceforge.net
> > https://lists.sourceforge.net/lists/listinfo/numpy-discussion
#include <float.h>
#ifndef FLT_DIGITS
#define FLT_DIGITS 24
#endif
#ifndef DBL_DIGITS
#define DBL_DIGITS 53
#endif
float
slamch(char *cmach) {
char ch = cmach[0];
float sfmin, small, one = 1.0, zero = 0.0;
if ('B' == ch || 'b' == ch) {
return FLT_RADIX;
} else if ('E' == ch || 'e' == ch) {
return FLT_EPSILON;
} else if ('L' == ch || 'l' == ch) {
return FLT_MAX_EXP;
} else if ('M' == ch || 'm' == ch) {
return FLT_MIN_EXP;
} else if ('N' == ch || 'n' == ch) {
return FLT_DIGITS;
} else if ('O' == ch || 'o' == ch) {
return FLT_MAX;
} else if ('P' == ch || 'p' == ch) {
return FLT_EPSILON * FLT_RADIX;
} else if ('R' == ch || 'r' == ch) {
return FLT_ROUNDS < 2 ? one : zero;
} else if ('S' == ch || 's' == ch) {
/* Use SMALL plus a bit, to avoid the possibility of rounding causing overflow
when computing 1/sfmin. */
sfmin = FLT_MIN;
small = one / FLT_MAX;
if (small >= sfmin) sfmin = small * (one + FLT_EPSILON);
return sfmin;
} else if ('U' == ch || 'u' == ch) {
return FLT_MIN;
}
return zero;
}
double
dlamch(char *cmach) {
char ch = cmach[0];
double sfmin, small, one = 1.0, zero = 0.0;
if ('B' == ch || 'b' == ch) {
return FLT_RADIX;
} else if ('E' == ch || 'e' == ch) {
return DBL_EPSILON;
} else if ('L' == ch || 'l' == ch) {
return DBL_MAX_EXP;
} else if ('M' == ch || 'm' == ch) {
return DBL_MIN_EXP;
} else if ('N' == ch || 'n' == ch) {
return DBL_DIGITS;
} else if ('O' == ch || 'o' == ch) {
return DBL_MAX;
} else if ('P' == ch || 'p' == ch) {
return DBL_EPSILON * FLT_RADIX;
} else if ('R' == ch || 'r' == ch) {
return FLT_ROUNDS < 2 ? one : zero;
} else if ('S' == ch || 's' == ch) {
/* Use SMALL plus a bit, to avoid the possibility of rounding causing overflow
when computing 1/sfmin. */
sfmin = DBL_MIN;
small = one / DBL_MAX;
if (small >= sfmin) sfmin = small * (one + DBL_EPSILON);
return small;
} else if ('U' == ch || 'u' == ch) {
return DBL_MIN;
}
return zero;
}
int
main() {
int i, j;
char cmach[] = "belmnprosu";
double v;
for (j = 0; j < 2; ++j)
for (i = 0; cmach[i]; ++i) {
if (j)
v = dlamch( cmach + i );
else
v = (double)slamch( cmach + i );
printf( "%c %23.15e\n", cmach[i], v );
}
return 0;
}
program lamch
! implicit none
integer i
external tlamch
do i = 1, 2
call tlamch('B', i)
call tlamch('E', i)
call tlamch('L', i)
call tlamch('M', i)
call tlamch('N', i)
call tlamch('P', i)
call tlamch('R', i)
call tlamch('O', i)
call tlamch('S', i)
call tlamch('U', i)
end do
end program lamch
subroutine tlamch(cmach, idx)
character, intent(in) :: cmach
integer, intent(in) :: idx
double precision x
real slamch
double precision dlamch
external dlamch, slamch
intrinsic dble
if (idx .eq. 1) then
x = dble(slamch(cmach))
write( *, '(a e25.15)') cmach, x
else
write( *, '(a e25.15)') cmach, dlamch(cmach)
end if
end subroutine tlamch
real function slamch(cmach)
character, intent(in) :: cmach
real zero, one
parameter (zero = 0.0e+0, one = 1.0e0)
real sfmin, small
intrinsic digits, epsilon, radix, real
if (cmach .eq. 'B' .or. cmach .eq. 'b') then
slamch = radix(zero)
else if (cmach .eq. 'E' .or. cmach .eq. 'e') then
slamch = epsilon(zero)
else if (cmach .eq. 'L' .or. cmach .eq. 'l') then
slamch = maxexponent(zero)
else if (cmach .eq. 'M' .or. cmach .eq. 'm') then
slamch = minexponent(zero)
else if (cmach .eq. 'N' .or. cmach .eq. 'n') then
slamch = real(digits(zero))
else if (cmach .eq. 'O' .or. cmach .eq. 'o') then
slamch = huge(zero)
else if (cmach .eq. 'P' .or. cmach .eq. 'p') then
slamch = epsilon(zero) * radix(zero)
else if (cmach .eq. 'R' .or. cmach .eq. 'r') then
slamch = one
else if (cmach .eq. 'S' .or. cmach .eq. 's') then
sfmin = tiny(zero)
small = one / huge(zero)
if (small .ge. sfmin) sfmin = small * (one+epsilon(zero))
slamch = sfmin
else if (cmach .eq. 'U' .or. cmach .eq. 'u') then
slamch = tiny(zero)
else
slamch = zero
end if
return
end function slamch
double precision function dlamch(cmach)
character, intent(in) :: cmach
double precision zero, one
parameter (zero = 0.0d+0, one = 1.0d0)
double precision :: small, sfmin
intrinsic digits, dble, epsilon, huge, maxexponent, minexponent
intrinsic radix, tiny
if (cmach .eq. 'B' .or. cmach .eq. 'b') then
dlamch = radix(zero)
else if (cmach .eq. 'E' .or. cmach .eq. 'e') then
dlamch = epsilon(zero)
else if (cmach .eq. 'L' .or. cmach .eq. 'l') then
dlamch = maxexponent(zero)
else if (cmach .eq. 'M' .or. cmach .eq. 'm') then
dlamch = minexponent(zero)
else if (cmach .eq. 'N' .or. cmach .eq. 'n') then
dlamch = dble(digits(zero))
else if (cmach .eq. 'O' .or. cmach .eq. 'o') then
dlamch = huge(zero)
else if (cmach .eq. 'P' .or. cmach .eq. 'p') then
dlamch = epsilon(zero) * radix(zero)
else if (cmach .eq. 'R' .or. cmach .eq. 'r') then
dlamch = one
else if (cmach .eq. 'S' .or. cmach .eq. 's') then
sfmin = tiny(zero)
small = one / huge(zero)
if (small .ge. sfmin) sfmin = small * (one+epsilon(zero))
dlamch = sfmin
else if (cmach .eq. 'U' .or. cmach .eq. 'u') then
dlamch = tiny(zero)
else
dlamch = zero
end if
return
end function dlamch
-------------------------------------------------------------------------
Take Surveys. Earn Cash. Influence the Future of IT
Join SourceForge.net's Techsay panel and you'll get the chance to share your
opinions on IT & business topics through brief surveys -- and earn cash
http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=DEVDEV
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion
