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
