Bug in MKL-Lapack95 GESDD

[mecej4]

The routines GESVD and GESDD perform the same task (computing the SVD) but use different algorithms and have slightly different allowed values of JOB (for GESVD) and JOBZ (for GESDD). When both the U and VT arguments are present, the F95-F77 wrapper apparently errs in the case of GESDD and returns an error from the F77 routine.

The code (run on Suse 11.1, Intel64, using Intel l_cprof_p_11.1.072 compiler and bundled MKL 10.2), compile/link line included at top of source code):

[fortran]!ifort -I/opt/IFC11/mkl/include/em64t/lp64 -g gesddex.f90 
!-L$MKLPATH -lmkl_lapack95_lp64 $MKLPATH/libmkl_solver_lp64_sequential.a 
!-Wl,--start-group -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread
!$MKLPATH is set to /opt/IFC11/mkl/lib/em64t
program GESDDtst

 use mkl95_lapack, only : GESDD,GESVD
 use mkl95_precision, only : wp=>dp
 implicit none

 real (kind=wp) :: A(3,5), U(3,3), S(3), VT(3,5)
 integer :: i,j, INFO
 data A/-2d0,2d0,6d0,1d0,-2d0,-8d0,-5d0,0d0,0d0,11d0,1d0,6d0,2d0,4d0,0d0/

! call GESVD( A, S, U, VT, JOB='U', INFO=INFO)       !works correctly
 call GESDD( A, S, U, VT, JOBZ='A',INFO=INFO)       ! fails to return U, VT
! call GESDD( A, S,INFO=INFO)                                  ! returns S correctly
 write(*,5)INFO
 write(*,10)((U(i,j),j=1,3),i=1,3)
 write(*,*)
 write(*,10)S
 write(*,*)
 write(*,20)((VT(i,j),j=1,5),i=1,3)
  5 format(' Info = ',I4)
 10 format(1x,3F10.4)
 20 format(1x,5F10.4)
end program GESDDtst
[/fortran]

gives the incorrect result

[bash]MKL ERROR: Parameter 10 was incorrect on entry to DGESDD
 Info =  -10
     0.0000    0.0000    0.0000
     0.0000    0.0000    0.0000
     0.0000    0.0000    0.0000

     0.0000    0.0000    0.0000

     0.0000    0.0000    0.0000    0.0000    0.0000
     0.0000    0.0000    0.0000    0.0000    0.0000
     0.0000    0.0000    0.0000    0.0000    0.0000
[/bash]

The corresponding call to GESVD works correctly. If GESDD is called with (A,S,INFO=INFO), the singular values are returned correctly. These alternative calls are commented out in the source code.

In the MKL documentation for GESDD, under "Fortran 95 Interface Notes" the optional argument JOBZ is incorrectly listed as JOB.

[Chao_Y_Intel]

Hi,

Thanks. We will check the code. What is the hardward platfrom you tested?

Thanks,
Chao

[Michael_C_Intel4]

Hi,

parameter 10 incorrect means LDVT size is inappropriate. The MKL documentation requires first dimension of VT to be at least n, so that VT would keep n-by-n orthogonal matrix. Just declare VT(5,5) instead of VT(3,5) and it'll be OK:

Info = 0

-0.7409 0.6716 0.0002

-0.1824 -0.2009 -0.9625

-0.6463 -0.7132 0.2713

14.0828 10.1124 3.9260

-0.1960 0.3404 0.2631 -0.8671 -0.1570

-0.5957 0.6704 -0.3321 0.2875 0.0534

-0.0758 -0.0624 -0.0003 0.1701 -0.9805

Michael.

[mecej4]

Thanks for the quick reply. The message about parameter 10 being incorrect is from the F77 routine DGESDD, which my code did not call directly. I called the F95 routine GESDD, for which the MKL Ref. manual says:

vt Holds the matrix VT of size (min(m, n),n).

This is inconsistent with what you wrote as the needed size for VT.

Perhaps, the documentation needs revision. At any rate, the wrapper routine should check the parameters before calling the F77 routine. and any message/error INFO would be better from the wrapper routine than the underlying F77 routine. In this case, the wrapper routine has all the information needed to return an error message without calling the F77 routine.

The F95 wrapper routine GESDD allocates the workspace arryas needed by DGESDD, and translates the F95 call with optional arguments to an F77 call with a fixed number of arguments. Typical usage of GESDD in F95 code should not require decoding the argument specifications of the F77 routine -- after all, the F95 routine is the middle-man.

[Michael_C_Intel4]

Hi,

in the current MKL manual ldvt is described as the following:

"If jobz = 'A' or jobz = 'O' and m < n, then ldvt n;

If jobz = 'S', ldvt min(m, n)."

So if you call DGESDD with jobz='A' ldvt should be larger or equal to n (that is, 5). If your documentation says something different it's a bug,and it has been fixed already.

I agree with you that the number of the wrong parameter associated with different interface looks strange, but there'sthe rationale behind it. Suppose, you get a message about wrong array parameter concerned F95 interface, then you hardly can guess that it's about a leading dimension. Referring to F77 error messages would help a lot, as in your testcase. I don't think it's going to be changed.

Probably we shoulddescribe the possible error messagemore clearly, for now we only have the statement about INFO parameter, whichvalue refers to F77 interface (at Fortran 95 Interface Conventions):

"Argument info is declared as optional in the Fortran 95 interface. If it is present in the calling sequence, the value assigned to info is interpreted as follows:

If this value is more than -1000, its meaning is the same as in the FORTRAN 77 routine."

I think there should be a guidance how to interprete the error message in case of F95 interface. Is it OK for you?

[mecej4]

Thanks for the reply.

I understand that it is troublesome to write wrapper routines in Fortran 95/2003 to interface to Fortran 77 computational routines.

However, once an F95 interface is provided and documented, a user should not have to be bothered with F77-level details such as leading-array-dimension, number of arguments, etc. The MKL F95 wrapper routine can extract all such information by applying the F95 intrinsic array inquiry functions PRESENT, SHAPE, SIZE, LBOUND and UBOUND to the arguments, before calling the F77 routine.

The F95 wrapper should take responsibility for (i) checking that the F95 arguments are correct, (ii) using the F95 arguments, produce a corresponding (usually longer) list and pass that pass to the F77 routine, and (iii) convert any error codes and restate messages given by the F77 routines, in the context of the F95 documentation. What you wrote amounts to doing only (i) and (ii).

Please consider the following. Assume that I do not know or care about the F77 interface to MKL. In the following example, which is a simplified version of my original example, I have only F95 code and F95 MKL subroutines called. The comments at the end of the program are from the current MKL documentation. Let me emphasize that I have no criticism about the F77 interface, and do not claim to have seen errors in the documentation of the F77 calls to ?GESVD. Nor is the issue a 'show-stopper' for me.

In the code, the array sizes declared in lines 6-7 satisfy the argument rules as stated in the current MKL documentation, included in lines 32-40 of the code for your convenience and taken from the URL given on line-23.

Yet, this short example code fails, and only gives the F77 error message. Based on what has been written in this thread, I think that the rule on line-38 is incorrect for the case of JOBZ='A' (the documentation incorrectly lists JOB in place of JOBZ).

The error goes away if line-7 is replaced by

real (kind=wp) :: A(M,N), U(M,MIN(M,N)), S(MIN(M,N)), VT(N,N)

but this contains a declaration of VT that is at variance with the documentation.

[fortran]program GESDDtst 
 
 use mkl95_lapack, only : GESDD
 use mkl95_precision, only : wp=>dp 
 implicit none 
 integer, parameter :: M=3, N= 5 
 real (kind=wp) :: A(M,N), U(M,MIN(M,N)), S(MIN(M,N)), VT(MIN(M,N),N) 
 integer :: i,j, INFO 
 data A/-2d0,2d0,6d0,1d0,-2d0,-8d0,-5d0,0d0,0d0,11d0,1d0,6d0,2d0,4d0,0d0/ 
 
 call GESDD( A, S, U, VT, JOBZ='A',INFO=INFO)       ! fails to return correct U, VT 
 write(*,5)INFO 
 write(*,10)((U(i,j),j=1,3),i=1,3) 
 write(*,*) 
 write(*,10)S 
 write(*,*) 
 write(*,20)((VT(i,j),j=1,5),i=1,3) 
  5 format(' Info = ',I4) 
 10 format(1x,3F10.4) 
 20 format(1x,5F10.4) 
end program GESDDtst 
 
!http://software.intel.com/sites/products/documentation/hpc/compilerpro/en-us/cpp/win/mkl/refman/index.htm

!Fortran 95 Interface Notes

!Routines in Fortran 95 interface have fewer arguments in the calling sequence than 
!their FORTRAN 77 counterparts. For general conventions applied to skip redundant or restorable arguments, see Fortran 95 Interface Conventions.

!Specific details for the routine gesdd interface are the following:

!a    Holds the matrix A of size (m, n).

!s    Holds the vector of length min(m, n).

!u    Holds the matrix U of size (m,min(m, n)).

!vt   Holds the matrix VT of size (min(m, n),n).

!job  Must be 'N', 'A', 'S', or 'O'. The default value is 'N'. <<=== error, 'jobz' is the name, not 'job'
[/fortran]

[Michael_C_Intel4]

Thanks for pointing at the documentation bugs. We will surely fix them as soon as possible.

I agree with you again that the error reports in F95 should be different from F77.We will address the issue, but can't change it quickly, sorry for inconvenience.

[Gennady_F_Intel]

This issue has been submitted to our internal development tracking database for further investigation, we will inform you once a new update becomes available.Here is a bug tracking number for your reference: DPD200192952.

--Gennady

[mecej4]

Michael, Gennady:

Thanks to both of you for recognizing the documentation/F95 interface issues.

Is it possible for me to check the status using the DPD number, without bothering readers of the forum by posting a message with just a request for status?