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

Showing results for

- Intel Community
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library & Intel® Math Kernel Library
- Cholesky factorization, traspose and inversion of sparse matrix.

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

Highlighted

Fiori

Beginner

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

10-27-2015
03:38 PM

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

Accepted Solutions

Highlighted

Ying_H_Intel

Employee

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

10-29-2015
12:36 AM

122 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

4 Replies

Highlighted

Ying_H_Intel

Employee

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

10-27-2015
07:35 PM

122 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

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

mecej4

Black Belt

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

10-28-2015
07:31 AM

122 Views

Highlighted

Fiori

Beginner

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

10-28-2015
10:53 AM

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

Highlighted

Ying_H_Intel

Employee

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

10-29-2015
12:36 AM

123 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

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