I am getting bad output from the band-diagonal LU solver, Bandec, that I coded myself from NR3 (pp. 59-60). When I compare the book's listing with that in the earlier Numerical Recipes in C, 2nd Edition, there are some confusing differences. For example, where NR3 has the expression:

for(j=m1-i;j<mm;j++), the earlier code reads:

for(j=m1+2-i;j<=mm;j++).

Part of the difference is due to the different method of indexing, but in addition, the two statements yield a different number of iterations (mm-m1+i for NR3, mm-m1+i-1 in the earlier book). Is this correct? There are other differences that seem odd, as well.

Has anyone been successful using the algorithm as written in NR3? Thanks for any help! -- JAW

...bad output...Has anyone been successful ...
Using unmodified code from the NR3 CD for banmul() and Bandec:

// Solve system of equations defined with banded matrix

#include "../code/nr3.h"
#include "../code/banded.h"

int main()
{
int i, j;

VecDoub x(7), b(7);
MatDoub a(7, 4);


a[0][0] = 0.0000000; // arbitrary: not used
a[0][1] = 0.0000000; // arbitrary: not used
a[0][2] = 0.5297000;
a[0][3] = 0.9304360;
a[1][0] = 0.0000000; // arbitrary: not used
a[1][1] = 0.0668422;
a[1][2] = 0.7226600;
a[1][3] = 0.6711490;
a[2][0] = 0.6316350;
a[2][1] = 0.8847070;
a[2][2] = 0.5194160;
a[2][3] = 0.6515190;
a[3][0] = 0.2624530;
a[3][1] = 0.7621980;
a[3][2] = 0.7533560;
a[3][3] = 0.9092080;
a[4][0] = 0.2727100;
a[4][1] = 0.8976560;
a[4][2] = 0.2749070;
a[4][3] = 0.5162920;
a[5][0] = 0.2470390;
a[5][1] = 0.4865170;
a[5][2] = 0.8461670;
a[5][3] = 0.8309650;
a[6][0] = 0.9910370;
a[6][1] = 0.6792960;
a[6][2] = 0.7664950;
a[6][3] = 0.0000000; // arbitrary: not used

x[0] = 0.4159990;
x[1] = 0.3835020;
x[2] = 0.3834160;
x[3] = 0.2377740;
x[4] = 0.0726859;
x[5] = 0.3592650;
x[6] = 0.0345721;

cout << "Banded matrix [a]:" << endl;
cout << fixed << setprecision(7);
for (i = 0; i < a.nrows(); i++) {
for (j = 0; j < a.ncols(); j++) {
cout << setw(12) << a[i][j];
}
cout << endl;
}
cout << endl;

cout << "Vector [x]:" << endl;
for (i = 0; i < x.size(); i++) {
cout << setw(18) << x[i] << endl;
}
cout << endl;

// multiply a times x, giving b
banmul(a, 2, 1, x, b);

cout << "After multiplying [a] times [x] this is [b]:" << endl;

for (i = 0; i < b.size(); i++) {
cout << setw(18) << b[i] << endl;
}
cout << endl;

//
// First, save original value of x
//
VecDoub xsave(x);
//
// The constructor does the decomposition
//
Bandec banded(a, 2, 1);

//
// Now solve for x
//
banded.solve(b, x);

cout << "After solving [a] times [x] = [b] with bandec.solve():" << endl;
cout << " Original Solved" << endl
<< " x x" << endl;
for (i = 0; i < x.size(); i++) {
cout << setw(12) << xsave[i] << setw(12) << x[i] << endl;
}

return 0;
}



Output


Banded matrix [a]:
0.0000000 0.0000000 0.5297000 0.9304360
0.0000000 0.0668422 0.7226600 0.6711490
0.6316350 0.8847070 0.5194160 0.6515190
0.2624530 0.7621980 0.7533560 0.9092080
0.2727100 0.8976560 0.2749070 0.5162920
0.2470390 0.4865170 0.8461670 0.8309650
0.9910370 0.6792960 0.7664950 0.0000000

Vector [x]:
0.4159990
0.3835020
0.3834160
0.2377740
0.0726859
0.3592650
0.0345721

After multiplying [a] times [x] this is [b]:
0.5771787
0.5622771
0.9561131
0.6381052
0.5234681
0.4268288
0.3425810

After solving [a] times [x] = [b] with bandec.solve():
Original Solved
x x
0.4159990 0.4159990
0.3835020 0.3835020
0.3834160 0.3834160
0.2377740 0.2377740
0.0726859 0.0726859
0.3592650 0.3592650
0.0345721 0.0345721


Does your output agree with mine for this test program?

Regards,

Dave

Thanks, Dave -

I ran your test data and got the following:

After solving [a] times [x] = [b] with bandec.solve():
Original Solved
x x
0.4159990 -0.3577917
0.3835020 1.0896333
0.3834160 -0.2998465
0.2377740 0.5738064
0.0726859 0.1632076
0.3592650 0.0877249
0.0345721 0.1581812

Since you ran with the provided code, it looks like that is OK. I may have made an error transcribing, but based on the listings I quoted in my first post, I wonder if there's a disconnect between the listing in the book and the implemented code?

Hi,

The code in the book is produced from the same file that is distributed as the source code. This is to guarantee no typos.The only time things might get out of sync is if there is a bug fix that is incorporated in a reprinting of the book and you have an older printing. But NR3 is still in its first printing.

All the best,
Saul Teukolsky

The code in the book is produced from the same file that is distributed as the source code...NR3 is still in its first printing.

I suspected as much, and am really (really) glad to hear it. Thank you for confirming it.

To the Original Poster: If you would actually print out the values of l, i, and j for the specific line you are questioning for this example (or go through it with pencil and paper) I think you might see something like:

From Version 2.1 bandec.c (with one-based arrays), in the loop where it is setting a[i][j-l] = a[i][j]

l = 2, i = 1, j = 3
l = 2, i = 1, j = 4
l = 1, i = 2, j = 2
l = 1, i = 2, j = 3
l = 1, i = 2, j = 4


From Version 3 banded.h (with zero-based arrays), in the loop where it is setting au[i][j-l] = au[i][j]

l = 2, i = 0, j = 2
l = 2, i = 0, j = 3
l = 1, i = 1, j = 1
l = 1, i = 1, j = 2
l = 1, i = 1, j = 3



So you see it is actually going through the loop the same number of times and setting corresponding elements in the respective arrays.

Now, if the test program gave you the same answer for for the b vector (after multiplying [a] times the original [x]), but the solution from solve() is wrong then either

The decomposition in the constructor is wrong.

or

The back substitution in the solve() function is wrong.

Here's what I got from the al and au matrices after the constructor had decomposed the a matrix:

The test program was changed by adding:

Bandec banded(a, 2, 1);
cout << "Matrix al" << endl;
for (i = 0; i < banded.al.nrows(); i++) {
for (j = 0; j < banded.al.ncols(); j++) {
cout << "al[" << i << "][" << j << "] = "
<< setw(12) << banded.al[i][j] << endl;
}
}
cout << endl;
for (i = 0; i < banded.au.nrows(); i++) {
for (j = 0; j < banded.au.ncols(); j++) {
cout << "au[" << i << "][" << j << "] = "
<< setw(12) << banded.au[i][j] << endl;
}
}
cout << endl;


Output from this part:

al[0][0] = 0.1058241
al[0][1] = 0.8386172
al[1][0] = 0.2996732
al[1][1] = 0.4172300
al[2][0] = -0.8143695
al[2][1] = -0.4396814
al[3][0] = 0.5311212
al[3][1] = 0.3706459
al[4][0] = 0.3881024
al[4][1] = 0.7701015
al[5][0] = -0.4905923
al[5][1] = 0.0000000
al[6][0] = 0.0000000
al[6][1] = 0.0000000

au[0][0] = 0.6316350
au[0][1] = 0.8847070
au[0][2] = 0.5194160
au[0][3] = 0.6515190
au[1][0] = 0.6290367
au[1][1] = 0.6161823
au[1][2] = -0.0689464
au[1][3] = 0.0000000
au[2][0] = -0.6202445
au[2][1] = -0.5257137
au[2][2] = 0.0000000
au[2][3] = 0.0000000
au[3][0] = 0.6665095
au[3][1] = 0.2749070
au[3][2] = 0.5162920
au[3][3] = 0.0000000
au[4][0] = 0.9910370
au[4][1] = 0.6792960
au[4][2] = 0.7664950
au[4][3] = 0.0000000
au[5][0] = -0.7973405
au[5][1] = -0.5902789
au[5][2] = 0.0000000
au[5][3] = 0.0000000
au[6][0] = 0.2439002
au[6][1] = 0.0000000
au[6][2] = 0.0000000
au[6][3] = 0.0000000


If your decomposition is the same, then look at the back substitution in solve().

Regards,

Dave

Footnote: Note that for this example the following elements of bandec.al are not initialized and they are not used (so their values on the printout are irrelevant):

Ignore values shown for al[5][1], al[6][0], al[6][1]

Thanks, Dave, for the info. It helped me isolate my error. Seems I wrote "1" for "l" when copying the listing. BTW, shouldn't there be a law against mixing these in code?

...It helped me isolate my error... Great!!!
...shouldn't there be a law ...
I would vote for it. Punishment would be summary execution by public impalement.

Regards,

Dave

C'mon, guys. If we had been publicly executed for this, we could never have written the Third Edition. But seriously, we don't really intend that people should be typing in the code from the book, even though it is (barely) allowed by the rules on page xix, where it also says:
..experience has shown that virtually all reported bugs in such cases are typing errors! (italics in original)

You're supposed to read the code in the book for the ideas behind it. If you want to execute the code, you're supposed to buy the code download :)