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

solve system when matrix is banded, symmetric, positive definite matrix.

Fiori
Beginner
1,147 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.

 

0 Kudos
1 Solution
mecej4
Honored Contributor III
1,147 Views

How about this:

  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.

View solution in original post

0 Kudos
4 Replies
mecej4
Honored Contributor III
1,148 Views

How about this:

  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.

0 Kudos
Fiori
Beginner
1,147 Views

 Thank you very much.

0 Kudos
Fiori
Beginner
1,147 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.

0 Kudos
mecej4
Honored Contributor III
1,147 Views

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

0 Kudos
Reply