Solving a banded system from C using LAPACK's DGBSV

Question

I have tried to search for a similar thread related to this topic, but nobody seems to care for banded systems (...).

I am interested in solving a banded matrix using LAPACK/ScaLAPACK from a C code. First, I want to achieve a sequential solution with LAPACK, before attempting anything with ScaLAPACK.

Problem: The row-major/column-major difference between both languages seems to be affecting my solution process. Here's the system I intend to solve:

The following code, translates that matrix into LAPACK's banded data structure, specified in here.

  int rr = 6;     // Rank.
  int kl = 2;     // Number of lower diagonals.
  int ku = 1;     // Number of upper diagonals.
  int nrhs = 1;   // Number of RHS.
  double vals[36] = {0.0,   0.0,   0.0,   0.0,   0.0,   0.0,  // Req. ex. space.
                     0.0,   0.0,   0.0,   0.0,   0.0,   0.0,  // Req. ex. space.
                   666.0,   0.0,   0.0,   0.0,   0.0,  22.5,  // First up diag.
                     1.0, -50.0, -50.0, -50.0, -50.0,  -2.6,  // Main diagonal.
                    27.5,  27.5,  27.5,  27.5,   4.0, 666.0,  // First low diag.
                     0.0,   0.0,   0.0,  -1.0, 666.0, 666.0}; // 2nd low diag.

  int lda = rr;   // Leading dimension of the matrix.
  int ipiv[6];    // Information on pivoting array.
  double rhs[] = {1.0, 1.0, 1.0, 1.0, 1.0, 0.0};  // RHS.
  int ldb = lda;  // Leading dimension of the RHS.
  int info = 0;   // Evaluation variable for solution process.
  int ii;         // Iterator.
  int jj;         // Iterator.

  dgbsv_(&rr, &kl, &ku, &nrhs, vals, &lda, ipiv, rhs, &ldb, &info);

  printf("info = %d\n", info);
  for (ii = 0; ii < ldb; ii++) {
    printf("%f\n", rhs[ii]);
  }
  putchar('\n');

As I said, I am worried that the way I am translating my matrix is incorrect given the col-major nature, as well as the indexing nature of Fortran, since my solution yields:

[ejspeiro@node01 lapack-ex02]$ make runs
`pwd`/blogs < blogs.in
info = 1
1.000000
1.000000
1.000000
1.000000
1.000000
0.000000

The return value from Fortran info = 1 implies that factorization was completed, but U(1,1) = 0 in the LU factorization of A = LU.

Answer 1

All right, I will answer this in order to sort of call it "solved".

These files, represent the functioning code. I am making this available in case anybody has the same problem: solving a banded system of equations, using LAPACK from C.

1. https://dl.dropbox.com/u/5432016/banded/cmtk.h
2. https://dl.dropbox.com/u/5432016/banded/blogs.c

And if anybody has any extra suggestions as far as the implementation goes, I will gladly welcome them!

My next step is to distribute the core matrices, in order to solve with ScaLAPACK.

Thanks!

Answer 2

As you noticed, your matrix was entered in row-major format. Writing it in column major, the rows we have visually will correspond to the columns:

  double vals[36] = {0.0,   0.0,   0.0,   1.0,  27.5,   0.0,
                     0.0,   0.0,   0.0, -50.0,  27.5,   0.0,
                     0.0,   0.0,  22.5, -50.0,  27.5,   0.0,
                     0.0,   0.0,  22.5, -50.0,  27.5,  -1.0,
                     0.0,   0.0,  22.5, -50.0,   4.0,  0.0,
                     0.0,   0.0,  22.5,  -2.6,   0.0,  0.0};