MATLAB: Problem using BLAS routine ztbsv.c

blascmexcomplex numbersdoublecomplexztbsv

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 2×2 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!

Best Answer

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.
Related Question