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

## sparse Ax=b with different right sides

New Contributor I
549 Views

Hi,

A have an iteration process.

At each step I need to solve Ax=b with the same sparse big matrix A, but with different b.

What is more effective? 1. To solve Ax=b using pardiso at each step or 2. find invesrion A^{-1} and multiply at each step.

The size of A ~ 50 000x50 000.

thanks a lot,

Alex

5 Replies
Honored Contributor III
538 Views

If what you want to do is to solve linear equations, do not form the inverse! Forming the inverse is equivalent to solving A.X = B, where B = I, the identity matrix. In your case, that means solving with 50,000 right hand sides. In certain cases, it can also be less accurate to form and multiply by the inverse.

Call Pardiso with phase=12 to analyze and factor the matrix. Then, call Pardiso with phase = 33. If you know m values of b, you set nrhs = m in this call, after entering b1, b2, ... into columns-1, 2, ... of B. If you only know one value of b at a time, call with nrhs = 1 when you have that b ready and wish to solve. Remember that you can do this only if the matrix itself remains unchanged between calls.

New Contributor I
532 Views
Thanks a lot.

The problem is that phase=12  and phase = 33 will be realized in different functions. Where can I store the results of factorization (phase=12) and how to transfer it to further functions realizing phase=33.

Namely

func1
{
phase=12
PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
}

and

func2
{
phase=33
PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error);
}

As I understand, func1 I will called 1 time, when func2 will be called many times. Matrix a is constant, but right side b will be calculated in iteration process.

What arguments should I store and transfer from func1 into func2 in my case.

Thank you.
Honored Contributor III
525 Views

As long as you do not change arguments such as pt, n, a, ia, ja, between successive calls to Pardiso, the solver state will be preserved. Arguments with "dum" are not used in the calculation and are just place holders. Pardiso manages storing the factors of the matrix, the pivot array, etc., internally, and you do not need to participate in those tasks.

How do you preserve the arguments to Pardiso that do need to be saved? That depends on the programming language that you use and your programming style. In Fortran, for example, those could be variables that are placed in a module that is used in the various functions func1, func2, etc. In C, they could be global variables. Or, they could be added to the argument lists of func1, func2, etc. Or, you could put all the calls to Pardiso in one function. In between those calls, you could call your func1, func2, etc. to do their specific tasks.

New Contributor I
520 Views

Thaks a lot. The problem is solved.

Moderator
507 Views

Hi,