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

Question

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 n = dim;    // from your "my_matrix.h" header
    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;
}

Answer 1

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;
        }
    }
}   

3. 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 n = a.nrows();    // from your "my_matrix.h" header
    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];
}

2. 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;

3. 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