Solving a banded system of equations from C using LAPACK's DGBSV
I have tried to search for a similar thread related to this topic, but nobody seems to care for banded systems (...)
I am interested in solving a banded matrix using LAPACK/ScaLAPACK from a C code. First, I want to achieve a sequential solution with LAPACK, before attempting anything with ScaLAPACK. Problem: The row-major/col-major difference between both languages seems to be affecting my solution process. Here's the system I intend to solve: Code:
[ 1, 0, 0, 0, 0, 0] [1] Code:
int rr = 6; // Rank. Code:
[ejspeiro@node01 lapack-ex02]$ make runs Any help is more than welcome. Thanks in advanced! |
Hi.
I solved your system of equations using Maxima, and it gives: Code:
matrix([1.0],[0.99999999999998],[1.000000000000018],[1.0],[0.99999999999999],[0.0]) It is quite close to what you have. I believe you actually have exactly the same solution, which is rounded by printf. Try something like Code:
printf("%.16f\n", rhs[ii]); |
Dear firstfire:
I must say, I am surprised... I have solved the system many times, in both MATLAB and Maple. They both agree on this: Code:
[1.0000000, 1.027346753, 1.105215007, 1.244831762, 1.459918908, 1.767247642] https://dl.dropbox.com/u/5432016/misc/plot.png Which, agrees with the actual analytical solution of the Partial Differential Equation. Furthermore, the value info being equal to 1, implies that the factorization was achieved, but a solution couldn't have been performed. Thank you very much for your reply. |
Quote:
|
Hi.
I figured out what is wrong with your code. First, your upper diagonal is entered incorrectly ([0,0,0,0,22.5] instead of [0, 22.5, 22.5, 22.5, 22.5]). Second, you must transpose the matrix before passing it to fortran wrapper routine. Here is the code with my modifications: Code:
#include <stdio.h> Code:
$ gcc test.c -llapack && ./a.out |
Excellent! Yes! I had come into that solution of transposing the "stack of arrays". However, I wonder, how would the transposition occur if the resulting CORE matrix were not to be squared. Take for example:
Code:
Core of the matrix: Code:
for (ii = 0; ii < in->co_num_rows; ii++) { Thanks! |
Hi.
Can you provide complete example which yields an error and the error message? Your code looks good to me and should work equally well on non-square matrix, provided that arrays in and out are of correct dimensions. |
Sight u.u' you are right... it would be helpful... let us see... I have created an API for this data type, so I apologize for any incorrect programming practice... although, hopefully, I will improve my code writing by means of your suggestions! :)
My project code is comprise of 4 files total. I am still working on these files! Feel free to highlight any error, but hopefully, I would have also spotted it ;) File 1: https://dl.dropbox.com/u/5432016/banded/cmtk.h File 2: https://dl.dropbox.com/u/5432016/banded/blogs.c File 3: https://dl.dropbox.com/u/5432016/banded/makefile (Optional if library paths are known?) File 4: https://dl.dropbox.com/u/5432016/banded/blogs.in Once everything is in one folder, this should build and execute: Code:
$ make Thanks! |
Hi.
Here are my observations: 1) To fix segfault in Transpose() you should swap dimensions of the transposed matrix (but see item 4 below): Code:
out->co_num_rows = out->fa_rank; //out->fa_bandwidth; Code:
for (ii = 0; ii < in->co_num_rows; ii++) { 2) This line Code:
rhs = (CMTK_Number *) calloc(bl->fa_rank, sizeof(CMTK_Number *)); 3) In my (not very big) experience 2D arrays (arrays of arrays) are not very convenient for representing two dimensional (and more) matrices in C. They are PITA to allocate, deallocate and carry around. I would use linear representation instead. Moreover, LAPACK functions also seem to prefer linear representation. 4) If Transpose() operates on general band matrices, then it may do two a little bit different things: A) Transpose a band matrix (basically exchange its diagonals) B) Transpose the core of a band matrix (i.e. transpose a compact representation of a band matrix) to be able pass such matrix to LAPACK routine. It seems that in its current state your code does not make distinction between these two cases and even mix them a little bit. 5) It is a matter of taste, I know, but why not omit these redundant curly braces? :) |
Dear firstfire: Thanks! This has been an incredibly helpful post!
1) Done! Gratefully, I had already spotted that! 2) You are right! In fact, I have made the adequate changes throughout the code! 3) Again, you are right! I have changed the entire implementation to rely on a 1-d array! Actually I have updated the entire code: https://dl.dropbox.com/u/5432016/banded/cmtk.h https://dl.dropbox.com/u/5432016/banded/blogs.c 4) ... 4) ... 4) Sight! u.u' This one is still giving me problems :( I am not able to get a solution for cases in where I get a non-squared CORE matrix /* Notice that the FACADE matrices WILL always be squared, since we are solving for a PDE, but the cores can be non-squared. */ :( Say 1 PDE, step size of 0.1, as in https://dl.dropbox.com/u/5432016/banded/blogs.in I have grown aware of the differences between transposing the banded matrix and actually transposing its core, (thanks for your help with that), and I have developed these two methods (under debugging): Code:
/*! Code:
/*! Thanks in advanced! |
I believe the issue arises when accessing the values in the soon-to-be transposed matrices, since I can only do it within the scope of the band (because the core is entirely pre-allocated), and for some reason, there are cases in where I add in locations outside of the scope of the band.
May there be something funky with the way I define the limits of my cores? \m/ |
All right! I think I have solved all of my issues with this code by now! :D Thank you firstfire! I will call this thread as solved, post the final files:
https://dl.dropbox.com/u/5432016/banded/cmtk.h https://dl.dropbox.com/u/5432016/banded/blogs.c And if anybody has any extra suggestions as far as the implementation goes, I will gladly welcome them! My next step is to distribute the core matrices, in order to solve with ScaLAPACK. Thanks in advanced! \m/ |
Hi.
I'm glad to see you solved the problem! New questions are welcome, we'll try to help. |
All times are GMT -5. The time now is 01:09 PM. |