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

Cholesky factorization, traspose and inversion of sparse matrix.

Fiori
Beginner
1,843 Views

Hello!

Let Sigma be a sparse matrix. I would like to compute the Cholesky factorization of Sigma (the Upper(Lt) or lower triangular (L)), transpose it, and compute the folowing terms


w = inv(L)*mu;
m = inv(Lt)*w;
v = inv(Lt)*b;


where mu, b are known.

The problem I face is that I can't find the routines (and examples) when the matrix is sparse. Any help would be really appreciated.
Thank you very much.

 

 

 

0 Kudos
1 Solution
Ying_H_Intel
Employee
1,843 Views

Hi GF,

As I understand, both Lapack function with band format and pardiso can do your computation. You may try both.

But if your matrix format is band obviously, I evaluate a  band matrix and lapack function call may better (easy to use and better performance). 

 Pardiso support separate Forward and Backward Substitution, so you should be able to call phase 331 for different & multiply z . 

The solver execution step (see parameterphase= 33below) can be divided into two or three separate

substitutions: forward, backward, and possible diagonal. This separation can be explained by the examples of
solving systems with different matrix types.
A real symmetric positive definite matrix A(mtype= 2) is factored by Intel MKL PARDISO as A= L*L
T
. In
this case the solution of the system A*x=bcan be found as sequence of substitutions: L*y=b(forward
substitution, phase=331) andL
T
*x=y(backward substitution, phase=333)

If iparm[35] = 2, phases 331, 332, and 333 perform a different

decomposition:
You can supply a custom implementation for phase 332 instead of calling
pardiso. For example, it can be implemented with dense LAPACK
functionality. Custom implementation also allows you to substitute the
matrix Swith your own.
NOTE
For very large Schur complement matrices use LAPACK functionality to
compute the Schur complement vector instead of the Intel MKL
PARDISO phase 332 implementation.

Best Regards,

Ying

View solution in original post

0 Kudos
4 Replies
Ying_H_Intel
Employee
1,843 Views

Hello, 

What are your problem size and sparse matrix structure (like Non-zero number ratio, storage format) etc? 

I guess, the sparse solver and sparse BLAS in MKL may help.  Please refer to the manual of MKL  or

 https://software.intel.com/en-us/articles/intel-mkl-sparse-solvers-training-material  ( i will attach the pdf file there). 

and you can find code sample in MKL install forder : solverc or solverf (fortran). 

Sparse Solver Routines=> Intel MKL PARDISO

Best Regards,

Ying 

 

0 Kudos
mecej4
Honored Contributor III
1,843 Views

If all that is required is to solve A.x = b given a sparse matrix A and vectors x and b, it is neither necessary nor efficient to form and use the explicit inverse of A. When you use a sparse linear solver, the usual sequence of operations is (i) analyze (ii) factorize and (iii) solve. If the matrix is banded, you can use the band solver routines in MKL/Lapack. Otherwise, use the sparse solver Pardiso directly or through the DSS interface. For details, consult the MKL documentation and the example source and data files in the MKL examples/solverf directory.

0 Kudos
Fiori
Beginner
1,843 Views

Hello. Thank you very much for your help.

I want to sample from a Multivariate Normal distribution(Q^{-1}b, Q^{-1}) , where Q is a triagonal, symmetric matrix:

the elements of the main diagonal are {a1, a2, a2, ...., a2, a1} and the elements of the upper (down) diagonal are {a3, a3, a3, .... a3}.

The algorithm in order to do this is the following:

1. Compute the Cholesky factorization, Q=LU, where U=L^{T}

2. Solve Lw = b

3. Solve Um=w

4. Sample z~N(0,1)

5. Solve Uv=z

6.Compute x=m+y

7.Return x

The number of rows of Q is 5e+5. 

Steps 2 and 3 gives the solution of Qm=b. The function PARDISO computes it. How can I compute step 5? Is it better to write it as a band or as a sparse matrix?   

 

0 Kudos
Ying_H_Intel
Employee
1,844 Views

Hi GF,

As I understand, both Lapack function with band format and pardiso can do your computation. You may try both.

But if your matrix format is band obviously, I evaluate a  band matrix and lapack function call may better (easy to use and better performance). 

 Pardiso support separate Forward and Backward Substitution, so you should be able to call phase 331 for different & multiply z . 

The solver execution step (see parameterphase= 33below) can be divided into two or three separate

substitutions: forward, backward, and possible diagonal. This separation can be explained by the examples of
solving systems with different matrix types.
A real symmetric positive definite matrix A(mtype= 2) is factored by Intel MKL PARDISO as A= L*L
T
. In
this case the solution of the system A*x=bcan be found as sequence of substitutions: L*y=b(forward
substitution, phase=331) andL
T
*x=y(backward substitution, phase=333)

If iparm[35] = 2, phases 331, 332, and 333 perform a different

decomposition:
You can supply a custom implementation for phase 332 instead of calling
pardiso. For example, it can be implemented with dense LAPACK
functionality. Custom implementation also allows you to substitute the
matrix Swith your own.
NOTE
For very large Schur complement matrices use LAPACK functionality to
compute the Schur complement vector instead of the Intel MKL
PARDISO phase 332 implementation.

Best Regards,

Ying

0 Kudos
Reply