x3 - compactly storing banded matrices

• considering a banded matrix, $A:$

  double amat[4][4] = {
    {-2,1,0,0},
    {2,-2,2,0},
    {0,3,-2,3},
    {0,0,4,-2}
  };

• it is a tridiagonal matrix, ie. it has a 'band' above and below the main diagonal

• now, defining a structure with helper function to hold the band matrix information:

struct band_mat{
	long ncol;                  // matrix size (n × n)
	long nbrows;            // number of stored diagonals
	long nbands_up;    // diagonals above main
	long nbands_low;  // diagonals below main
	double *array;        // compact band storage
	
	// LAPACK-specific storage for solving
	long nbrows_inv;
	double *array_inv;
	int *ipiv;
};

tyedef struct band_mat band_mat;

int init_band_mat(band_mat *bmat, long nbands_lower, long nbands_upper, long n_columns) {
	bmat->nbrows = nbands_lower + nbands_upper + 1;
	bmat->ncol   = n_columns;
	bmat->nbands_up = nbands_upper;
	bmat->nbands_low= nbands_lower;
	bmat->array      = (double *) malloc(sizeof(double)*bmat->nbrows*bmat->ncol);
	bmat->nbrows_inv = bmat->nbands_up*2 + bmat->nbands_low + 1;
	bmat->array_inv  = (double *) malloc(sizeof(double)*(bmat->nbrows+bmat->nbands_low)*bmat->ncol);
	bmat->ipiv       = (int *) malloc(sizeof(int)*bmat->ncol);
	if (bmat->array==NULL||bmat->array_inv==NULL) {
		return 0;
  }  
  /* Initialise array to zero */
  long i;
  for (i=0;i<bmat->nbrows*bmat->ncol;i++) {
	bmat->array[i] = 0.0;
  }
  return 1;
};

void finalise_band_mat(band_mat *bmat) {
	free(bmat->array);
	free(bmat->array_inv);
	free(bmat->ipiv);
}

• now the banded matrix can be stored in the structure:

band_mat bmat;
long ncols = 4;
long nbands_low = 1;
long nbands_up = 1;
init_band_mat(&bmat, nbands_low, nbands_up, ncols);

• the matrix is compactly stored as:

 0   1   2   3  <- upper band
-2 -2 -2 -2  <- main diagonal
 2   3   4   0 <- lower band

• helpers needed to set/get elements using the original indices

/* get pointer to a location in the band matrix */
double *getp(band_mat *bmat, long row, long column) {
	int bandno = bmat->nbands_up + row - column;
	if(row<0 || column<0 || row>=bmat->ncol || column>=bmat->ncol ) {
		printf("Indexes out of bounds in getp: %ld %ld %ld \n",row,column,bmat->ncol);
		exit(1);
	}
	return &bmat->array[bmat->nbrows*column + bandno];
}

double setv(band_mat *bmat, long row, long column, double val) {
	*getp(bmat,row,column) = val;
	return val;
}

double getv(band_mat *bmat, long row, long column) {
	return *getp(bmat,row,column);
}

long i, j;
for(i=0; i<ncols; i++) {
	if(i > 0){setv(&bmat, i, i-1, 1.0 * (i + 1));};
	setv(&bmat, i, i, -2.0);
	if(i<ncols-1) {setv(&bmat,i,i+1,1.0*(i+1));}
	/* Uniform source term in heat equation */
	b[i] = 1.0;   
}

• now solve the matrix

int solve_Ax_eq_b(band_mat *bmat, double *x, double *b) {
	/* copy bmat array into the temporary store */
	int i, bandno;
	for(i=0; i<bmat->ncol; i++) { 
		for (bandno=0; bandno<bmat->nbrows; bandno++) {
			  bmat->array_inv[bmat->nbrows_inv*i+(bandno+bmat->nbands_low)] = bmat->array[bmat->nbrows*i+bandno];
		}
		x[i] = b[i];
	}
	
	long nrhs = 1;
	long ldab = bmat->nbands_low*2 + bmat->nbands_up + 1;
	int info = LAPACKE_dgbsv( LAPACK_COL_MAJOR, bmat->ncol, bmat->nbands_low, bmat->nbands_up, nrhs, bmat->array_inv, ldab, bmat->ipiv, x, bmat->ncol);
	return info;
}