How to use MKL Routines for Matrix Inversion

[Jeferson_Vanderlinde]

Hi everyone!

I never used libraries in my programs, usually I do all my routines on my own, but now I'm having efficiency problems and I decided to use libraries to improve the efficiency of my program.

I need an efficient routine to invert a sparse matrix. I implemented the Gauss-Jordan elimination algorithm using the "do concurrent" statement to make it parallelizable.

The problem I need to solve has a sparse rectangular matrix with m rows and n columns where n is a large number of columns. I use on each iteration of my program a sparse square matrix with dimension m x m formed by the inverse matrix of m columns of the sparse rectangular matrix mentioned above. In each iteration a new column of the rectangular matrix is selected to be part of the square matrix and my program does all the necessary calculations to enter this new column.

The problem appears after several iterations where appears several numerical errors in the square matrix due to the iterative process. One way to clean these errors is to calculate the inverse matrix of all columns of the rectangular matrix that are being used in the square matrix.

The university computer that I am using has installed Microsoft Visual Studio Professional 2013 Update 3, Intel Parallel Studio XE 2016 Update 2 Cluster Edition for Windows with Intel Math Kernel Library 11.3.

I found the routine P?GETRI from ScaLAPACK in the MKL documentation* that can solve my problem, but I did not find examples of how to use this routine and I'm in trouble of how to use the routine and declare variables. I found only one example of routine DGETRI from LAPACK* that has helped me to understand some of how to use the routine P?GETRI.

Does anyone know another routine that solves my problem or give me an example of how to use the routine P?GETRI?

I activated the MKL library on project properties of my of project in Visual Studio following the steps of this video tutorial https://www.youtube.com/watch?v=E1n8-azFvmg. Is the procedure on the video correct or do I need to perform other procedures to use the MKL library on Visual Studio?

*I add the links in the text

[Ying_H_Intel]

Hi Jeferson, 

Only the function P?GETRI can't solve linear equation, you may consider several factors and try them out: 

1) Cluster vs. Single 

MKL provide Cluster lapack and lapack API, for example, Pdgesv or dgesv.  The "P" is for cluster distributed version, which require MPI cluster.  In single machine, dgesv was used often. So if you really need the scalapack, you may need install MPI on your windows machines cluster. 

2) API choose for dense matrix

 MKL provide several functions for solve linear equation. For dense matrix, like you mentioned P?GETRI ,  but as documentation

Routines for Matrix Inversion
This sections describes ScaLAPACK routines that compute the inverse of a matrix based on the previously
obtained factorization. Note that it is not recommended to solve a system of equations Ax = b by first
computing A-1 and then forming the matrix-vector product x = A-1b. Call a solver routine instead (see
Solving Systems of Linear Equations); this is more efficient and more accurate.

Like pdgesv,  you may find some scalapack sample in MKL install directory for example,

C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2017.0.072\windows\mkl\examples

or  in the forum, for example:  like https://software.intel.com/en-us/forums/intel-math-kernel-library/topic/558359

3) for sparse matrix,  MKL have pardiso solver and cluster sparse solver
you can find the c or fortran sample under  MKL install folder/example  

4) Regarding the using MKL in MSVS2013 IDE.  You may refer to the online article. 

https://software.intel.com/en-us/articles/intel-math-kernel-library-intel-mkl-compiling-and-linking-with-microsoft-visual-cc/

or you can use the way of command line build as showed in cluster sample in  C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2016.xxx \windows\mkl\examples.  (if you build MPI code in  MSVC IDT, you may also refer to Intel MPI documentation) 

https://software.intel.com/en-us/articles/intel-math-kernel-library-parallel-direct-sparse-solver-for-clusters

Best Regards,

Ying