writing a banded matrix in a row major layout for lapack solver dgbsv

Question

I want to solve this linear system Ax=b, where:

A =  
    [-0.23  2.54  -3.66  0;  
     -6.98  2.46  -2.73  -2.73;  
      0     2.56   2.46   4.07;  
      0     0     -4.78   3.82]

b = [4.42 27.13 -6.14  10.5]

and the solution should be

x = [-2  3  1  -4]

which is a banded matrix with lower band equals to 1 and upper band equals to 2

It is solved using DGBSV solver as follows

#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"
#define N 4
int main() { 
    int i;
    MKL_INT ipiv[N];
    double a[16] = { 0,    -0.23,  2.54, -3.66,
                    -6.98,  2.46, -2.73, -2.13,
                     2.56,  2.46,  4.07,  0,
                    -4.78,  3.82,  0,     0};

    double b[4] =  {4.42, 27.13, -6.14, 10.5};

    LAPACKE_dgbsv( LAPACK_ROW_MAJOR, N, 1, 2, 1, a, N, ipiv, b, 1);

    for(i=0;i<N;i++)
        printf("%f\n", b[i]);
}

The code is aborting at the dgbsv solver. When I write matrices a and b pointers it gives the values of the address.

ewcz ewcz · Accepted Answer · 2017-03-13T11:10:22

For the input stated in your question, i.e.,

A =  
[-0.23  2.54  -3.66  0;  
 -6.98  2.46  -2.73  -2.73;  
  0     2.56   2.46   4.07;  
  0     0     -4.78   3.82]

b = [4.42 27.13 -6.14  10.5]

the solution (solving the system as a dense one) I get is:

[-3.77599, -1.28156, -1.85975, 0.421568]

However, as for the code, there are several things worth mentioning. What the function LAPACKE_dgbsv essentially does is that it checks the validity of the input and in turn calls the function LAPACKE_dgbsv_work. If this function detects that the supplied layout is LAPACK_ROW_MAJOR, it merely transposes its input (matrix, right-hand sides) and passes all this to LAPACKE_dgbsv which expects column-major version of the packed matrix.

So, if your code specifies LAPACK_ROW_MAJOR, the array a should contain row-major version of the packed matrix as specified in the link above. Also, perhaps more importantly, LAPACKE_dgbsv requires additional space in the array a so that it can store the result of the LU decomposition. Specifically, there have to be additional kl rows.

Thus

#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"
#define N 4
int main() {

    int i;
    MKL_INT ipiv[N];

    double b[4] =  {4.42, 27.13, -6.14, 10.5};

    double a[] = {
        0, 0, 0, 0,
        0, 0, -3.66, -2.73,
        0, 2.54, -2.73, 4.07,
        -0.23, 2.46, 2.46, 3.82,
        -6.98, 2.56, -4.78, 0};

    MKL_INT info = LAPACKE_dgbsv(LAPACK_ROW_MAJOR, N, 1, 2, 1, a, N, ipiv, b, 1);
    //printf("%d\n", info);

    for(i=0;i<N;i++) printf("%f\n", b[i]);
}

then produces:

which agrees with the solution obtained with the dense solver.

writing a banded matrix in a row major layout for lapack solver dgbsv

2 Answers