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