Question on calling spbstf in VC++

zhouxingchi — Thu, 20 Mar 2008 06:53:59 GMT

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?

Re: Question on calling spbstf in VC++

TimP — Thu, 20 Mar 2008 17:08:09 GMT

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.

topic Re: Question on calling spbstf in VC++ in Intel® oneAPI Math Kernel Library

Question on calling spbstf in VC++

Re: Question on calling spbstf in VC++