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

## solve system when matrix is banded, symmetric, positive definite matrix. Beginner
330 Views

Hello.

Let Q be a banded, symmetric, positive definite matrix. (The number of rows of Q is 2e+5) . I want to do 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+v
7. Return x

Steps 2 and 3 give the solution of Qm=b.

I have asked this question here, but then Q was a triagonal, symmetric matrix. I used fuctions  'LAPACKE_dpbtrf, cblas_dtbsv' and solved the problem,

```   /* Colesky factorization */
info = LAPACKE_dpbtrf(LAPACK_COL_MAJOR, 'U', dim+1, 1, Sigmab, 2 );
if(info!= 0){mexPrintf( "C++ error: Cholesky failed");  }

/* step 2*/
cblas_dtbsv(CblasColMajor, CblasUpper, CblasTrans, CblasNonUnit, dim, 1, Sigmab, 2, y1, 1);

/* step 3*/
cblas_dtbsv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, dim, 1, Sigmab, 2, y1, 1);

/* step 5 */
cblas_dtbsv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, dim, 1, Sigmab 2, y2, 1);  ```

I found functions 'LAPACKE_dpbtrf, LAPACKE_dpbtrs' that give the solution of steps 2 and 3.

```   info = LAPACKE_dpbtrf(LAPACK_COL_MAJOR, 'L', dim, p, Sigma, p+1);
if(info!= 0){mexPrintf( "C++ error: Cholesky failed");  }

info = LAPACKE_dpbtrs(LAPACK_COL_MAJOR, 'L', dim, p, NRHS, Sigma, p+1, y1, dim);
if(info!= 0){mexPrintf( "C++ error: the execution is not successful");  }```

Firstly, I would like to ask if there is a better way to solve this and secondly  I don't know how to find the solution of step 5.

Thank you very much.

1 Solution Black Belt
330 Views

1. Factorize by calling ?PBTRF.
2. Call ?TBTRS with UPLO='L' to solve Lw = b.
3. Sample z~N(0,1).
4. Compute y = w+z.
5. Call ?TBTRS with UPLO='U' to solve Ux = y and obtain the solution x.

​Note that I have collapsed the two calls in your Steps 3 and 5 into a single call in my Step 5. The computational savings are modest but it is the spirit that counts.

My Step 5 also contains the answer to  your statement "I don't know how to find the solution of step 5".

P.S.: I notice that you are coding in C and using the LapackE interface. Please use the LapackE versions of the Lapack routines that I named.

4 Replies Black Belt
331 Views

1. Factorize by calling ?PBTRF.
2. Call ?TBTRS with UPLO='L' to solve Lw = b.
3. Sample z~N(0,1).
4. Compute y = w+z.
5. Call ?TBTRS with UPLO='U' to solve Ux = y and obtain the solution x.

​Note that I have collapsed the two calls in your Steps 3 and 5 into a single call in my Step 5. The computational savings are modest but it is the spirit that counts.

My Step 5 also contains the answer to  your statement "I don't know how to find the solution of step 5".

P.S.: I notice that you are coding in C and using the LapackE interface. Please use the LapackE versions of the Lapack routines that I named. Beginner
330 Views

Thank you very much. Beginner
330 Views

I would like to ask you, if you mean at step 2 and 5 to use the function  'LAPACKE_dtbtrs' and not '?PBTRS'.?

Thank you very much. Black Belt
330 Views

You are correct on both counts. I have corrected #2 accordingly. 