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

Question on calling spbstf in VC++

zhouxingchi
Beginner
742 Views

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?

0 Kudos
1 Reply
TimP
Honored Contributor III
742 Views
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.
0 Kudos
Reply