Problem using BLAS routine ztbsv.c

2 visualizzazioni (ultimi 30 giorni)
Hi there, I have a problem using the BLAS routine ztbsv.c in a c-mex file to solve a linear equation system Ax=b, with A being a quadratic, triangular, band matrix stored in a "packed band storage" as follows:
In the real case, the triangular band matrix
1 5 9 0 0
0 2 6 10 0
0 0 3 7 11
0 0 0 4 8
0 0 0 0 5
Is stored as
0 0 9 10 11
0 5 6 7 8
1 2 3 4 5
with the top left 2x2 triangle not being referenced.
The following C-code works perfectly:
void mexFunction( int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[] )
{
double br[5]={1,2,3,4,5};
double Ar[15]={0,0,1,0,5,2,9,6,3,10,7,4,11,8,5};
char UH = 'U';
char EN = 'N';
ptrdiff_t n = 5;
ptrdiff_t k = 2;
ptrdiff_t l = k+1;
ptrdiff_t one = 1;
dtbsv_(&UH,&EN,&EN,&n,&k,Ar,&l,br,&one);
return;
}
However, the complex case is not working. According to the Matlab documentation, a complex vector of length n is stored as a real vector of length 2n, with real and imaginary part alternating, for example:
1+2i
3+4i
5+6i
is stored as
{1,2,3,4,5,6}
All complex BLAS routines are working fine this way, except for ztbsv. I tried the same storage scheme as in the real case, with real and imaginary parts alternating. If I solve Ax=b this way, some components of x (the last ones) are correct, but not all.
I haven't found a explanation how to properly store the matrix. The official BLAS documentation describes a different implementation of the complex routines. Here, a structure "doublecomplex", containing the two doubles real and imaginary part, is defined. But the blas.h file in my Matlab directory defines complex arrays as described above: as a double array of two times the length, alternating real and imaginary part. And again: All routines I tried work, except for ztbsv.
Any help is highly appreciated!
  2 Commenti
Steven Lord
Steven Lord il 6 Lug 2015
Modificato: Steven Lord il 6 Lug 2015
"According to the Matlab documentation, a complex vector of length n is stored as a real vector of length 2n, with real and imaginary part alternating"
Where in the documentation do you see this? If that is in the documentation for MATLAB, it needs to be corrected as that is incorrect. The real part and the imaginary part are stored separately. In a MEX file, you use the mxGetPr function to access the pointer to the real data and mxGetPi to access the pointer to the imaginary data.
James Tursa
James Tursa il 6 Lug 2015
Modificato: James Tursa il 6 Lug 2015
@Steven: Christoph is correct. He is talking specifically about the BLAS/LAPACK interface for complex variables, which always use the Fortran standard of interleaving the real and imaginary parts in memory. If one is passing a complex MATLAB variable to a BLAS/LAPACK routine via mex, you have to copy the data into this format first, and then copy the results back to MATLAB format after.

Accedi per commentare.

Risposta accettata

James Tursa
James Tursa il 5 Lug 2015
Modificato: James Tursa il 5 Lug 2015
All of the BLAS and LAPACK complex storage is based on an original Fortran standard, which has real and imaginary parts interleaved. The C structure doublecomplex is intended to mimic this storage scheme, and the blas.h header file should as well. Early versions of MATLAB did not have correct headers for the complex arguments of BLAS/LAPACK functions. I think that has been corrected but have not actually checked lately.
Your example doesn't make sense, and I assume it is a typo. This:
1+2i
2+3i
4+6i
Should be stored in memory as:
{1,2,2,3,4,6}
Which is not what you have written.
I think you will need to post your complex code so we can look at it. Do you have both Ar and br as doublecomplex storage memory? I would also look at the blas.h and lapack.h header files for your version of MATLAB to see if they are as expected (possibly could be incorrect as noted above). There is always the chance that the library function itself has bugs in it, but I wouldn't assume that without checking everything else first. You could also download a C source code for ztbsv and call that routine (i.e., bypass the library function) to see if things work properly.
  1 Commento
Christoph Nicolai
Christoph Nicolai il 5 Lug 2015
I use Matlab R2014b and the blas.h file is not using doublecomplex but a double array of length 2n for a complex vector of length n.
But I have found a little mistake I made in the packed store scheme when I was preparing the complex code to post it here. I double checked it before, but didn't catch this mistake.
I guess it works now, thank you for your time and help!

Accedi per commentare.

Più risposte (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by