DGBTRF fails for large matrices

Open discussion regarding features, bugs, issues, vendors, etc.

DGBTRF fails for large matrices

Postby LucGoe » Wed Apr 03, 2013 5:20 am

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
 
Posts: 5
Joined: Wed Apr 03, 2013 2:38 am

Re: DGBTRF fails for large matrices

Postby LucGoe » Fri Apr 05, 2013 1:21 am

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.
LucGoe
 
Posts: 5
Joined: Wed Apr 03, 2013 2:38 am

Re: DGBTRF fails for large matrices

Postby rodney » Fri Apr 05, 2013 4:20 pm

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

LucGoe wrote:[/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
rodney
 
Posts: 48
Joined: Thu Feb 10, 2011 8:20 pm
Location: Colorado College

Re: DGBTRF fails for large matrices

Postby LucGoe » Sat Apr 06, 2013 3:57 am

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
LucGoe
 
Posts: 5
Joined: Wed Apr 03, 2013 2:38 am

Re: DGBTRF fails for large matrices

Postby rodney » Sat Apr 06, 2013 8:45 am

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
rodney
 
Posts: 48
Joined: Thu Feb 10, 2011 8:20 pm
Location: Colorado College

Re: DGBTRF fails for large matrices

Postby LucGoe » Sun Apr 07, 2013 12:08 pm

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
Attachments
matrix_original.txt
Pentadiagonal matrix [A]
(2.94 KiB) Downloaded 29 times
LucGoe
 
Posts: 5
Joined: Wed Apr 03, 2013 2:38 am

Re: DGBTRF fails for large matrices

Postby rodney » Sun Apr 07, 2013 5:27 pm

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
Attachments
B.txt
banded input matrix
(1.02 KiB) Downloaded 26 times
rodney
 
Posts: 48
Joined: Thu Feb 10, 2011 8:20 pm
Location: Colorado College

Re: DGBTRF fails for large matrices

Postby LucGoe » Sun Apr 07, 2013 6:13 pm

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
LucGoe
 
Posts: 5
Joined: Wed Apr 03, 2013 2:38 am

Re: DGBTRF fails for large matrices

Postby rodney » Sun Apr 07, 2013 9:25 pm

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
rodney
 
Posts: 48
Joined: Thu Feb 10, 2011 8:20 pm
Location: Colorado College


Return to User Discussion

Who is online

Users browsing this forum: Google [Bot] and 2 guests