topic How about this: in Intel® oneAPI Math Kernel Library

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

Fiori — Fri, 20 May 2016 15:25:43 GMT

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:

Compute the Cholesky factorization, Q=LU, where U=L^{T}
Solve Lw = b
Solve Um=w
Sample z~N(0,1)
Solve Uv=z
Compute x=m+v
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.

How about this:

mecej4 — Fri, 20 May 2016 15:57:00 GMT

How about this:

Factorize by calling ?PBTRF.
Call ?TBTRS with UPLO='L' to solve Lw = b.
Sample z~N(0,1).
Compute y = w+z.
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.

Thank you very much for your

Fiori — Fri, 20 May 2016 20:08:00 GMT

Thank you very much.

I would like to ask you, if

Fiori — Mon, 23 May 2016 15:27:58 GMT

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.

You are correct on both

mecej4 — Mon, 23 May 2016 16:25:38 GMT

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