2 Replies Latest reply on Jun 24, 2014 10:42 AM by lewisrd

    Different result from zgesdd() with ifort and gfortran

    lewisrd

      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.