Problems using ZGBSV

[Diaboflo]

Hi,
I'm trying to solve a general banded complex matrix with LAPACK. I use C=++ and cross compile it with f77. I managed to save the matrix in band storage & column major order but I still get wrong results.
Can anyone help me out what I'm doing wrong?
I was trying it with just a simple matrix to get it running.
I'm working with linux and g++ compiler

Compilation:

Code: Select all

slioch:temp> g++ -c zgbsv.cpp
slioch:temp> f77 -o zgbsv zgbsv.cpp /home/andrew/lib/lapack_LINUX.a -lstdc++ -lm /home/andrew/lib/blas_LINUX.a -lm


I'm filling the band elements just with the values A(1,3)=13, A(3,2)=32 etc. For my easier understanding I first create the matrix then transpose it and then convert it to band storage. My rhs Vector is just zero so I would expect the trivial solution which I don't get. INFO takes the value three which means U(3,3) is zero (I don't really know what to do with this information).
Can anyone help me out with that?
Thanks, Florian

My code:

Code: Select all


#include <iostream>
#include <iomanip>
#include <string>
#include <fstream>
#include <math.h>

using namespace std;

struct cplxnag{double re; double im;};
extern "C" void zgbsv_(int *, int *, int *, int *, cplxnag *,int *, int *, cplxnag *, int *, int *);


int main(){
   cplxnag *A,*Ad, *AB, *B;
   int *IPIV;
   cplxnag eigVecL[9][9],eigVecR[9][9], eigVal[9];
   double lWork[6],rWork[6];

   int N = 8;
   int KL  = 2;
   int KU = 1;
   int NRHS = 1;   
   int LDAB = 2*KL+KU+1;
   AB = (cplxnag *)malloc(sizeof(cplxnag)*N*LDAB);
   IPIV = (int *)malloc(sizeof(int)*N);
   int LDB = N;
   B = (cplxnag *)malloc(sizeof(cplxnag)*LDB*NRHS);
   int INFO = 0;

   for (int i=0;i<LDB;i++) {B[i].re =0.; B[i].im =0.;}
   for (int i=0;i<LDB;i++) {cout << B[i].re<<"   ";}
   cout <<endl;   
   
   A = (cplxnag *)malloc(sizeof(cplxnag)*N*N);
   Ad = (cplxnag *)malloc(sizeof(cplxnag)*N*N);


// fill the matrix
   for (int i=0;i<N;i++){
      for (int j=0; j<N;j++){
         if (i==j){
            A[j*N+i].re=(j+1)*10+i+1;
            A[j*N+i].im=0;
         }
         else if((i-j)==1.){
               A[j*N+i].re=(1+j)*10+i+1;
               A[j*N+i].im=0;   
            }
         else if(i-j==-1){
               A[j*N+i].re=(1+j)*10+i+1;
               A[j*N+i].im=0;   
            }
         else if(i-j==-2){
               A[j*N+i].re=(1+j)*10+i+1;
               A[j*N+i].im=0;   
            }
         else{
            A[j*N+i].re=0;
            A[j*N+i].im=0;            
         }
      }
   }

//matrix output
   cout <<"Matrix A\n";
   for (int i=0;i<N;i++){
      for (int j=0; j<N;j++){
         cout  << setw(2)<<setfill('0')<<A[i*N+j].re<< "   ";
         Ad[i*N+j].re = 0.;
      }
   cout << endl;
   }

//transpose matrix
   cout <<"Matrix A transposed\n";
   for (int i=0;i<N;i++){
      for (int j=0; j<N;j++){
       Ad[i*N+j]=A[j*N+i];
       if (j<LDAB) AB[i*N+j].re=0.;
       cout  << setw(2)<<setfill('0')<< Ad[i*N+j].re<< "   ";
      }
   cout << endl;
   }
   cout <<endl<<endl;

// band storage
   for (int j=0;j<=N;j++){
      for (int i=max(j-KU,0);i<=min(j+KL,N-1);i++){
         AB[(KL+KU+i-j)+j*N].re=Ad[j*N+i].re;
         AB[(KL+KU+i-j)+j*N].im=Ad[j*N+i].im;
      }
   }
   cout<<"Matrix AB\n";
   for (int i=0;i<LDAB;i++){
      for (int j=0;j<N;j++){
         cout << setw(2)<<setfill('0')<< AB[i+j*N].re<<"   ";
      }
   cout << endl;
   }
         
//solve matrix
   zgbsv_(&N, &KL, &KU, &NRHS, AB, &LDAB, IPIV, B, &LDB, &INFO);
   cout<<"Matrix AB after solving \n";
   for (int i=0;i<LDAB;i++){
      for (int j=0;j<N;j++){
         cout << setw(12)<<setfill(' ')<< AB[i+j*N].re<<"   ";
      }
   cout << endl;
   }
   cout <<"Solution:\n";
   for (int i=0;i<N;i++) cout << B[i].re<<"  "<<B[i].im<<endl;
   cout <<"Info = "<<INFO<<endl;
   
}




My output:

Code: Select all

B
0   0   0   0   0   0   0   0
Matrix A
11   12   00   00   00   00   00   00
21   22   23   00   00   00   00   00
31   32   33   34   00   00   00   00
00   42   43   44   45   00   00   00
00   00   53   54   55   56   00   00
00   00   00   64   65   66   67   00
00   00   00   00   75   76   77   78
00   00   00   00   00   86   87   88
Matrix A transposed
11   21   31   00   00   00   00   00
12   22   32   42   00   00   00   00
00   23   33   43   53   00   00   00
00   00   34   44   54   64   00   00
00   00   00   45   55   65   75   00
00   00   00   00   56   66   76   86
00   00   00   00   00   67   77   87
00   00   00   00   00   00   78   88


Matrix AB
00   00   00   00   00   00   00   00
00   00   00   00   00   00   00   00
00   12   23   34   45   56   67   78
11   22   33   44   55   66   77   88
21   32   43   54   65   76   87   00
31   42   53   64   75   86   00   00
Matrix AB after solving
           0             12              0             34             45              0              0              0
           0             22              0        6.58065             55              0   8.48798e-314              0
           0      -0.193548             43         12.563        -0.5625             76    1.6976e-313             78
          31      -0.369501             53             64      -0.690341             86             77             88
    0.677419             32             -5        0.84375              0         -3.875             87              0
    0.354839              0       -9.54545         0.6875              0       -4.75568              0              0
Solution:
77  0
87  0
0  0
0  0
0  0
0  0
0  0
78  0
Info = 3


[Julien Langou]

Almost ...
There are two mistakes in your code.

1. This is j*LDAB not j*N when your moving in column in
   AB. So want to write:
   

   Code: Select all

   // band storage
   for (int j=0;j<=N;j++){
      for (int i=max(j-KU,0);i<=min(j+KL,N-1);i++){
         AB[(KL+KU+i-j)+j*LDAB].re=Ad[j*N+i].re;
         AB[(KL+KU+i-j)+j*LDAB].im=Ad[j*N+i].im;
      }
   }


   
   And the same when you print:
   

   Code: Select all

      cout<<"Matrix AB\n";
   for (int i=0;i<LDAB;i++){
      for (int j=0;j<N;j++){
         cout << setw(2)<<setfill('0')<< AB[i+j*LDAB].re<<"   ";
      }
   cout << endl;
   }


   
   And again at:
   

   Code: Select all

      cout<<"Matrix AB after solving \n";
   for (int i=0;i<LDAB;i++){
      for (int j=0;j<N;j++){
         cout << setw(12)<<setfill(' ')<< AB[i+j*LDAB].re<<"   ";
      }
   cout << endl;
   }


2. Second when you are initializing AB in the loop, the statement (j<LDAB) should
   be (i<LDAB) or something like this. This gives
   

   Code: Select all

      cout <<"Matrix A transposed\n";
   for (int i=0;i<N;i++){
      for (int j=0; j<N;j++){
       Ad[i*N+j]=A[j*N+i];
       if (i<LDAB) AB[i*N+j].re=0.;
       cout  << setw(2)<<setfill('0')<< Ad[i*N+j].re<< "   ";
      }
   cout << endl;
   }
   cout <<endl<<endl;


   
   Otherwise you're messing up your memory ...
   Actually you do not need to initialize AB with 0.0e+00 unless to have nice print.
   LAPACK will overwrite those values, or do not touch them.
   
3. Thrid, well, testing with the zero RHS is fine but you won't see much except that the
   matrix is singular or not .... So put a nontrivial RHS.
   

I have edited your code and tested, this works fine. Try:

Code: Select all



#include <iostream>
#include <iomanip>
#include <string>
#include <fstream>
#include <math.h>

using namespace std;

struct cplxnag{double re; double im;};
extern "C" void zgbsv_(int *, int *, int *, int *, cplxnag *,int *, int *, cplxnag *, int *, int *);


int main(){
   cplxnag *A,*Ad, *AB, *B;
   int *IPIV;
   cplxnag eigVecL[9][9],eigVecR[9][9], eigVal[9];
   double lWork[6],rWork[6];

   int N = 8;
   int KL  = 2;
   int KU = 1;
   int NRHS = 1;   
   int LDAB = 2*KL+KU+1;
   AB = (cplxnag *)malloc(sizeof(cplxnag)*N*LDAB);
   IPIV = (int *)malloc(sizeof(int)*N);
   int LDB = N;
   B = (cplxnag *)malloc(sizeof(cplxnag)*LDB*NRHS);
   int INFO = 0;

   for (int i=0;i<LDB;i++) {B[i].re =1.0e+00; B[i].im =0.;}
   for (int i=0;i<LDB;i++) {cout << B[i].re<<"   ";}
   cout <<endl;   
   
   A = (cplxnag *)malloc(sizeof(cplxnag)*N*N);
   Ad = (cplxnag *)malloc(sizeof(cplxnag)*N*N);

// fill the matrix
   for (int i=0;i<N;i++){
      for (int j=0; j<N;j++){
         if (i==j){
            A[j*N+i].re=(j+1)*10+i+1;
            A[j*N+i].im=0;
         }
         else if((i-j)==1.){
               A[j*N+i].re=(1+j)*10+i+1;
               A[j*N+i].im=0;   
            }
         else if(i-j==-1){
               A[j*N+i].re=(1+j)*10+i+1;
               A[j*N+i].im=0;   
            }
         else if(i-j==-2){
               A[j*N+i].re=(1+j)*10+i+1;
               A[j*N+i].im=0;   
            }
         else{
            A[j*N+i].re=0;
            A[j*N+i].im=0;           
         }
      }
   }

//matrix output
   cout <<"Matrix A\n";
   for (int i=0;i<N;i++){
      for (int j=0; j<N;j++){
         cout  << setw(2)<<setfill('0')<<A[i*N+j].re<< "   ";
         Ad[i*N+j].re = 0.;
      }
   cout << endl;
   }

//transpose matrix
   cout <<"Matrix A transposed\n";
   for (int i=0;i<N;i++){
      for (int j=0; j<N;j++){
       Ad[i*N+j]=A[j*N+i];
       if (i<LDAB) AB[i*N+j].re=0.;
       cout  << setw(2)<<setfill('0')<< Ad[i*N+j].re<< "   ";
      }
   cout << endl;
   }
   cout <<endl<<endl;

// band storage
   for (int j=0;j<=N;j++){
      for (int i=max(j-KU,0);i<=min(j+KL,N-1);i++){
         AB[(KL+KU+i-j)+j*LDAB].re=Ad[j*N+i].re;
         AB[(KL+KU+i-j)+j*LDAB].im=Ad[j*N+i].im;
      }
   }
   cout<<"Matrix AB\n";
   for (int i=0;i<LDAB;i++){
      for (int j=0;j<N;j++){
         cout << setw(2)<<setfill('0')<< AB[i+j*LDAB].re<<"   ";
      }
   cout << endl;
   }
         
//solve matrix
   zgbsv_(&N, &KL, &KU, &NRHS, AB, &LDAB, IPIV, B, &LDB, &INFO);
   cout<<"Matrix AB after solving \n";
   for (int i=0;i<LDAB;i++){
      for (int j=0;j<N;j++){
         cout << setw(12)<<setfill(' ')<< AB[i+j*LDAB].re<<"   ";
      }
   cout << endl;
   }
   cout <<"Solution:\n";
   for (int i=0;i<N;i++) cout << B[i].re<<"  "<<B[i].im<<endl;
   cout <<"Info = "<<INFO<<endl;
   
}


Solution x to your Ax=b system when b has just ones should be:

Code: Select all

        0.010006405414537
        0.074160795036674
       -0.036594434978788
       -0.013991989844831
        0.001654552571558
        0.064358430608231
       -0.036712163988438
       -0.015237122333111


-j

[Diaboflo]

Great! Thanks for the quick help, it works fine. Now I can apply this on my 8000*8000 matrix :-)
Florian