# C++ code Thomas algorithm for solving a pentadiagonal Ax=b [closed]

I am looking to solve $$Ax=b$$ for $$x$$ where $$A$$ is pentadiagonal square matrix (elements on the upper and lower diagonals can however equal to zero) and $$x$$, $$b$$ two vectors of the same size.

The dimension of the problem can be large. $$A$$ can be 1000x1000 or more. I couldn't find a C++ implementation of a Thomas like algorithm ([Thomas Algorithm][1]). Has anyone come across one?

EDIT :

I am running this test and the solution does not equal to the product $$Ax$$ :

#include <iostream>

#include <cmath>
//#include <vector>
//#include "nr.h"   // this is a header provided with the book
//#include "my_matrix.h"    // substitute nr.h with your custom matrix class

using namespace std;

void bandec(double** a, const int m1, const int m2, double** al, int* indx, double& d, int dim)
// a is input matrix A_c and upper triangular matrix on output, m1 m2 dimensions of
//subdiagonals, al lower triangular matrix on output, indx records the row
//permutations, d = 1 or -1 depending on whether row interchanges are even or odd
{
const double TINY = 1.e-20;

int n = dim;    // custom function from your "my_matrix.h" class
int mm = m1 + 1 + m2;
int l = m1;
int i = 0;
for (i = 0; i < m1; i++) {
for (int j = m1 - i; j < mm; j++)
a[i][j-l] = a[i][j];
l--;
for (int j = mm - l - 1; j < mm; j++)
a[i][j] = 0.0;
}
d = 1.0;
l = m1;
double dum = 0;
int k = 0;
// finds pivot for each row in a
for (k = 0; k < n; k++) {
dum = a[k][0];
i = k;
if (l < n)
l++;
for (int j = k + 1; j < l; j++) {
if (fabs(a[j][0]) > fabs(dum)) {
dum = a[j][0];
i = j;
}
}

indx[k] = i + 1;
// in case input matrix a is singular
if (dum == 0.0)
a[k][0] = TINY;
if (i != k) {        // swap values in rows i and k
d = -d;
for (int j = 0; j < mm; j++) {
double aux = a[k][j];
a[k][j] = a[i][j];
a[i][j] = aux;
}
}
for (i = k + 1; i < l; i++) {   // perform elimination
dum = a[i][0] / a[k][0];
al[k][i-k-1] = dum;
for (int j = 1; j < mm; j++)
a[i][j-1] = a[i][j] - dum * a[k][j];
a[i][mm-1] = 0.0;
}
}
}

void banbks(double** a, const int m1, const int m2, double** al, int* indx, double* b, int dim)
// a, al , and indx inputs from bandec;
// m1 and m2 are the number of subdiagonals in starting system matrix A;
// b is the rhs of linear system Ax = b on input, solution vector x on output
{
int mm = m1 + 1 + m2;
int l = m1;
for (int k = 0; k < n; k++) {
int j = indx[k] - 1;
if (j != k) {
double aux = b[k];
b[k] = b[j];
b[j] = aux;
}
if (l < n)
l++;
for (j = k + 1; j < l; j++)
b[j] -= al[k][j-k-1] * b[k];
}
l = 1;
for (int i = n - 1; i >= 0; i--) {
double dum = b[i];
for (int k = 1; k < l; k++)
dum -= a[i][k] * b[k+i];
b[i] = dum / a[i][0];
if (l < mm)
l++;
}
}

void printC (double ** A, int locdim){
for (int ni = 0; ni < locdim; ++ni){
for (int nj = 0; nj < locdim; ++nj){
if (ni < 3 ){
if ( nj < 5)
std::cout << A[ni][nj] << " ";
else
std::cout << 0 << " ";
}
else if ( ni >= locdim - 3 ){
if ( nj >= locdim - 5 )
std::cout << A[ni][nj - ( locdim - 5 ) ] << " ";
else
std::cout << 0 << " ";
}
else{
if ( abs(nj-ni) <= 2 )
if (ni > nj)
std::cout << A[ni][ 2-( ni-nj) ] << " ";
else
std::cout << A[ni][ 2 + nj-ni ] << " ";
else
std::cout << 0 << " ";
}
}
std::cout << std::endl;
}
return;
}
void print (double ** A, int locdim){

for (int ni = 0; ni < locdim; ++ni){
for (int nj = 0; nj < locdim; ++nj){
if (abs (ni-nj)> 2) A[ni][nj] = 0.;
std::cout << A[ni][nj] << " ";
}
std::cout << std::endl;
}
return;
}

int main()
{
int locdim= 19;
double** A = new double*[locdim];
double** Ac = new double*[locdim];
double** al = new double*[locdim];
double* mult = new double[locdim];
double* b = new double[locdim];
int* indx = new int[locdim];

for (int ni = 0; ni < locdim; ++ni) {
A[ni] = new double[locdim];
Ac[ni] = new double[5];
al[ni] = new double[2];
mult[ni] = 0.;
}

A[0][0]=0.367456883431542;
A[1][0]=-0.160364379614644;
A[2][0]=0;
A[3][0]=0;
A[4][0]=0;
A[5][0]=0;
A[6][0]=0;
A[7][0]=0;
A[8][0]=0;
A[9][0]=0;
A[10][0]=0;
A[11][0]=0;
A[12][0]=0;
A[13][0]=0;
A[14][0]=0;
A[15][0]=0;
A[16][0]=0;
A[17][0]=0;
A[18][0]=0;
A[0][1]=-0.167129422166972;
A[1][1]=0.252958341387691;
A[2][1]=-0.0959599256411007;
A[3][1]=0.00950885432646034;
A[4][1]=0;
A[5][1]=0;
A[6][1]=0;
A[7][1]=0;
A[8][1]=0;
A[9][1]=0;
A[10][1]=0;
A[11][1]=0;
A[12][1]=0;
A[13][1]=0;
A[14][1]=0;
A[15][1]=0;
A[16][1]=0;
A[17][1]=0;
A[18][1]=0;
A[0][2]=0;
A[1][2]=-0.101034468246643;
A[2][2]=0.254833792976566;
A[3][2]=-0.0959599256411005;
A[4][2]=0;
A[5][2]=0;
A[6][2]=0;
A[7][2]=0;
A[8][2]=0;
A[9][2]=0;
A[10][2]=0;
A[11][2]=0;
A[12][2]=0;
A[13][2]=0;
A[14][2]=0;
A[15][2]=0;
A[16][2]=0;
A[17][2]=0;
A[18][2]=0;
A[0][3]=0;
A[1][3]=0.0107774899778459;
A[2][3]=-0.101034468246643;
A[3][3]=0.200836586771544;
A[4][3]=-0.10374883813066;
A[5][3]=0.0109311103569249;
A[6][3]=0;
A[7][3]=0;
A[8][3]=0;
A[9][3]=0;
A[10][3]=0;
A[11][3]=0;
A[12][3]=0;
A[13][3]=0;
A[14][3]=0;
A[15][3]=0;
A[16][3]=0;
A[17][3]=0;
A[18][3]=0;
A[0][4]=0;
A[1][4]=0;
A[2][4]=0;
A[3][4]=-0.101645109227889;
A[4][4]=0.263233346447371;
A[5][4]=-0.10374883813066;
A[6][4]=0;
A[7][4]=0;
A[8][4]=0;
A[9][4]=0;
A[10][4]=0;
A[11][4]=0;
A[12][4]=0;
A[13][4]=0;
A[14][4]=0;
A[15][4]=0;
A[16][4]=0;
A[17][4]=0;
A[18][4]=0;
A[0][5]=0;
A[1][5]=0;
A[2][5]=0;
A[3][5]=0.010405178131232;
A[4][5]=-0.101645109227888;
A[5][5]=0.210577303742067;
A[6][5]=-0.105506349217122;
A[7][5]=0.0110095501419037;
A[8][5]=0;
A[9][5]=0;
A[10][5]=0;
A[11][5]=0;
A[12][5]=0;
A[13][5]=0;
A[14][5]=0;
A[15][5]=0;
A[16][5]=0;
A[17][5]=0;
A[18][5]=0;
A[0][6]=0;
A[1][6]=0;
A[2][6]=0;
A[3][6]=0;
A[4][6]=0;
A[5][6]=-0.105662605927611;
A[6][6]=0.269008354233555;
A[7][6]=-0.105506349217121;
A[8][6]=0;
A[9][6]=0;
A[10][6]=0;
A[11][6]=0;
A[12][6]=0;
A[13][6]=0;
A[14][6]=0;
A[15][6]=0;
A[16][6]=0;
A[17][6]=0;
A[18][6]=0;
A[0][7]=0;
A[1][7]=0;
A[2][7]=0;
A[3][7]=0;
A[4][7]=0;
A[5][7]=0.0110486143195262;
A[6][7]=-0.105662605927611;
A[7][7]=0.216345569946779;
A[8][7]=-0.107440961096779;
A[9][7]=0.0110670374338971;
A[10][7]=0;
A[11][7]=0;
A[12][7]=0;
A[13][7]=0;
A[14][7]=0;
A[15][7]=0;
A[16][7]=0;
A[17][7]=0;
A[18][7]=0;
A[0][8]=0;
A[1][8]=0;
A[2][8]=0;
A[3][8]=0;
A[4][8]=0;
A[5][8]=0;
A[6][8]=0;
A[7][8]=-0.110546644894689;
A[8][8]=0.27582700508029;
A[9][8]=-0.107440961096779;
A[10][8]=0;
A[11][8]=0;
A[12][8]=0;
A[13][8]=0;
A[14][8]=0;
A[15][8]=0;
A[16][8]=0;
A[17][8]=0;
A[18][8]=0;
A[0][9]=0;
A[1][9]=0;
A[2][9]=0;
A[3][9]=0;
A[4][9]=0;
A[5][9]=0;
A[6][9]=0;
A[7][9]=0.0118434583833746;
A[8][9]=-0.110546644894688;
A[9][9]=0.218222694534443;
A[10][9]=-0.107440961096779;
A[11][9]=0.0110670374338972;
A[12][9]=0;
A[13][9]=0;
A[14][9]=0;
A[15][9]=0;
A[16][9]=0;
A[17][9]=0;
A[18][9]=0;
A[0][10]=0;
A[1][10]=0;
A[2][10]=0;
A[3][10]=0;
A[4][10]=0;
A[5][10]=0;
A[6][10]=0;
A[7][10]=0;
A[8][10]=0;
A[9][10]=-0.110546644894689;
A[10][10]=0.27582700508029;
A[11][10]=-0.107440961096779;
A[12][10]=0;
A[13][10]=0;
A[14][10]=0;
A[15][10]=0;
A[16][10]=0;
A[17][10]=0;
A[18][10]=0;
A[0][11]=0;
A[1][11]=0;
A[2][11]=0;
A[3][11]=0;
A[4][11]=0;
A[5][11]=0;
A[6][11]=0;
A[7][11]=0;
A[8][11]=0;
A[9][11]=0.0118434583833746;
A[10][11]=-0.110546644894689;
A[11][11]=0.218222694534443;
A[12][11]=-0.107440961096779;
A[13][11]=0.0110670374338972;
A[14][11]=0;
A[15][11]=0;
A[16][11]=0;
A[17][11]=0;
A[18][11]=0;
A[0][12]=0;
A[1][12]=0;
A[2][12]=0;
A[3][12]=0;
A[4][12]=0;
A[5][12]=0;
A[6][12]=0;
A[7][12]=0;
A[8][12]=0;
A[9][12]=0;
A[10][12]=0;
A[11][12]=-0.110546644894689;
A[12][12]=0.27582700508029;
A[13][12]=-0.107440961096779;
A[14][12]=0;
A[15][12]=0;
A[16][12]=0;
A[17][12]=0;
A[18][12]=0;
A[0][13]=0;
A[1][13]=0;
A[2][13]=0;
A[3][13]=0;
A[4][13]=0;
A[5][13]=0;
A[6][13]=0;
A[7][13]=0;
A[8][13]=0;
A[9][13]=0;
A[10][13]=0;
A[11][13]=0.0118434583833746;
A[12][13]=-0.110546644894689;
A[13][13]=0.218222694534443;
A[14][13]=-0.107440961096779;
A[15][13]=0.0110670374338972;
A[16][13]=0;
A[17][13]=0;
A[18][13]=0;
A[0][14]=0;
A[1][14]=0;
A[2][14]=0;
A[3][14]=0;
A[4][14]=0;
A[5][14]=0;
A[6][14]=0;
A[7][14]=0;
A[8][14]=0;
A[9][14]=0;
A[10][14]=0;
A[11][14]=0;
A[12][14]=0;
A[13][14]=-0.110546644894689;
A[14][14]=0.27582700508029;
A[15][14]=-0.107440961096779;
A[16][14]=0;
A[17][14]=0;
A[18][14]=0;
A[0][15]=0;
A[1][15]=0;
A[2][15]=0;
A[3][15]=0;
A[4][15]=0;
A[5][15]=0;
A[6][15]=0;
A[7][15]=0;
A[8][15]=0;
A[9][15]=0;
A[10][15]=0;
A[11][15]=0;
A[12][15]=0;
A[13][15]=0.0118434583833746;
A[14][15]=-0.110546644894689;
A[15][15]=0.218222694534443;
A[16][15]=-0.107440961096779;
A[17][15]=0.0110670374338972;
A[18][15]=0;
A[0][16]=0;
A[1][16]=0;
A[2][16]=0;
A[3][16]=0;
A[4][16]=0;
A[5][16]=0;
A[6][16]=0;
A[7][16]=0;
A[8][16]=0;
A[9][16]=0;
A[10][16]=0;
A[11][16]=0;
A[12][16]=0;
A[13][16]=0;
A[14][16]=0;
A[15][16]=-0.110546644894689;
A[16][16]=0.27582700508029;
A[17][16]=-0.107440961096779;
A[18][16]=0;
A[0][17]=0;
A[1][17]=0;
A[2][17]=0;
A[3][17]=0;
A[4][17]=0;
A[5][17]=0;
A[6][17]=0;
A[7][17]=0;
A[8][17]=0;
A[9][17]=0;
A[10][17]=0;
A[11][17]=0;
A[12][17]=0;
A[13][17]=0;
A[14][17]=0;
A[15][17]=0.0118434583833746;
A[16][17]=-0.110546644894689;
A[17][17]=0.218222694534443;
A[18][17]=-0.107440961096779;
A[0][18]=0;
A[1][18]=0;
A[2][18]=0;
A[3][18]=0;
A[4][18]=0;
A[5][18]=0;
A[6][18]=0;
A[7][18]=0;
A[8][18]=0;
A[9][18]=0;
A[10][18]=0;
A[11][18]=0;
A[12][18]=0;
A[13][18]=0;
A[14][18]=0;
A[15][18]=0;
A[16][18]=0;
A[17][18]=-0.110546644894689;
A[18][18]=0.27582700508029;

std::cout << std::endl << "A: " << std::endl;
print (A, locdim);

b[0]=0;
b[1]=0;
b[2]=0;
b[3]=0;
b[4]=0;
b[5]=0;
b[6]=0;
b[7]=0;
b[8]=0;
b[9]=0;
b[10]=0;
b[11]=0;
b[12]=0;
b[13]=0;
b[14]=0;
b[15]=-6.61491687424037e-05;
b[16]=0.000198392687107573;
b[17]=0.000792280471656408;
b[18]=0.00258639026622243;

//Compact A
double first_row [5] = {0, 0, A[0][0], A[0][1], A[0][2]};
Ac[0]= first_row;
double second_row [5] = {0, A[1][0], A[1][1], A[1][2], A[1][3]};
Ac[1]= second_row;
for (int nj = 2; nj < locdim - 2; ++nj)
for (int nk = 0; nk < 5; ++nk)
Ac[nj][nk] = A[nj][nk];
double second_to_last_row [5] = {A[17][15], A[17][16], A[17][17], A[17][18], 0};
Ac[17]= second_to_last_row;
double last_row [5] = {A[18][16], A[18][17], A[18][18], 0, 0};
Ac[18]= last_row;

// initialise to 0s lower triangular al
for (int nj = 0; nj < locdim; ++nj)
for (int nk = 0; nk < 2; ++nk)
al[nj][nk] = 0;

double d;

std::cout << std::endl << "b: " << std::endl;
for (int nj = 0; nj < locdim; ++nj)
std::cout << b[nj] << std::endl;

bandec(Ac, 2, 2, al, indx, d, locdim);
banbks(Ac, 2, 2, al, indx, b, locdim);

//std::cout << std::endl << "Ac: " << std::endl;
//printC(Ac, locdim);
//std::cout << std::endl << "al: " << std::endl;
//printC(al, locdim);

for(int i = 0; i < locdim; ++i){
for(int k = 0; k < locdim; ++k){
mult[i] += A[i][k] * b[k];
}
}

std::cout << std::endl;

for (int ni = 0; ni < locdim; ++ni) {
std::cout << "mult["<<ni<<"] = " << mult[ni]<<std::endl;
}

std::cout << std::endl;

for (int ni = 0; ni < locdim; ++ni) {
std::cout << "solution["<<ni<<"] = " << b[ni]<<std::endl;
}

for (int ni = 0; ni < locdim; ++ni)
{
//delete [] A[ni] ;
//delete [] Ac[ni] ;
delete [] al[ni] ;
}

//delete [] A ;
delete [] Ac;
delete [] al;
delete [] b;
delete [] indx;
delete [] mult ;

return 0;
}

• The question is perhaps more suited for maths SE (if you don't understand the algorithm) or stack overflow (if you struggle to implement it)? Dec 29 '20 at 15:20
• Besides @Kevin ‘s comment - did you have a look at the numerical recipes book by Press et al.? They have a section on sparse linear systems, including tridiag and banded tridiag. Dec 29 '20 at 17:30
• @Kermittfrog I picked up the book as soon as I read the question, in fact. Great book! They say it's badly coded but corrections are easy to make IMO. Dec 29 '20 at 18:10
• The $A$ matrix you present in the cut&paste excel bit (first code section in your question) is different from the input in the C++ program from the second code bit! Look at the last lines they are filled with nonzero value eg A[18][9] = 0.666666666666667;, A[18][10] = 0.52704627669473; they are supposed to be 0 for the matrix to be banded. Dec 30 '20 at 16:09
• I have also omitted to say that $A_c$ is supposed to store the elements in the principal diagonal in its $m_1$-th column. Its dimension is $N \times (m_1 + 1 + m_2)$, in this case $m_1 = m_2 = 2$ so diagonal elements $A_{ii} \, , \quad i = 1, \dots N$ should go in the central column ($m_1 = 2$) of $A_c$. Additionally, lower triangular matrix $L$ should be set of dimensions $N \times m_1$ but does not need to be initialised from values of $A$ as it is an output of bandec, leave it as a $19 \times 2$, 0-filled matrix. Dec 30 '20 at 16:18

You are looking to solve a linear system with a band diagonal system matrix $$A$$ of dimensions $$N \times N$$ and having the property $$m_1 = 2, \, m_2 = 2$$, respectively the number of sub-diagonals under and over the main diagonal of $$A$$. $$A$$ is empty elsewhere:

$$a_{ij} = 0 \qquad \text{for } j > i + m_2 \quad \text{ or } \quad i > j + m_1$$

In order to solve a linear system of equations $$A \, x = b$$ you may want to proceed as follows:

1. Resort to a compact representation of $$A$$, that is, a new matrix $$A_c$$ that includes only the elements in the nonzero diagonals, and as such is not square any more but with dimensions $$N \times (m_1 + 1 + m_2)$$. For instance, if $$m_1 = 2$$, $$m_2 = 2$$ as in the OP question, then $$A$$ and $$A_c$$ can be [EDIT: C++ subroutines below require that $$A_c$$ store the principal diagonal of $$A$$ in column $$m_1 + 1$$ (column m1 if index of first element is 0 as in C++)]: \begin{align} A = \begin{bmatrix} 1 & 2 & 3 & 0 & 0 & 0\\ 4 & 5 & 6 & 7 & 0 & 0\\ 8 & 9 & 1 & 2 & 3 & 0\\ 0 & 4 & 5 & 6 & 7 & 8\\ 0 & 0 & 9 & 1 & 2 & 3\\ 0 & 0 & 0 & 4 & 5 & 6 \end{bmatrix} \qquad & \text{then} \qquad A_c = \begin{bmatrix} 0 & 0 & 1 & 2 & 3\\ 0 & 4 & 5 & 6 & 7\\ 8 & 9 & 1 & 2 & 3\\ 4 & 5 & 6 & 7 & 8\\ 9 & 1 & 2 & 3 & 0\\ 4 & 5 & 6 & 0 & 0 \end{bmatrix} \end{align}

2. Perform a LU decomposition of $$A_c$$. You can use subroutine bandec from the book Numerical recipes in C++, that can be found here (not for free though). The subroutine overwrites the input matrix a($$A_c$$) with upper triangular matrix $$U$$ (EDIT: of dimensions $$N \times (m_1 + 1 + m_2)$$), while the lower triangular matrix $$L$$ (EDIT: of dimensions $$N \times m_1$$) is returned in al. Here below my adaptation of this subroutine:

#include <cmath>
#include <vector>
//#include "nr.h"   // this is a header provided with the book
#include "my_matrix.h"    // substitute nr.h with your custom matrix class

using namespace std;

void bandec(matrix& a, const int m1, const int m2, matrix& al, vector<int>& indx, double& d)
// a is input matrix A_c and upper triangular matrix on output,
// m1 m2 dimensions of subdiagonals,
// al lower triangular matrix on output of dimensions (a.nrows(), m1),
// indx records the row permutations,
// d = 1 or -1 depending on whether row interchanges are even or odd
{
const double TINY = 1.e-20;

int n = a.nrows();    // custom function from your "my_matrix.h" class
int mm = m1 + 1 + m2;
int l = m1;
int i = 0;
for (i = 0; i < m1; i++) {
for (int j = m1 - i; j < mm; j++)
a[i][j-l] = a[i][j];
l--;
for (int j = mm - l - 1; j < mm; j++)
a[i][j] = 0.0;
}
d = 1.0;
l = m1;
double dum = 0;
int k = 0;
// finds pivot for each row in a
for (k = 0; k < n; k++) {
dum = a[k][0];
i = k;
if (l < n)
l++;
for (int j = k + 1; j < l; j++) {
if (fabs(a[j][0] > fabs(dum)) {
dum = a[j][0];
i = j;
}
}

indx[k] = i + 1;
// in case input matrix a is singular
if (dum == 0.0)
a[k][0] = TINY;
if (i != k) {        // swap values in rows i and k
d = -d;
for (int j = 0; j < mm; j++) {
double aux = a[k][j];
a[k][j] = a[i][j];
a[i][j] = aux;
}
}
for (i = k + 1; i < l; i++) {   // perform elimination
dum = a[i][0] / a[k][0];
al[k][i-k-1] = dum;
for (int j = 1; j < mm,; j++)
a[i][j-1] = a[i][j] - dum * a[k][j];
a[i][mm-1] = 0.0;
}
}
}

1. Use output upper a and lower al triangular matrices (EDIT: and output vector indx tracking pivoted rows in al as well) from the LU decomposition of $$A_c$$ performed above in bandec, as inputs for the subroutine solving linear systems, banbks, again from the Numerical recipes book:
#include <vector>
#include "my_matrix.h"

using namespace std;

void banbks(matrix& a, const int m1, const int m2, matrix& al, vector<int>& indx, vector<double>& b)
// a, al , and indx inputs from bandec;
// m1 and m2 are the number of subdiagonals in starting system matrix A;
// b is the rhs of linear system Ax = b on input, solution vector x on output
{
int mm = m1 + 1 + m2;
int l = m1;
for (int k = 0; k < n; k++) {
int j = indx[k] - 1;
if (j != k) {
double aux = b[k];
b[k] = b[j];
b[j] = aux;
}
if (l < n)
l++;
for (j = k + 1; j < l; j++)
b[j] -= al[k][j-k-1] * b[k];
}
l = 1;
for (int i = n - 1; i >= 0; i--) {
double dum = b[i];
for (int k = 1; k < l; k++)
dum -= a[i][k] * b[k+i];
b[i] = dum / a[i][0];
if (l < mm)
l++;
}
}


Final result $$x$$ is yielded by banbks as overwrite of vector b.

EDIT: Recently added OP's C++ coded MWE, where he employs a $$19 \times 19$$ matrix, should be modified as follows IMO:

1. initialise al to dimension $$19 \times 2$$:
for (int ni = 0; ni < locdim; ++ni) {
A[ni] = new double[locdim];
Ac[ni] = new double[5];
al[ni] = new double[2];
}

1. populate Ac with principal diagonal of A as central (m1-th) column, and initialise also al with 0s. Delete loops that were in place previously to initialise the matrices like Ac[nj][4-nk] = A[nj][locdim - nk - 1]; etc EDITED AGAIN CODE INSIDE STARS:
    //Compact A
double first_row [5] = {0, 0, A[0][0], A[0][1], A[0][2]};
Ac[0]= first_row;
double second_row [5] = {0, A[1][0], A[1][1], A[1][2], A[1][3]};
Ac[1]= second_row;
// ******************** EDITED AGAIN ****************************
for (int nj = 2; nj < locdim - 2; ++nj){
int ii = 0;
for (int nk = nj-2; nk <= nj+2; ++nk)
Ac[nj][ii++] = A[nj][nk];
}
// ******************** END EDIT ********************************
double second_to_last_row [5] = {A[17][15], A[17][16], A[17][17], A[17][18], 0};
Ac[17]= second_to_last_row;
double last_row [5] = {A[18][16], A[18][17], A[18][18], 0, 0};
Ac[18]= last_row;

// initialise to 0s lower triangular al
for (int nj = 0; nj < locdim; ++nj)
for (int nk = 0; nk < 2; ++nk)
al[nj][nk] = 0;

1. In place of A, feed Ac to the subroutines instead:
bandec(Ac, 2, 2, al, indx, d, locdim);
banbks(Ac, 2, 2, al, indx, b, locdim);


Modifying the code this way, plus feeding it the correct banded matrix input A, should yield the correct solution.

EDITED AGAIN: modified code in point 2 just above here.

mult=b now:

mult[0] = -2.88544e-024
mult[1] = 1.23754e-023
mult[2] = -2.02013e-023
mult[3] = -4.36564e-023
mult[4] = -1.03398e-025
mult[5] = 1.77391e-022
mult[6] = -4.65703e-022
mult[7] = 6.67948e-023
mult[8] = 2.9034e-022
mult[9] = -1.37002e-021
mult[10] = 1.23415e-021
mult[11] = 2.90175e-021
mult[12] = -5.00941e-021
mult[13] = -3.86459e-021
mult[14] = -5.16425e-020
mult[15] = -6.61492e-005
mult[16] = 0.000198393
mult[17] = 0.00079228
mult[18] = 0.00258639


and the solution vector $$x$$ is:

solution[0] = 1.11822e-007
solution[1] = 2.45856e-007
solution[2] = 5.65224e-007
solution[3] = 1.19213e-006
solution[4] = 2.34089e-006
solution[5] = 4.84546e-006
solution[6] = 9.50624e-006
solution[7] = 1.93638e-005
solution[8] = 3.72459e-005
solution[9] = 7.41132e-005
solution[10] = 0.000142386
solution[11] = 0.000283238
solution[12] = 0.000544113
solution[13] = 0.00108235
solution[14] = 0.00207922
solution[15] = 0.00413597
solution[16] = 0.00849962
solution[17] = 0.0153931
solution[18] = 0.0153728

• Nice explanation. I understood until the last step where I am supposed to solve 2 systems, the first with a as the coefficient matrix (what is the right hand side vector?), the second with al as the matrix (again what is on right hand side?). What do I do with these two solutions, how do I get the 'final answer'. Sorry if it is obvious. Dec 29 '20 at 18:30
• @noob2 yes you basically want to decompose $A x = b$ in 2 separate systems to solve, $A = L \cdot U \rightarrow (LU) \cdot x = b \rightarrow L (U \cdot x) = b$ so you have $U x = y$, that gets solved by the second loop in above banbks, and $L y = b$, that is instead solved by the first loop. Given that $U$ and $L$ are triangular matrices, the advantage of the method comes from the fact that $L y = b$ can be trivially solved by forward substitution (1st loop), while $U x = y$ by backward substitution (2nd loop). Dec 29 '20 at 18:46
• Incidentally, it would be nice to have C++ syntax highlighted in the code snippets, but there seems to be no tag available for computer languages specifically, only one generic programming tag. Dec 29 '20 at 18:50
• It would be nice indeed. I seem to remember something about it being disabled because of server side computational burden. Maybe we could try again. Dec 29 '20 at 19:16
• @BobJansen you have been advocating this for quite a while I see, quant.meta.stackexchange.com/questions/1508/… But it kind of fell flat. As I report in that post, we need to submit a request to the main Stack Exchange committee citing at least 10 posts where syntax highlight would have been convenient to have. This is one, 9 to go ... Dec 29 '20 at 21:19