topic Re: Calling mkl zgesvd from C in Intel® oneAPI Math Kernel Library

Calling mkl zgesvd from C

geerits — Wed, 15 Dec 2021 10:59:04 GMT

Dear reader,

So far I have been using the IMSL Fortran library for all of my computational work. My main code is in C and so I'm calling Fortran library from C. Currently, I'm testing to replace the double precision Singular Value Decomposition (SVD) of the IMSL Fortran library (DLSVCR()) by the MKL one (zgesvd()). I compile and link with the intel compiler (version 19) using the Visual Studio IDE. I get no errors but the result of the SVD is wrong. I have a strong impression this has something to do with the data type mapping. Below the call, as far as I think important to the problem:

######C-code snapshot#####################################

char jobu, jobvt;

jobu = 'A';
jobvt = 'A';

info = -1; // info should be zero after the call. Indeed it is!

zgesvd(&jobu, &jobvt, &NRA, &NCA, SCAT_TR, &LDA, S_SVD_2, U_SVD, &LDU, V_SVD, &LDV, work, lwork, rwork, &info);

####END of snapshot####################################################

As far as the Fortran double complex data type I have used my own structure, which is included in my header file (snapshot below):

typedef struct
{
double real;
double imag;
} COMPLEX16;

The prototype for zgesvd() is listed in my header file as follows:

extern "C"

{

void zgesvd(char *jobu, char *jobvt, int *NRA, int *NCA, COMPLEX16 *SCAT_TR, int *LDA, double *S_SVD, COMPLEX16 *U_SVD, int *LDU, COMPLEX16 *V_SVD, int *LDV, COMPLEX16 *work, int *lwork, double *rwork, int *info);

}

PS: With the IMSL routine all of this works just fine.

Any comments what could be the problem with my approach in calling zgesvd()??

With kind regards,

Tim.

Re: Calling mkl zgesvd from C

mecej4 — Wed, 15 Dec 2021 13:58:59 GMT

If you wish to call ZGESVD from C through the F77 interface, see the example (part of the MKL examples provided in Zip form in the .../mkl/examples directory, or zgesvd_ex.c.html (intel.com) .

Even if you get the calling sequence right, there can still be errors caused by not linking with the appropriate MKL libraries. The MKL Link Line Advisor can be of help.

Re: Calling mkl zgesvd from C

geerits — Wed, 15 Dec 2021 14:38:39 GMT

Many thanks for your prompt and useful reply! I double checked the calling sequence and the user-defined type. All seems to be correct. You refer to the MKL link-line advisor. However, I'm using MS VS IDE and do not know WHERE to invoke the link line (Linker>Command Line ???). The MKL documentation appears to be very vague on how to use the MKL libraries in conjunction with the VS IDE. The only thing I could find was changing the setting under "Use Intel MKL" under "Intel Performance Libraries" under the property pages. Looking forward to your insights.

Re: Calling mkl zgesvd from C

Gennady_F_Intel — Wed, 15 Dec 2021 14:53:29 GMT

You may try to check the developer guide following the link: https://www.intel.com/content/www/us/en/develop/documentation/onemkl-windows-developer-guide/top/programming-onemkl-in-ide/configure-integrated-dev-env-to-link-with-onemkl/config-microsoft-visual-c-dev-sys-to-link-onemkl.html

Re:Calling mkl zgesvd from C

Gennady_F_Intel — Wed, 15 Dec 2021 15:01:52 GMT

and also You might build the F77 example which meaej4 pointed out, from command line prompt and see if the computed results could be correct or not.

-Gennady

Re: Calling mkl zgesvd from C

geerits — Wed, 15 Dec 2021 15:32:16 GMT

Yes, I know that link. But that only mentions what I already mentioned to mecej4, i.e., changing the setting under "Use Intel MKL" under "Intel Performance Libraries" under the property pages. The problem is that you do not know what logic the IDE hides in the background. Yes, if I do not get this problem solved quickly, I might have to dive into the command line approach again (like in the old UNIX days :-)) as suggested by you.

Re: Calling mkl zgesvd from C

mecej4 — Wed, 15 Dec 2021 18:16:14 GMT

I took the MKL ZGESVD example at https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/zgesvd_ex.c.htm , and modified the source file to use the 6 X 3 matrix that is shown in the IMSL documentation for LSVCR. I compiled and linked the C source file using the command

icl /Qmkl xzgesvd.c

and ran the resulting EXE. The resulting singular values and right and left singular vectors agreed with those that are displayed in the IMSL manual.

Re: Calling mkl zgesvd from C

Gennady_F_Intel — Thu, 16 Dec 2021 06:57:58 GMT

Tim,

>>Yes, I know that link. But that only mentions what I already mentioned to mecej4, i.e., changing the setting under "Use Intel MKL" under "Intel Performance Libraries" under the property pages. The problem is that you do not know what logic the IDE hides in the background.

<< Then You choose "Use Intel MKL" button, You will have the following options: No, Parallel, Sequential, Cluster. This option will say compiler to link against the corresponding mkl's libraries

parallel: mkl_intel_lp64.lib mkl_intel_thread.lib mkl_core.lib libiomp5md.lib

sequential: mkl_intel_lp64.lib mkl_sequential.lib mkl_core.lib

here you could see some details:

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-windows-developer-guide/top/linking-your-application-with-onemkl/linking-quick-start/using-the-qmkl-compiler-option.html

Re: Calling mkl zgesvd from C

Gennady_F_Intel — Thu, 16 Dec 2021 07:01:07 GMT

thanks mecej4,

I see the same outputs with the current version of oneMKL 2021. I switched on the verbose mode to see all needed the mkl and zgesvd calls details:

ZGESVD Example Program Results
MKL_VERBOSE oneMKL 2021.0 Update 4 Product build 20210904 for Intel(R) 64 architecture Intel(R) Advanced Vector Extensions 2 (Intel(R) AVX2) enabled processors, Win 2.60GHz cdecl intel_thread
MKL_VERBOSE ZGESVD(A,A,3,6,00000078D4F4F2C0,3,00000078D4F4F3E0,00000078D4F4F400,3,00000078D4F4F510,6,00000078D4F4F2A8,-1,00000078D4F4F490,0) 1.36ms CNR:OFF Dyn:1 FastMM:1 TID:0 NThr:2
MKL_VERBOSE ZGESVD(A,A,3,6,00000078D4F4F2C0,3,00000078D4F4F3E0,00000078D4F4F400,3,00000078D4F4F510,6,000001B4ECC88ED0,12,00000078D4F4F490,0) 278.30us CNR:OFF Dyn:1 FastMM:1 TID:0 NThr:2

Singular values
11.77 9.30 4.99

Left singular vectors (stored columnwise)
( -0.66, 0.00) ( 0.27, 0.00) ( -0.70, 0.00)
( -0.74, 0.04) ( -0.38, -0.07) ( 0.55, -0.06)
( -0.05, -0.13) ( -0.17, 0.86) ( -0.02, 0.45)

Right singular vectors (stored rowwise)
( -0.20, -0.22) ( -0.34, 0.35) ( -0.15, -0.23) ( -0.30, 0.08) ( -0.23, 0.60) ( -0.29, 0.03)
( -0.50, -0.02) ( 0.29, -0.02) ( 0.54, -0.14) ( -0.22, -0.27) ( 0.13, -0.14) ( -0.44, 0.04)
( -0.20, -0.10) ( 0.12, -0.23) ( -0.44, -0.44) ( -0.05, -0.09) ( 0.32, -0.01) ( 0.05, -0.62)

Re: Calling mkl zgesvd from C

geerits — Thu, 16 Dec 2021 10:43:47 GMT

Hi Gennady_F_Intel,

I use the "sequential" option. But that doesn't solve my problem. However, your link triggered a question from my side with regards to "cdecl" (IA-32) versus "LP64". Since my target platform in Visual Studio is X64, should I use different data types in my C call to Fortran LAPACK routine zgesvd (i.e., in particular, replace int* by MKL_INT*)??? Weird though that I do not get any error message during link time. Your comments/suggestions are very much appreciated!

Re: Calling mkl zgesvd from C

mecej4 — Thu, 16 Dec 2021 12:08:33 GMT

If you are building 64-bit EXEs or DLLs, you may choose either the LP64 or the ILP64 libraries. Unless you are going to use very large arrays (large enough to need 64-bit subscript variables), the LP64 model should suffice. In short, the difference is that integer arguments passed to MKL routines are 32-bit integers in LP64 and are 64-bit integers in ILP64.

There is no reason to tie IA32 and CDECL.

If you link to LP64 libraries when you should have used ILP64, or vice versa, you will not see any link errors since the routine names do not change from LP64 to ILP64. The program, when run, will either crash or produce incorrect results.

Re: Calling mkl zgesvd from C

geerits — Thu, 16 Dec 2021 14:56:21 GMT

Hi mecej4,

I found the problem after a second inspection of the example you sent me in your very first reply. It turned out I didn't understand the documentation when it comes to the use of the "work" array in zgesvd(). Its original size "1" (while lwork=-1) cannot be kept at successive calls. I now understand what is meant with "If lwork = -1, then a workspace query is assumed; the routine only calculates the optimal size of the work array, returns this value as the first entry of the work array, and no error message related to lwork is issued by xerbla. See Application Notes for details.", in the reference manual. By the way, the idea of having to make a separate call to zgesvd() in order to determine the optimum size of "work" seems unnecessary. What is the exact function of the "work" array?

In any case, I have clicked "Accept solution" in your original reply. Many thanks for your assistance!

With kind regards,

Tim.

Re: Calling mkl zgesvd from C

mecej4 — Thu, 16 Dec 2021 15:52:08 GMT

Major parts of BLAS and Lapack were developed in the 1970s and 1980s, and the Fortran versions available were various versions of Fortran 77 and earlier. Those versions of Fortran had rather primitive and highly restricted ways of allocating memory for workspace (a.k.a. scratchpad) memory. In particular, one could not read problem data, calculate the required workspace size, and then allocate the exact size required. Thus, a user had to estimate the size needed, write the program with that estimate and make trial runs. If the Lapack library found the size inadequate, it would print an error message and abort, following which the programmer would adjust the memory size and try again.

We do not need those restrictions anymore, but the interfaces have been set and so widely used in software that we still use the old protocols for assigning workspace, etc.

Re:Calling mkl zgesvd from C

Gennady_F_Intel — Fri, 17 Dec 2021 03:54:34 GMT

This issue has been resolved and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only.