bitofbully
12-11-2008, 08:02 AM
Hi,
I am looking to solve a simple problem using numerical recipes code for the first time. I am a long time programmer but never before used Numerical Recipes - i was looking forward to it but it has not been a seamless use of the product so far ... hopefully someone can explain to me some issues I am having
Some important issues are the following
I am basically trying to test a simple direct solver for Ax = b where A represents the Laplacian in 2D so d2u/ dx2 + d2u/dy2 and I am using an LU type decomposition for this where i use packed storage and the 2D laplacian has obviously nice structure to it with large gap of zeros between L_x and L_y etc. For now, let us take a square domain with Nx = 10 and Ny 10 (*** but I will want this to work with anisotropic mesh so maybe Nx = 100 Ny = 30 etc )
The matrix for this operator is a banded matrix with a large bandwidth between dimension x and y of size Ny
(1) I am using row-major storage
(2) I am using C and C++ so array indexing starts at 0 not 1
(3) I am using full matrices that incorporate boundary conditions so the example matrix for the 1D second derivative with Dirichlet conditions is
1 0 0 0 ........0
-1 2 -1 0....... 0
0 -1 2 -1 0 ...... 0
-
-
-
0 ...........-1 2 -1
0 .............0 0 1
not the, smaller, usual
2 -1 0 0 ....... 0
-
-
0 .............0 -1 2
This matrix formatiing allows to work purely with matrices and use differing boundary conditions easily and not use boundary vectors adding onto right-hand side.
Some sample code for using bandec( ) and banbks( ) for my problem is the following
int nx = num_rows;
int ny = num_cols;
int packed_rows = num_rows * num_cols;
int packed_cols = 5 // for the 2D laplacian etc i.e m1 + m2 + 1
// The solution vector
double *u = new double [ packed_rows ];
double *A = new double [ packed_rows * packed_cols ];
double *L = new double [ packed_rows * num_cols ]; // matrices as according to chapter 2.4 so L is different to A. From bandec( ) comments
/* set_up and initialise matrix A using my own routine */
initialise_matrix( A, nx, ny );
bandec( A, packed_rows, num_cols, num_cols, L, indx, d );
banbks( A, packed_rows, num_cols, num_cols, L, indx, u );
where u is the solution vector and packed_rows is just a int that specifies how many rows there are in packed storage format i.e the diagonal of the original matrix.
Here is the sample output of the packed storage of the 2D laplacian for nx = 5 ny = 5 with boundary conditions built into the matrix
----------------------------------------------------------------------------------------------------------------------
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
------------------------------------------------------------------------------------------------------------------------
SO .... my questions are the following
1 - Can bandec( ) and banbks( ) handle this type of matrix ... which is just a simple 2D laplacian with banded storage according to 2.4 from the book and has boundary condiitons built into the matrix ???
Because my output for A (which according to the recipes book contains U ) and L are the following
Matrix A
1e-20 0 0 0 0 0 0 0 0 0 0
1e-20 1 0 0 0 0 0 0 0 0 0
1e+20 0 0 -1 4 -1 0 0 0 -1 0
1 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0
-1 4 -1 0 0 -1e-20 -1 0 0 0 0
-4 1 0 -1 4 0 0 0 0 -1 0
3.75 -1 0.25 -1 0 -1 0 0 0.25 0 0
1 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0
-2.067 4 -1.067 0 0 0.2667 -1 0 0 0 0
-1.935 0.5161 0 -1 3.871 -0.5161 0 0 0 -1
(where I have only shown half the matrix A due to storgage in this message)
Obviously this does not look right ? As A after bandec( ) should now be U which is upper diagonal matrix I also had to change by hand the array indexing from 1 to 0 and wondering if this is the cause of the issue ? (i still dont know why a Num Recipes in C book uses [1] to [N] and not [0] to [N-1] )
2 - Is there any need for banmul( ) if we want to solve Ax = b using LU decomposition where A = 2D Laplacian ?
3 - Why in Chapter 2 section 4 in the bandec( ) and banbks( ) algorithms is the indexing and loops from i = 1 ... i <= N .... surely this is numerical recipes in C not matlab etc so all indexing is from 0 to N-1 ? Here, I would echo the comment from the thread where he/she discusses changing indices from 1 to 0 for arrays etc http://www.nr.com/forum/showthread.php?t=1541&highlight=boundary+conditions
4 - Why are all matrices stored as literally a double array [ ][ ] ... is it not better to store them as 1 long array [ ] and stride through it since we just nees to say [ i*num_cols + j ] to get [i i ][ j ]? this will help avoid cache missing especially for large matrices ?
5 - Is bandec( ) and banbks( ) only useable for specific matirces ? i.e square, symmetric, interior matrices so do not include boundary conditions in the matrices ? if so it would be useful to either state the criteria in a reply or in the book itself.
6 - Isnt the 2D laplacian and obvious test-case for packed/banded storage ? Obviously a pentadiagonal matrix is another test case but a 2D laplacian would be more useful as it uses both the large bandwidth of zeros and the packed storage ..... i dont know why a simple example that was shown in 2.4 was used as it far less realistic than pent or laplace ?
Note: I have been coding for a long time and write most of my code myself but am now looking to make use of Numerical Recipes IF I can get this test code to work and like its functionality ( hence my maybe naive question 2 regarding indexing as this might be explained in the book somewhere but i havent seen it )
I hope my post is not too critical !! but I was just trying to point out some things and do hope I get a reply ... even if it is critical of my question/points.
Thanks for any help offered
Bbully
I am looking to solve a simple problem using numerical recipes code for the first time. I am a long time programmer but never before used Numerical Recipes - i was looking forward to it but it has not been a seamless use of the product so far ... hopefully someone can explain to me some issues I am having
Some important issues are the following
I am basically trying to test a simple direct solver for Ax = b where A represents the Laplacian in 2D so d2u/ dx2 + d2u/dy2 and I am using an LU type decomposition for this where i use packed storage and the 2D laplacian has obviously nice structure to it with large gap of zeros between L_x and L_y etc. For now, let us take a square domain with Nx = 10 and Ny 10 (*** but I will want this to work with anisotropic mesh so maybe Nx = 100 Ny = 30 etc )
The matrix for this operator is a banded matrix with a large bandwidth between dimension x and y of size Ny
(1) I am using row-major storage
(2) I am using C and C++ so array indexing starts at 0 not 1
(3) I am using full matrices that incorporate boundary conditions so the example matrix for the 1D second derivative with Dirichlet conditions is
1 0 0 0 ........0
-1 2 -1 0....... 0
0 -1 2 -1 0 ...... 0
-
-
-
0 ...........-1 2 -1
0 .............0 0 1
not the, smaller, usual
2 -1 0 0 ....... 0
-
-
0 .............0 -1 2
This matrix formatiing allows to work purely with matrices and use differing boundary conditions easily and not use boundary vectors adding onto right-hand side.
Some sample code for using bandec( ) and banbks( ) for my problem is the following
int nx = num_rows;
int ny = num_cols;
int packed_rows = num_rows * num_cols;
int packed_cols = 5 // for the 2D laplacian etc i.e m1 + m2 + 1
// The solution vector
double *u = new double [ packed_rows ];
double *A = new double [ packed_rows * packed_cols ];
double *L = new double [ packed_rows * num_cols ]; // matrices as according to chapter 2.4 so L is different to A. From bandec( ) comments
/* set_up and initialise matrix A using my own routine */
initialise_matrix( A, nx, ny );
bandec( A, packed_rows, num_cols, num_cols, L, indx, d );
banbks( A, packed_rows, num_cols, num_cols, L, indx, u );
where u is the solution vector and packed_rows is just a int that specifies how many rows there are in packed storage format i.e the diagonal of the original matrix.
Here is the sample output of the packed storage of the 2D laplacian for nx = 5 ny = 5 with boundary conditions built into the matrix
----------------------------------------------------------------------------------------------------------------------
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
-1 0 0 0 -1 4 -1 0 0 0 -1
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
------------------------------------------------------------------------------------------------------------------------
SO .... my questions are the following
1 - Can bandec( ) and banbks( ) handle this type of matrix ... which is just a simple 2D laplacian with banded storage according to 2.4 from the book and has boundary condiitons built into the matrix ???
Because my output for A (which according to the recipes book contains U ) and L are the following
Matrix A
1e-20 0 0 0 0 0 0 0 0 0 0
1e-20 1 0 0 0 0 0 0 0 0 0
1e+20 0 0 -1 4 -1 0 0 0 -1 0
1 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0
-1 4 -1 0 0 -1e-20 -1 0 0 0 0
-4 1 0 -1 4 0 0 0 0 -1 0
3.75 -1 0.25 -1 0 -1 0 0 0.25 0 0
1 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0
-2.067 4 -1.067 0 0 0.2667 -1 0 0 0 0
-1.935 0.5161 0 -1 3.871 -0.5161 0 0 0 -1
(where I have only shown half the matrix A due to storgage in this message)
Obviously this does not look right ? As A after bandec( ) should now be U which is upper diagonal matrix I also had to change by hand the array indexing from 1 to 0 and wondering if this is the cause of the issue ? (i still dont know why a Num Recipes in C book uses [1] to [N] and not [0] to [N-1] )
2 - Is there any need for banmul( ) if we want to solve Ax = b using LU decomposition where A = 2D Laplacian ?
3 - Why in Chapter 2 section 4 in the bandec( ) and banbks( ) algorithms is the indexing and loops from i = 1 ... i <= N .... surely this is numerical recipes in C not matlab etc so all indexing is from 0 to N-1 ? Here, I would echo the comment from the thread where he/she discusses changing indices from 1 to 0 for arrays etc http://www.nr.com/forum/showthread.php?t=1541&highlight=boundary+conditions
4 - Why are all matrices stored as literally a double array [ ][ ] ... is it not better to store them as 1 long array [ ] and stride through it since we just nees to say [ i*num_cols + j ] to get [i i ][ j ]? this will help avoid cache missing especially for large matrices ?
5 - Is bandec( ) and banbks( ) only useable for specific matirces ? i.e square, symmetric, interior matrices so do not include boundary conditions in the matrices ? if so it would be useful to either state the criteria in a reply or in the book itself.
6 - Isnt the 2D laplacian and obvious test-case for packed/banded storage ? Obviously a pentadiagonal matrix is another test case but a 2D laplacian would be more useful as it uses both the large bandwidth of zeros and the packed storage ..... i dont know why a simple example that was shown in 2.4 was used as it far less realistic than pent or laplace ?
Note: I have been coding for a long time and write most of my code myself but am now looking to make use of Numerical Recipes IF I can get this test code to work and like its functionality ( hence my maybe naive question 2 regarding indexing as this might be explained in the book somewhere but i havent seen it )
I hope my post is not too critical !! but I was just trying to point out some things and do hope I get a reply ... even if it is critical of my question/points.
Thanks for any help offered
Bbully