Código C++ Algoritmo Thomas para resolver una pentadiagonal Ax=b

Question

Estoy buscando resolver $Ax=b$ para $x$ donde $A$ es una matriz cuadrada pentadiagonal (los elementos de las diagonales superior e inferior pueden sin embargo ser iguales a cero) y $x$ , $b$ dos vectores del mismo tamaño.

La dimensión del problema puede ser grande. $A$ puede ser de 1000x1000 o más. No he podido encontrar una implementación en C++ de un algoritmo similar al de Thomas ([ Algoritmo Thomas ][1]). ¿Alguien se ha encontrado con uno?

EDITAR :

Estoy realizando esta prueba y la solución no es igual al producto $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

Se trata de resolver un sistema lineal con una matriz de sistema diagonal de banda $A$ de las dimensiones $N \times N$ y tener la propiedad $m_1 = 2, \, m_2 = 2$ respectivamente, el número de subdiagonales por debajo y por encima de la diagonal principal de $A$ . $A$ está vacío en otro lugar:

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

Para resolver un sistema lineal de ecuaciones $A \, x = b$ puede proceder de la siguiente manera:

1. Recurrir a un compacto representación de $A$ es decir, una nueva matriz $A_c$ que incluye sólo los elementos en las diagonales no nulas, y como tal ya no es cuadrado sino con dimensiones $N \times (m_1 + 1 + m_2)$ . Por ejemplo, si $m_1 = 2$ , $m_2 = 2$ como en la pregunta del OP, entonces $A$ y $A_c$ puede ser [EDIT: Las subrutinas C++ que aparecen a continuación requieren que $A_c$ almacenan la diagonal principal de $A$ en la columna $m_1 + 1$ (columna m1 si el índice del primer elemento es 0 como en 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. Realice una Descomposición LU de $A_c$ . Puede utilizar la subrutina bandec del libro Recetas numéricas en C++ que se puede encontrar aquí (aunque no de forma gratuita). La subrutina sobrescribe la matriz de entrada a ( $A_c$ ) con matriz triangular superior $U$ ( EDIT: de las dimensiones $N \times (m_1 + 1 + m_2)$ ), mientras que la matriz triangular inferior $L$ ( EDIT: de las dimensiones $N \times m_1$ ) se devuelve en al . A continuación mi adaptación de esta subrutina:

   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. Utilizar la salida superior a y más bajo al matrices triangulares ( EDIT: y vector de salida indx seguimiento de las filas pivotantes en al también ) de la descomposición LU de $A_c$ realizado arriba en bandec como entradas para la subrutina que resuelve los sistemas lineales, banbks , de nuevo desde el Recetas numéricas libro:

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

Resultado final $x$ se obtiene mediante banbks como sobreescritura del vector b .

---

EDITAR: Recientemente se ha añadido el MWE codificado en C++ de OP, donde emplea un $19 \times 19$ matriz, debería modificarse de la siguiente manera OMI:

1. inicializar al a la dimensión $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. rellenar Ac con la diagonal principal de A como central ( m1 -ésima) columna, e inicializar también al con 0 s. Elimina los bucles que había antes para inicializar las matrices como Ac[nj][4-nk] = A[nj][locdim - nk - 1]; etc. Editado de nuevo el código interior 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. En lugar de A , alimentación Ac a las subrutinas en su lugar:

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

Modificando el código de esta manera, además de alimentarlo con la entrada correcta de la matriz de bandas A debería dar la solución correcta.

---

EDITADO DE NUEVO: código modificado en el punto 2 justo por encima de aquí.

mult = b ahora:

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

y el solution vector $x$ es:

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