Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- How to find inverse of a sparse matrix using pardiso?

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Chandramohan_V

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-28-2018
04:52 AM

435 Views

How to find inverse of a sparse matrix using pardiso?

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.

Thanks in advance.

1 Solution

Spencer_P_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-28-2018
10:36 AM

435 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

Link Copied

7 Replies

Gennady_F_Intel

Moderator

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-28-2018
07:31 AM

435 Views

Gennady_F_Intel

Moderator

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-28-2018
07:32 AM

435 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

Chandramohan_V

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-28-2018
08:54 AM

435 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?

Spencer_P_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-28-2018
10:36 AM

436 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

JohnNichols

Valued Contributor I

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-30-2018
07:07 AM

435 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.

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-01-2018
04:47 AM

435 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.

Chandramohan_V

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-02-2018
09:12 PM

435 Views

Spencer P. (Intel) wrote:Thank you for your respose.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

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.