DGBTRF fails for large matrices

[LucGoe]

Dear reader

To compute the inverse of a square, non-symmetric, pentadiagonal matrix [A] of the following form

Code: Select all

[A] =
[  1 a12 a13   0   0     
 a21   1 a23 a24   0   
 a31 a32   1 a34 a35
   0 a42 a43   1 a45
   0   0 a53 a54   1]


I compactly store [A] as [AB]:

Code: Select all

[AB] =
[  0   0   0   0   0
   0   0   0   0   0
   0   0 a13 a24 a35
   0 a12 a23 a34 a45
   1   1   1   1   1
 a21 a32 a43 a54   0
 a31 a42 a53   0   0]


and pass it to DGBTRF to find the LU decomposition (for the exemplary matrix above):

Code: Select all

M = 5
N = 5
KL = 2
KU = 2
LDAB = 2*KL+ KU+ 1
call DGBTRF(M, N, KL, KU, AB, LDAB, IPIV, INFO)


This works as expected. If however the dimensions of [A] exceed 9x9, DGBTRF returns a value of 1 for INFO,
i.e. a singular factor U in the LU decomposition.
Using MATLAB's lu(A) gives sensible results for matrices of any size.

I would greatly appreciate your thoughts on where I am going wrong with LAPACK.
Thank you!
Luc

[LucGoe]

A side-note:

Trying to compute the inverse of [A] directly with the routine DGBSV() returns an INFO value
of 1 for matrices of any size.

[rodney]

The banded matrix below has two extra upper bands. You should remove these first two bands of zeros. Otherwise, you will need to set KU=4 instead of 2, as well as changing the values of M=5 and LDAB=5 to 7.

Rodney

[/code]
I compactly store [A] as [AB]:

Code: Select all

[AB] =
[  0   0   0   0   0
   0   0   0   0   0
   0   0 a13 a24 a35
   0 a12 a23 a34 a45
   1   1   1   1   1
 a21 a32 a43 a54   0
 a31 a42 a53   0   0]


and pass it to DGBTRF to find the LU decomposition (for the exemplary matrix above):

Code: Select all

M = 5
N = 5
KL = 2
KU = 2
LDAB = 2*KL+ KU+ 1
call DGBTRF(M, N, KL, KU, AB, LDAB, IPIV, INFO)


This works as expected. If however the dimensions of [A] exceed 9x9, DGBTRF returns a value of 1 for INFO,
i.e. a singular factor U in the LU decomposition.
Using MATLAB's lu(A) gives sensible results for matrices of any size.

I would greatly appreciate your thoughts on where I am going wrong with LAPACK.
Thank you!
Luc

[LucGoe]

Dear Rodney

To my understanding the two rows of zeros in [AB] are consistent with the LAPACK
band storage format as described in the documentation:

Note: when a band matrix is supplied for LU factorization, space must be allowed to store an
additional kl superdiagonals, generated by fill-in as a result of row interchanges. This means
that the matrix is stored according to the above scheme, but with kl + ku superdiagonals.

Source: [ http://www.netlib.org/lapack/lug/node124.html ]

Regards
Luc

[rodney]

Luc,

You are correct; my mistake, I had forgotten about the extra storage needed for the factorization. Could you post an example of a banded matrix with n>9 for which xGBTRF fails?

Rodney

[LucGoe]

Dear Rodney

I attached a 10x10 matrix [A] in the text-file below. Apart from the first three and last three rows
the entries in the band do not change from row to row.
For matrices of any larger size the only difference is that the interior band is longer, the entries
themselves stay the same (the matrix is part of a finite difference scheme).

Thank you!
Regards
Luc

[rodney]

Dear Luc,

I converted your general matrix to band storage form with M=N=10, KL=KU=2, LDAB=7, which is attached.

and successfully called DGBTRF (info=0). I compared the output to the general matrix form, which I got from calling DGETRF with your original matrix, and the answers were identical (the output of DGBTRF in band format, and DGETRF in general format).

I am not sure what could be happening where you are getting an info=1. I tried a few different versions of LAPACK and all were okay. What version of LAPACK are you using?

Rodney

[LucGoe]

Dear Rodney

THANK YOU for taking the time to look into this!

After examining your matrix in compact form (it only differs from the one in my code
in the data type/format) I made a few tests. The conclusion:
DGBTRF fails if the matrix passed to it is of type DOUBLE PRECISION. It runs successfully
if it is of type REAL.

This is quite peculiar - I will have a second look at the installation I'm using (LAPACK 3.4.2).

Regards
Luc

[rodney]

Dear Luc,

I did my tests using double precision; the data file I attached was not the one I used, which had full precision for double. I think that your problem should work in either single or double precision

Rodney