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

[Fiori]

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.

[mecej4 · Accepted solution]

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.

[Fiori]

 Thank you very much.

[Fiori]

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.

[mecej4]

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