AnsweredAssumed Answered

Different result from zgesdd() with ifort and gfortran

Question asked by lewisrd on Jun 23, 2014

I wish to compute the singular value decomposition (SVD) of a particular complex-valued 57-by-384 matrix. When I attempt this computation using zgesdd(), the calculation succeeds when I compile with gfortran and link against the gfortran64 version of ACML 5.3.1.  However, the routine fails with error code info=1 when I compile with ifort and link against the ifort64 version of ACML 5.3.1.  This is obviously unexpected and very distressing!  According to the man page of zgesdd(), the error code of 1 means that "[t]he updating process of DBDSDC did not converge."

 

The matrix that causes the problem is attached. It is stored as column major, little endian, double precision complex numbers, with real and imaginary parts interleaved.  As far as I can tell, there is nothing unusual about this matrix that should cause problems during the computation of its SVD.

 

Here is Fortran code that exhibits the problem:

 

PROGRAM tzsvd
  IMPLICIT NONE

  INTEGER*4, PARAMETER :: k=57, m=384, lwork=147456

  COMPLEX*16 :: B(k,m), U(k,k), VT(k,m), work(lwork)
  REAL*8     :: sval(k), &
                rwork(MIN(m,k)*MAX(5*MIN(m,k)+7,2*MAX(m,k)+2*MIN(m,k)+1))
  INTEGER*4  :: info, iwork(8*MIN(m,k))

 

  ! Load the matrix from disk

  CALL rzmtx(B, k, m)

 

  ! Compute the SVD

  CALL zgesdd('S', k, m, B, k, sval, U, k, VT, k, work, lwork, &
       rwork, iwork, info)
  IF (info /= 0) THEN
     PRINT '("zgesdd() failed, info=",I5)', info
     STOP
  END IF

END PROGRAM tzsvd


The routine rzmtx() simply reads the matrix from disk. Here is its C implementation (modify the hard-coded path as appropriate):


#include <stdio.h>
#include <stdint.h>

void rzmtx_(void *d, int32_t *r, int32_t *c)
{
    FILE *f;
    f = fopen("/path/to/mtx-57-384-zgesdd-fail.bin", "r");
    fread(d, 2*sizeof(double), (*r)*(*c), f);
    fclose(f);
    return;
}


Perhaps there is something wrong with my call to zgesdd(), but if so, I have been unable to spot my error. I would be very grateful for any assistance.

 

Here are a few additional notes:

OS: Linux

Architecture: x86-64

gfortran version: 4.8.2

ifort version: 14.0.3

I have verified that this behavior is the same for ACML 5.3.0 and 5.3.1.

Neither Intel's MKL nor Netlib's LAPACK not have any trouble with the call to zgesdd(). The issue appears to be unique to ACML.

 

Thank you in advance for your help,

--Ryan

 

UPDATE (20140624): I mistakenly attached the wrong data file to my original message, and have now corrected this error.  The sha1sum of the correct attachment is 56236c3135f9c34a22a0e1a15e8e60f0b17fc185.

Outcomes