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

Open discussion regarding features, bugs, issues, vendors, etc.

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

Postby bartoldeman » Mon Aug 16, 2010 11:51 am

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
Attachments
gesvj_gejsv.zip
Contains patch and example
(1.03 KiB) Downloaded 25 times
bartoldeman
 
Posts: 4
Joined: Mon Mar 12, 2007 5:29 pm

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

Postby Julien Langou » Wed Aug 25, 2010 7:36 pm

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
Julien Langou
 
Posts: 735
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA


Return to User Discussion

Who is online

Users browsing this forum: Yahoo [Bot] and 2 guests