dgbsv() and "double* ab" argument, banded matrix format in C++

[ozarbit]

Hello, I am now playing with the banded matrix form for the Ax=B system anddgbsv() driver function.

I have a 5 band matrix ( 2 upper + diagonal + 2lower )

I coulnd't find in the MKL manual the format expected, in C/C++, for banded matrices,

I couldn't even find the signature explicitely in the manual, but i found it in ...

the double* ab points to the banded matrix.

double ab[] = {
0., 0., 5., 6. , 7. , 8. , 9. , 10., 11., // most upper diagonal
0., 5., 6., 7. , 8. , 9. , 10., 11., 12., // upper diagonal
5., 6., 7., 8. , 9. , 10., 11., 12., 13., // diagonal
6., 7., 8., 9. , 10., 11., 12., 13., 0. , // lower diagonal
7., 8., 9., 10., 11., 12., 13., 0. , 0. }; // lowest diagonal

double b[] =
{ 0., 1., 2., 3., 4., 5., 6., 7. , 8. };

const MKL_INT n =sizeof(b)/sizeof(double);
MKL_INT nn = n;
MKL_INT nrhs =1;
MKL_INT kl =2;
MKL_INT ku =2;
MKL_INT ldab =5;
MKL_INT info;
MKL_INT ipiv;

LARGE_INTEGER t0, t1;
dgbsv( &nn, &kl, &ku, &nrhs, ab, &ldab, ipiv, b, &nn, &info );

the error is

MKL ERROR: Parameter 6 was incorrect on entry to DGBSV

is Parameter 6 ldab (value 5) ?or is it ipiv?

regards,

[TimP]

Quoting - ozarbit

const MKL_INT n =sizeof(b)/sizeof(double);
MKL_INT nn = n;
MKL_INT nrhs =1;
MKL_INT kl =2;
MKL_INT ku =2;
MKL_INT ldab =5;
MKL_INT info;
MKL_INT ipiv;

LARGE_INTEGER t0, t1;
dgbsv( &nn, &kl, &ku, &nrhs, ab, &ldab, ipiv, b, &nn, &info );

the error is

MKL ERROR: Parameter 6 was incorrect on entry to DGBSV

is Parameter 6 ldab (value 5) ?or is it ipiv?

I suggest referring to the reference source http://www.netlib.no/netlib/lapack/double/dgbsv.f

LDAB    (input) INTEGER
*          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.

so you have violated this documented constraint.

[ozarbit]

Quoting - tim18

I suggest referring to the reference source http://www.netlib.no/netlib/lapack/double/dgbsv.f

LDAB    (input) INTEGER
*          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.

so you have violated this documented constraint.

Hello, I fixed the value of LDAB.
However the function call crashes now.

I am very unsure about the ab argument

From the documentation

"Band storage: an m-by-n band matrix with kl sub-diagonals and ku superdiagonals is stored compactly in a two-dimensional array ab with kl+ku+1 rows and n columns. Columns of the matrix are stored in the corresponding columns of the array, and diagonals of the matrix are stored in rows of the array."

The compact storage is then a 5 rows 9 columns in my case.

However, the C function signature expects a double*.

What should the ab array layout be? I looked at netlib.no re any remarks for the C layout w/o success.

regards,

[TimP]

Quoting - ozarbit

Hello, I fixed the value of LDAB.
However the function call crashes now.

I am very unsure about the ab argument

From the documentation

"Band storage: an m-by-n band matrix with kl sub-diagonals and ku superdiagonals is stored compactly in a two-dimensional array ab with kl+ku+1 rows and n columns. Columns of the matrix are stored in the corresponding columns of the array, and diagonals of the matrix are stored in rows of the array."

The compact storage is then a 5 rows 9 columns in my case.

However, the C function signature expects a double*.

What should the ab array layout be? I looked at netlib.no re any remarks for the C layout w/o success.

regards,

No, it seems unlikely anyone has posted a C example of the usage. As you are probably aware, the C view of a 2d Fortran array is a 1d array which runs by columns of the Fortran array. By my reading of the diagram in the source code, it begins with the first column, the first 3 elements of which are "unused," then runs down the second column, and so forth. ending with the last column, with 2 "unused" elements at the end, in that example. This agrees with my reading of your text quote.
You probably should avoid attempting to interpret the resulting LU decomposition, if at all possible.
While this organization is backwards from the usage I remember in IMSL from years ago, it is better for data locality.

[ozarbit]

Quoting - tim18

No, it seems unlikely anyone has posted a C example of the usage. As you are probably aware, the C view of a 2d Fortran array is a 1d array which runs by columns of the Fortran array. By my reading of the diagram in the source code, it begins with the first column, the first 3 elements of which are "unused," then runs down the second column, and so forth. ending with the last column, with 2 "unused" elements at the end, in that example. This agrees with my reading of your text quote.
You probably should avoid attempting to interpret the resulting LU decomposition, if at all possible.
While this organization is backwards from the usage I remember in IMSL from years ago, it is better for data locality.

I tested on Excel. That initial matrix was singular, I changed it to the following problem.

| 10 5 5 0 0 0 0 0 0 |
| 1 5 56 6 0 0 0 0 0 |
| 7 7 7 7 7 0 0 0 0 |
| 0 8 8 8 8 8 0 0 0 |
| 0 0 9 9 10 9 9 0 0 | = A
| 0 0 0 10 10 10 10 10 0 |
| 0 0 0 0 11 11 11 11 11 |
| 0 0 0 0 0 12 12 12 12 |
| 0 0 0 0 0 0 13 13 13 |

This a band matrix, 2 upper diagonals and 2 lower diagonals. This actually has determinant !=0 and

A*x = b

where

b = 0

1

2

3

4

5

6

7

8

has a solution.

Now, the signature in C of the driver function for Ax=b where A is banded is

dgbsv( MKL_INT *n, MKL_INT *kl, MKL_INT *ku, MKL_INT *nrhs, double *ab, MKL_INT *ldab, MKL_INT *ipiv, double *b, MKL_INT *ldb, MKL_INT *info );

MKL_INT n = 9;

MKL_INT kl = 2;
MKL_INT ku = 2;
MKL_INT nrhs =1;
MKL_INT ldab = 2*kl + ku +1 = 7;
MKL_INT ipiv;
MKL_INT ldb = n;

MKL_INT info;

In fortran, the compact storage of A in ab would be:

x x 5 6 7 8 9 10 11

x 5 56 7 8 9 10 11 12

10 5 7 8 9 10 11 12 13

1 7 8 9 10 11 12 13 x

7 8 9 10 11 12 13 x x

But, from the docs, "The second dimension of ab must be at least max(1, n). ", i.e. 9

"The first dimension of the array ab. (ldab 2kl + ku +1)", i.e. 7

The 1stdimis the number of rows in Fortran, so 7 rows

The 2nd dim is the number of columns in Fortran so 9 columns

So i change the compact storage to include 2 additional rows (probably required for the LU decomposition)

x x 5 6 7 8 9 10 11
x 5 56 7 8 9 10 11 12
10 5 7 8 9 10 11 12 13
1 7 8 9 10 11 12 13 x
7 8 9 10 11 12 13 x x

x x x x x x x x x

x x x x x x x x x

Then I transform this 2d array to 1d array in C, based on your reply" which runs by columns of the Fortran array",

to give me

double ab[] = {
0. , 0. , 10., 1. , 7. , 0., 0.,
0. , 5. , 5. , 7. , 8. , 0., 0.,
5. , 56., 7. , 8. , 9. , 0., 0.,
6. , 7. , 8. , 9. , 10., 0., 0.,
7. , 8. , 9. , 10., 11., 0., 0.,
8. , 9. , 10., 11., 12., 0., 0.,
9. , 10., 11., 12., 13., 0., 0.,
10., 11., 12., 13., 0. , 0., 0.,
11., 12., 13., 0. , 0. , 0., 0.
};

and

double b[] = { 0., 1., 2., 3., 4., 5., 6., 7. , 8. };

callingdgbsv( &nn, &kl, &ku, &nrhs, ab, &ldab, ipiv, b, &nn, &info );

returns info =8

which means U(8,8) is 0 according to docs, however matlab and excel concord in the solution.

I deduce that my transposition to the C array is wrong.

thank you for any advice,

best regards,

[TimP]

In the dgbsv source code listing comments, the expansion rows are shown above the entry data rows (preceding, in C), opposite to what you show here. I think your correspondence between Fortran and C array layout is OK.

[ozarbit]

Quoting - tim18

In the dgbsv source code listing comments, the expansion rows are shown above the entry data rows (preceding, in C), opposite to what you show here. I think your correspondence between Fortran and C array layout is OK.

finally, that was it.... thanks very much.

I think the MKL manual should definitely describe the placement of the expansion rows.

thanks again

[Michael_C_Intel4]

Just locate the matrix differently:

double ab[] = {
0., 0., 0. , 0. , 10., 1. , 7. ,
0., 0., 0. , 5. , 5. , 7. , 8. ,
0., 0., 5. , 56., 7. , 8. , 9. ,
0., 0., 6. , 7. , 8. , 9. , 10.,
0., 0., 7. , 8. , 9. , 10., 11.,
0., 0., 8. , 9. , 10., 11., 12.,
0., 0., 9. , 10., 11., 12., 13.,
0., 0., 10., 11., 12., 13., 0. ,
0., 0., 11., 12., 13., 0. , 0.
};

Extra kl elements are added to the upper factor U during factorization, therefore the space should be left above the matrix, not below.

There's a description of the banded matrix arguments in Appendix B. Matrix Arguments of MKL manual mklman.pdf, you may refer to it.

Do you feel that the MKL documentation is ambiguous or not sufficient in this case?

Michael.

[ozarbit]

Quoting - Michael Chuvelev (Intel)

Just locate the matrix differently:

double ab[] = {
0., 0., 0. , 0. , 10., 1. , 7. ,
0., 0., 0. , 5. , 5. , 7. , 8. ,
0., 0., 5. , 56., 7. , 8. , 9. ,
0., 0., 6. , 7. , 8. , 9. , 10.,
0., 0., 7. , 8. , 9. , 10., 11.,
0., 0., 8. , 9. , 10., 11., 12.,
0., 0., 9. , 10., 11., 12., 13.,
0., 0., 10., 11., 12., 13., 0. ,
0., 0., 11., 12., 13., 0. , 0.
};

Extra kl elements are added to the upper factor U during factorization, therefore the space should be left above the matrix, not below.

There's a description of the banded matrix arguments in Appendix B. Matrix Arguments of MKL manual mklman.pdf, you may refer to it.

Do you feel that the MKL documentation is ambiguous or not sufficient in this case?

Michael.

Thanks Michael.

As a matter of a fact, the example in Appendix B Matrix arguments. The expansion rows are visible in the example. It is clear. I apologize.

The Fortan to C translation for arrays is also probably mentioned somewhere in the manual, but if it isn't, it may be useful to list as well. Thanks again,

regards,