Solve a banded matrix system of equations

2
votes

I need to solve a 2D Poisson equation, that is, a system of equations in the for AX=B where A is an n-by-n matrix and B is a n-by-1 vector. Being A a discretization matrix for the 2D Poisson problem, I know that only 5 diagonals will be not null. Lapack doesn't provide functions to solve this particular problem, but it has functions for solving banded matrix system of equations, namely DGBTRF (for LU factorization) and DGBTRS. Now, the 5 diagonals are: the main diagonal, the first diagonals above and below the main and two diagonals above and below by m diagonals wrt the main diagonal. After reading the lapack documentation about band storage, I learned that I have to create a (3*m+1)-by-n matrix to store A in band storage format, let's call this matrix AB. Now the questions:

1) what is the difference between dgbtrs and dgbtrs_? Intel MKL provides both but I can't understand why

2) dgbtrf requires the band storage matrix to be an array. Should I linearize AB by rows or by columns?

3) is this the correct way to call the two functions?

int n, m;
double *AB;
/*... fill n, m, AB, with appropriate numbers */

int *pivots;
int nrows = 3 * m + 1, info, rhs = 1;
dgbtrf_(&n, &n, &m, &m, AB, &nrows, pivots, &info);
char trans = 'N';
dgbtrs_(&trans, &n, &m, &m, &rhs, AB, &nrows, pivots, B, &n, &info);
2
csystemlapackequationspoisson

2 Answers

0
votes

  1. It also provides DGBTRS and DGBTRS_. Those are fortran administrativa that you should not care about. Just call dgbtrs (reason is that on some architectures, fortran routine names have underscore appended, on other not, and names may be either upper or lower case -- Intel MKL #defines the right one to dgbtrs).

  2. LAPACK routines expects column major matrices (ie. Fortran style): store columns one after the others. The banded storage you must use is not hard : http://www.netlib.org/lapack/lug/node124.html.

  3. It seems good to me, but please try it on small problems beforehand (always a good idea by the way). Also make sure you handle non-zero info (this is the way errors are reported).

Better style is to use MKL_INT instead of plain int, this is a typedef to the right type (may be different on some architectures).

Also make sure to allocate memory for pivots before calling dgbtrf.

0
votes

This might be off topic. But for Poisson equation, FFT based solution is much faster. Just do 2D FFT of your potential field, divided by -(k^2+lambda^2) then do IFFT. lambda is a small number to avoid divergence for k=0. The 5-diagonal equation is a band-limited approximation of the Poisson equation, which approximate the differential operator by finite difference.

http://en.wikipedia.org/wiki/Screened_Poisson_equation