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

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.

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

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.

Link Copied

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

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.

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

Thank you very much.

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

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.

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

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

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