Page 1 of 1

Division by 0 in *GESVJ & *GEVSJ with zero columns

PostPosted: Mon Aug 16, 2010 11:51 am
by bartoldeman
Hi,

the Jacobi SVD routines, when floating point trapping is enabled crash if a column in the input matrix consists of just zeros. I am not sure if these routines are supposed to work in an environment where this happens -- from what I have read that depends on an ILAENV 10/11 query which *GESVJ/*GEVSJ don't use.

A trivial example is to just pass the zero 1x1 matrix as follows:

Code: Select all
  double precision a(1,1), s(1), u(1), v(1,1), iwork(8), work(100)
  integer info
  a(1,1) = 0
  call dgesvj('G','N','V',1,1,a,1,s,1,v,1,work,100,info)


compiling with
Code: Select all
gfortran -ffpe-trap=zero -llapack svd.f90

a floating point exception happens when run.

I attached a patch:
the code calls *LASSQ on every column with parameters so that
AAQQ will contain the L2-norm of the column divided by the max norm, or 0 if the column contains all zeros and
AAPP will contain the max norm.

Normally, AAQQ.GE.1, so BIG/AAQQ is well defined.
Only in the case of an all zeros column a division by zero occurs.

Regards,
Bart Oldeman

Re: Division by 0 in *GESVJ & *GEVSJ with zero columns

PostPosted: Wed Aug 25, 2010 7:36 pm
by Julien Langou
Dear Bart Oldeman,
Thanks a lot for the bug report and the patch, Zlatko just provided us with corrected xGESVJ and xGEVSJ, the fixed routines are currently available in the svn repository.
(Bug report listed as bug #0063 in lapack Errata http://www.netlib.org/lapack/Errata/.)
Zlatko fix is not exactly the same as yours but I believe it should work just as fine.
Cheers,
Julien