MKL PARDISO, get L U factors

[DmitrySmi]

Hello!

Is there a way, maybe 'semi-official', to extract L and U factors from PARDISO after the numerical factorization is done?

Here is what I am trying to do.
I need to compute the following product
A'*inv(Y)*A
Where Y is a huge sparse positive definite matrix. n-by-n, where n is several million, sometimes ten million, sometimes more.
A is also sparse, n-by-k, where k is several thousand.
(Just in case this is from the circuit analysis, Y is the nodal admittance matrix, A is the adjacency matrix and this product is the impedance matrix of the selected subset of the circuit edges).
The most efficient way to compute this product (one I can think of) is to decompose Y
Y = L*L'
So we have
A'*inv(L*L')*A
or
A'*inv(L')*inv(L)*A
Then compute
C=inv(L)*A
And, finally,
C'*C

Once Y has been decomposed, this
C=inv(L)*A
Can be computed efficiently by sparse triangular solve.
The efficient sparse triangular solve algorithm is known (described for example in Direct Sparse Methods by Davis), and I know that C is going to be sparse (I checked it for a number of examples)
The problem however it that PARDISO neither implements sparse triangular solve, nor gives access to L factor.

Any suggestions/ideas?

Thanks!

[mecej4]

From the documentation of the Basel/Lugano Pardiso 7.2, it seems to me, that version may have the capability that you seek -- see page 17 of the manual, under "IPARM (26) — Splitting of Forward/Backward Solve." 

I have not actually tried this suggestion myself, but you can judge for yourself whether it offers what you seek.

[DmitrySmi]

Thanks.

Yes it allows to split the forward and the backward solve (and scaling too), but still it can only do dense forward/backward solve.

I need to do C=L\A, where a is large (say 1,000,000-by-10,000) but is very sparse -- it is an incidence matrix and has 2 or 1 entries per row, thus up to 2,000,000 nonzeros in total. C is known to be sparse as well. Dense A and C do not even fit in memory, while the sparse forward solve can be done efficiently, but PARDISO neither supports it nor gives access to the L factor.

So far my hope is that L is available in the solver data (pt array of pointers) in some common format, so I can extract it...

[Kirill_V_Intel]

Hi Dmitry,

First of all,  (especially for @mecej4) oneMKL PARDISO also has the feature of splitting the solving phase into multiple phases, there is no need to switch to the other solver to use this feature. In oneMKL PARDISO one can do only one of the solving phases by specifying phase = 331,332 or 333 (for two triangular solves and the diagonal one).

Second, there is no functionality currently to extract factors from oneMKL PARDISO (but it certainly could be done in case we get a feature request).

However, this feature ( = exporting factors in CSR format) is already enabled for the distributed cluster sparse solver via cluster_sparse_solver_export functionality. As a workaround, I suggest you try it. You can use the cluster sparse solver with a single MPI process and whichever #OpenMP threads you want and you will likely get similar (maybe even better in some cases) performance compared to non-distributed PARDISO solver.

There are certain limitations on when this feature can be used. For more information, see https://software.intel.com/content/www/us/en/develop/articles/export-functionality-in-the-cluster-sparse-solver.html, example cl_solver_export_c.c and documentation https://software.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/parallel-direct-sparse-solver-for-clusters-interface/cluster-sparse-solver-export.html

Lastly, I'd like to ask a general question: do you really need the matrix C as a sparse matrix? Often it is enough to be able to compute matrix-vector product and then you can comput v = C * u as w = A*u, L v = w (and then just use the backward solve phase of PARDISO).

Best,
Kirill

[DmitrySmi]

Kirill, awesome, thanks!

This is definitely something to play with!

A quick question. I can see the following two export operations are defined: SPARSE_PTLUQT, SPARSE_DPTLUQT.

Does it mean that if I do a symmetric factorization (hermitian, or complex symmetric) the export functionality is not available?

And yes, it appears that I really do need the matrix C as the sparse matrix. More details are above, but briefly at the end I need to compute the product Z = C'*C.

Thanks again, Dmitry

[Kirill_V_Intel]

Dmitry,

Sorry, I've missed the details about your case when I asked my question about the use case. Now I understand it better.

The symmetric case should work (think of it as U = L^T is a speciic case of U). 

Let us know in case the export feature helps/doesn't help/doesn't work, your feedback will be very welcome.

Best,
Kirill

[DmitrySmi]

Kirill, thanks, I will get back once I try this.

[DmitrySmi]

Hi @Kirill_V_Intel and Everybody,

Thanks again for the advice, I finally got a chance to try this.

I seem to be able to get it work -- exported L and U factors in a non-symmetric case are correct, but when factorizing symmetric matrices the results are a bit unexpected.

To demonstrate, I modified the cl_solver_export_c.c from the MKL distribution. I replaced the matrix with a really simple 2-by-2 one, and added printing of the exported factors.

The modified example is enclosed to this post, and here is the output I am getting (U data is all zeros so it is not being printed):

Matrix A (upper triangle):
A(0,0) = 12.000000
A(0,1) = -4.000000
A(1,1) = 5.000000
Matrix L:
L(0,0) = 12.000000
L(1,0) = -0.333333
L(1,1) = 3.666667
Permutation vector P:
P(0) = 0
P(1) = 1
Permutation vector Q:
Q(0) = 0
Q(1) = 1

It can be easily seen that L*L' does not equal A (both permutations are identities so we ignore them).

This is the case for all the examples I tried so far, so I wonder what am I doing wrong...

Thanks,

Dmitry

[Gennady_F_Intel]

Dmitry,

The fix of the issue available in the official release of MKL 2021.3 which has happened yesterday. You could take to try and let us know the results.

[Gennady_F_Intel]

validation of Intel MKL v.2021.3 results:

Allocating data for the factors L and U and permutations P and Q...

Saving allocated pointers inside the cluster sparse solver...

Calling the main export functionality...

Main export functionality has been caled successfully...

Matrix A (upper triangle): 

 A(0,0) = 12.000000 

 A(0,1) = -4.000000 

 A(1,1) = 5.000000 

Matrix L: 

 L(0,0) = 3.464102 

 L(1,0) = -1.154701 

 L(1,1) = 1.914854 

Permutation vector P: 

 P(0) = 0 

 P(1) = 1 

Permutation vector Q: 

 Q(0) = 0 

 Q(1) = 1 

 TEST PASSED

As we can see L*L^t == A,

the issue is fixed.

[Gennady_F_Intel]

The issue is closing 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.

[DmitrySmi]

Thanks a lot for the fix!