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.
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:
-3.775989
-1.281561
-1.859751
0.421568
which agrees with the solution obtained with the dense solver.