Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
7036 Discussions

## How to find inverse of a sparse matrix using pardiso?

Beginner
2,642 Views

I have a sparse matrix of A 17000 x 17000 (real data). Need to find the inverse of A, I am new to intel math library.

Anyone could help me to solve this problem by using pardiso?, pardiso has take 16 variables but I have only a single matrix and if use pardiso where will get the corresponding result?

My system configuration are as follows,

Architecture : x86_64

CPUs : 12

Model : Intel(R) Core(TM) i7-3960x CPU @ 3.30Ghz

RAM : 16 Gb

ISS: Intel system studio 2019.

1 Solution
Employee
2,642 Views

Hi,

I can't help but wonder why it is you are wanting the inverse of the sparse matrix?  This is a very ill conditioned problem. If it truly is needed, the best way to do it is by assembling each column or set of columns by solving Ax = e_i   where e_i is the elemental vector with zeros everywhere but a 1 at the ith spot.  You can do this in batches of rhs, but it may be required to use a number of passes to get them all.  Direct sparse solvers like pardiso try to create a Cholesky (L L^T)  or LU factorization of the matrix A which can then be used to solve any system. (Note: so do direct dense solvers like in LAPACK)   (Note: of course to maintain stability there are a lot of scaling and permutation matices involved in the final solution so that it could look more like PAP^T = LDL^T  for Cholesky or many other variants of LU for symmetric, non symmetric, structurally symmetric matrices ... See https://software.intel.com/en-us/mkl-developer-reference-c-intel-mkl-pardiso-parallel-direct-sparse-solver-interface  for some details and references for more details...  These are all used in the solve step to produce the final solution of your system...)

If you are going to solve for the inverse with pardiso, you can perform a single time the phase 12 (start with phase 1 = symbolic factorization, end on phase 2 = numerical factorization) to obtain the desired factorizations.  Then you can run phase 33 (start and end on 3 = solve system)  many times with the various rhs or sets of rhs reusing the matrix factorization.  I would start with 8 or 16 rhs at a time, and then play around with various powers of 2 to see what gives the best performance... It doesn't have to be powers of 2 but those will be best performance :)  It will be solved for in dense 17k x (8 or 16 ...) chunks and sparse matrices become dense inverse matrices so this is probably the best format for storing the inverse...  Copy each solution as you solve it into the right place in a separate 17kx17k memory chunk and eventually you will get your full inverse.

Now, the real question is why you would ever actually need the inverse of a matrix.  I am sure there are reasons and I would be interested to hear what your reasonings are, but it is rare to find a problem that is best solved by first getting the inverse and then doing something with it.  It is way more efficient and numerically stable to use the LU factorizations to do your operations in place of the inverse...

Anyway,  I hope this helps a little and I look forward to finding out more.

As a side note, if you are not using CSR format, then what are you using?  Pardiso is a for sparse matrices, and CSR or BSR matrix formats are the standards for that domain.  If you have a sparse matrix but are storing it in dense format then as was mentioned above, the getri type LAPACK functions are appropriate.  There are both C and Fortran interfaces for that functionality.  It was originally written in fortran so fortran interfaces are simpler for LAPACK, but not by much.  Here is the link for the C interface https://software.intel.com/en-us/mkl-developer-reference-c-getri. The C interface adds a LAPACKE_ to the front... You should also check out the lapack function finding advisor https://software.intel.com/en-us/articles/intel-mkl-function-finding-advisor to select the right function for your matrix and operation...

Best Regards,

Spencer

7 Replies
Moderator
2,642 Views

17K is too small case. do you really need to use Pardiso in that case? Are your data is represent into CSR format? if not and your input data representation in the form of dense matrix then you may try to look at getri function from Lapack.

Moderator
2,642 Views

https://software.intel.com/en-us/mkl-developer-reference-fortran-getri

and there are some example which you may find into standard example directory

Beginner
2,642 Views

@ Gennady F, Thanks for your response.I am doing 'C' code,but getri is fortran. How can i use getri in my C code?

I am not using CSR format, If I like to proceed with pardiso, How can I proceed?

Employee
2,643 Views

Hi,

I can't help but wonder why it is you are wanting the inverse of the sparse matrix?  This is a very ill conditioned problem. If it truly is needed, the best way to do it is by assembling each column or set of columns by solving Ax = e_i   where e_i is the elemental vector with zeros everywhere but a 1 at the ith spot.  You can do this in batches of rhs, but it may be required to use a number of passes to get them all.  Direct sparse solvers like pardiso try to create a Cholesky (L L^T)  or LU factorization of the matrix A which can then be used to solve any system. (Note: so do direct dense solvers like in LAPACK)   (Note: of course to maintain stability there are a lot of scaling and permutation matices involved in the final solution so that it could look more like PAP^T = LDL^T  for Cholesky or many other variants of LU for symmetric, non symmetric, structurally symmetric matrices ... See https://software.intel.com/en-us/mkl-developer-reference-c-intel-mkl-pardiso-parallel-direct-sparse-solver-interface  for some details and references for more details...  These are all used in the solve step to produce the final solution of your system...)

If you are going to solve for the inverse with pardiso, you can perform a single time the phase 12 (start with phase 1 = symbolic factorization, end on phase 2 = numerical factorization) to obtain the desired factorizations.  Then you can run phase 33 (start and end on 3 = solve system)  many times with the various rhs or sets of rhs reusing the matrix factorization.  I would start with 8 or 16 rhs at a time, and then play around with various powers of 2 to see what gives the best performance... It doesn't have to be powers of 2 but those will be best performance :)  It will be solved for in dense 17k x (8 or 16 ...) chunks and sparse matrices become dense inverse matrices so this is probably the best format for storing the inverse...  Copy each solution as you solve it into the right place in a separate 17kx17k memory chunk and eventually you will get your full inverse.

Now, the real question is why you would ever actually need the inverse of a matrix.  I am sure there are reasons and I would be interested to hear what your reasonings are, but it is rare to find a problem that is best solved by first getting the inverse and then doing something with it.  It is way more efficient and numerically stable to use the LU factorizations to do your operations in place of the inverse...

Anyway,  I hope this helps a little and I look forward to finding out more.

As a side note, if you are not using CSR format, then what are you using?  Pardiso is a for sparse matrices, and CSR or BSR matrix formats are the standards for that domain.  If you have a sparse matrix but are storing it in dense format then as was mentioned above, the getri type LAPACK functions are appropriate.  There are both C and Fortran interfaces for that functionality.  It was originally written in fortran so fortran interfaces are simpler for LAPACK, but not by much.  Here is the link for the C interface https://software.intel.com/en-us/mkl-developer-reference-c-getri. The C interface adds a LAPACKE_ to the front... You should also check out the lapack function finding advisor https://software.intel.com/en-us/articles/intel-mkl-function-finding-advisor to select the right function for your matrix and operation...

Best Regards,

Spencer

Valued Contributor III
2,642 Views

Hello:

There are examples on this forum of the code to turn a matrix in your form into the form needed for the PARDISO solver, mecej4 provided the much needed clues and code when I had this problem.  You will need to hunt a bit but the code is present and it does work.  Once you code it, it is very easy to use.

Honored Contributor III
2,642 Views

V, chandramohan wrote:
I am doing 'C' code,but getri is fortran. How can i use getri in my C code?

Almost all MKL routines have C and Fortran interfaces, and reading the MKL documentation will help to avoid forming incorrect opinions. See https://software.intel.com/en-us/mkl-developer-reference-c-getri .

V, chandramohan wrote:
I am not using CSR format, If I like to proceed with pardiso, How can I proceed?

You cannot proceed until you put your data in CSR format. Failing that, look elsewhere for a library that accepts whatever format you insist on using. Conversion routines, either in MKL or user-provided, can be used as needed to put the matrix data into CSR format before calling Pardiso.

Spencer P. gave you an excellent and complete response. Consider reading it again, and do so carefully.

Your questions and statements indicate that you need to do considerable reading before attempting to use Pardiso in particular and MKL in general.

Beginner
2,642 Views

Spencer P. (Intel) wrote:

Hi,

I can't help but wonder why it is you are wanting the inverse of the sparse matrix?  This is a very ill conditioned problem. If it truly is needed, the best way to do it is by assembling each column or set of columns by solving Ax = e_i   where e_i is the elemental vector with zeros everywhere but a 1 at the ith spot.  You can do this in batches of rhs, but it may be required to use a number of passes to get them all.  Direct sparse solvers like pardiso try to create a Cholesky (L L^T)  or LU factorization of the matrix A which can then be used to solve any system. (Note: so do direct dense solvers like in LAPACK)   (Note: of course to maintain stability there are a lot of scaling and permutation matices involved in the final solution so that it could look more like PAP^T = LDL^T  for Cholesky or many other variants of LU for symmetric, non symmetric, structurally symmetric matrices ... See https://software.intel.com/en-us/mkl-developer-reference-c-intel-mkl-par...  for some details and references for more details...  These are all used in the solve step to produce the final solution of your system...)

If you are going to solve for the inverse with pardiso, you can perform a single time the phase 12 (start with phase 1 = symbolic factorization, end on phase 2 = numerical factorization) to obtain the desired factorizations.  Then you can run phase 33 (start and end on 3 = solve system)  many times with the various rhs or sets of rhs reusing the matrix factorization.  I would start with 8 or 16 rhs at a time, and then play around with various powers of 2 to see what gives the best performance... It doesn't have to be powers of 2 but those will be best performance :)  It will be solved for in dense 17k x (8 or 16 ...) chunks and sparse matrices become dense inverse matrices so this is probably the best format for storing the inverse...  Copy each solution as you solve it into the right place in a separate 17kx17k memory chunk and eventually you will get your full inverse.

Now, the real question is why you would ever actually need the inverse of a matrix.  I am sure there are reasons and I would be interested to hear what your reasonings are, but it is rare to find a problem that is best solved by first getting the inverse and then doing something with it.  It is way more efficient and numerically stable to use the LU factorizations to do your operations in place of the inverse...

Anyway,  I hope this helps a little and I look forward to finding out more.

As a side note, if you are not using CSR format, then what are you using?  Pardiso is a for sparse matrices, and CSR or BSR matrix formats are the standards for that domain.  If you have a sparse matrix but are storing it in dense format then as was mentioned above, the getri type LAPACK functions are appropriate.  There are both C and Fortran interfaces for that functionality.  It was originally written in fortran so fortran interfaces are simpler for LAPACK, but not by much.  Here is the link for the C interface https://software.intel.com/en-us/mkl-developer-reference-c-getri. The C interface adds a LAPACKE_ to the front... You should also check out the lapack function finding advisor https://software.intel.com/en-us/articles/intel-mkl-function-finding-adv... to select the right function for your matrix and operation...

Best Regards,

Spencer