[Numpy-discussion] Fwd: Collection.findTransformation() never stops

Piotr Luszczek S luszczek at cs.utk.edu
Tue Sep 19 09:21:35 CDT 2006


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 <oliphant.travis at ieee.org> 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 at lists.sourceforge.net
> > https://lists.sourceforge.net/lists/listinfo/numpy-discussion
-------------- next part --------------
A non-text attachment was scrubbed...
Name: lamch.c
Type: text/x-csrc
Size: 2371 bytes
Desc: not available
Url : http://projects.scipy.org/pipermail/numpy-discussion/attachments/20060919/5d0d9d48/attachment-0001.bin 
-------------- next part --------------
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


More information about the Numpy-discussion mailing list