- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi
I am quite newto MKL andLapack, so I was trying to call spbstf in c++
but I am not sure how I should store the lower lower triangular part of the original matrix
Say I have the following matrix
1 0.3 0.3 0.3
0.3 10.3 0.3
0.30.3 10.3
0.30.3 0.31
If I want to call the following to get the cholesky split
dpbstf(uplo, n, kb, bb, ldbb, info)
Is bb a one dimensional or two dimensional array? how should lower triangular part stored in array bb?
Anyone can give me a runnable example?
Link Copied
1 Reply
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
http://www.netlib.no/netlib/lapack/double/dpbstf.f
(or same page from USA netlib)
shows this as band storage scheme, a 2-d Fortran array with each column starting at the main diagonal. /1.,.3,.3,.3,1,.3,.3,,1.,.3,.... if ldbb==4
Of course, that's a 1-d array in the C/C++ view, with the ldbb (by reference) used to mark the partition into columns.
Cholesky square root method is not the most efficient for this application. You may be able to beat MKL with your own optimized code, and you are entitled to use the public code as a starting point, provided you give credit.
The source code would be a little easier to follow if converted to f90 with named loops.
Unless someone can find a specific description of the symmetric band storage format used in netlib, a look at the source code is indispensable.
Simply saying band storage format doesn't do the job; IMSL Fortran band storage format was backwards and not properly optimizable, because of mis-translation from original published Algol code.
(or same page from USA netlib)
shows this as band storage scheme, a 2-d Fortran array with each column starting at the main diagonal. /1.,.3,.3,.3,1,.3,.3,,1.,.3,.... if ldbb==4
Of course, that's a 1-d array in the C/C++ view, with the ldbb (by reference) used to mark the partition into columns.
Cholesky square root method is not the most efficient for this application. You may be able to beat MKL with your own optimized code, and you are entitled to use the public code as a starting point, provided you give credit.
The source code would be a little easier to follow if converted to f90 with named loops.
Unless someone can find a specific description of the symmetric band storage format used in netlib, a look at the source code is indispensable.
Simply saying band storage format doesn't do the job; IMSL Fortran band storage format was backwards and not properly optimizable, because of mis-translation from original published Algol code.
Reply
Topic Options
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page