Creating a matrix that saves storage

Question

Suppose $A\in\mathbb{R}^{n\times n}$ is a banded matrix, i.e., a matrix with all of its nonzero elements on the main diagonal, i.e., $\alpha_{i,i}\neq 0$, the first superdiagonal, i.e., $\alpha_{i,i+1}\neq 0$, through the $k$-th superdiagonal, i.e., $\alpha_{i,i+k}\neq 0$, the first subdiagonal, i.e., $\alpha_{i-1,i}\neq 0$, through the $k$-th subdiagonal, i.e., $\alpha_{i,i+k}\neq 0$. All elements not on these diagonal are $0$. For $k = 4$ and $ n = 15$ the pattern is

\begin{pmatrix} * & * & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ * & * & * & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ * & * & * & * & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ * & * & * & * & * & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ * & * & * & * & * & * & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & * & * & * & * & * & * & * & * & * & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & * & * & * & * & * & * & * & * & * & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & * & * & * & * & * & * & * & * & * & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & * & * & * & * & * & * & * & * & * & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & * & * & * & * & * & * & * & * & * & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & * & * & * & * & * & * \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & * & * & * & * & * \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & * & * & * & * \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & * & * & * \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & * & * \\ \end{pmatrix}

Now we need to store this matrix as efficiently as possible. I wrote my code in C++ and the output is correct and the way I did it I believe saves the storage since I do not store any of the zeros.

Here is my code that I did:

double** banded_matrix(int n) {
    double **data = new double *[n];
    data[0] = new double[3];
    for(int i = 1; i <= 3; i++){
        data[i] = new double[4];
    }
    for(int i = 4; i <= n-3; i++){
        data[i] = new double[5];
    }
    for(int i = n-4; i <= n-2; i++){
        data[i] = new double[4];
    }
    data[n-1] = new double[3];
    return data;
}

double get(double** A, int n, int i, int j) {
    double result;
    // main tridiagonal
    if(abs(i-j) < 2){
        if(i == 0){
            result = A[i][j];
        } else if(i <= 3){
            result = A[i][j-(i-2)];
        } else if(i >= 3){
            result = A[i][j-(i-1)];
        }
    // 4th superdiagonal
    } else if (j - i == 4){
        if(i == 0 || i == n-1) {
            result = A[i][2];
        } else if((i >= 1 && i <= 3) || (i >= n - 4 && i <= n - 2)) {
            result = A[i][3];
        } else {
            result = A[i][4];
        }
    // 4th subdiagonal
    } else if (i - j == 4){
        result = A[i][0];
    // the rest of the matrix    
    } else {
        result = 0.;
    }
    return result;
}

void set(double** A, int n, int i, int j, double val) {
    // main tridiagonal
    if(abs(i-j) < 2){
        if(i == 0){
            A[i][j] = val;
        } else if(i <= 3){
            A[i][j-(i-2)] = val;
        } else if(i >= 3){
            A[i][j-(i-1)] = val;
        }
    // 4th superdiagonal
    } else if (j - i == 4){
        if(i == 0 || i == n-1) {
            A[i][2] = val;
        } else if((i >= 1 && i <= 3) || (i >= n - 4 && i <= n - 2)) {
            A[i][3] = val;
        } else {
            A[i][4] = val;
        }
    // 4th subdiagonalf
    } else if (i - j == 4){
        A[i][0] = val;
    // the rest of the matrix    
    } else {
        cout << "cannot set element (" << i <<"," << j <<") in matrix" << endl;
    }
}

int main()
{
    int N = 15;

    //Initialize A
    double **A = banded_matrix(N);
    create(A,N,0);
//Print A
    cout << "A = " << endl;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
           cout << get(A,N,i,j) << "   ";        
        }
       cout << endl;
    }
  cout << endl;


}

I just want to know if my method is efficient or if I should have done it a different way.

Answer 1

Reiterating some of points made by Tobias and Wolfgang Bangerth:

• I see no delete to free memory! I highly suggest running valgrind on your code to find any memory leaks.
• Storing your matrix with two dimensional arrays probably isnt the best idea especially for large systems. Instead I suggest looking into other sparse matrix storage formats, e.g. Compressed Row Storage (CRS), Compressed Column Storage, ...
• using std::vector is also a good idea
• I think writing your functions to return double pointer variables is needlessly messy.

An example of using the CRS format:

#include <iostream>

// Sample test matrix
// [2,1,0,0,1,0,0,0
//  1,2,1,0,0,1,0,0
//  0,1,2,1,0,0,1,0
//  0,0,1,2,1,0,0,1
//  1,0,0,1,2,1,0,0
//  0,1,0,0,1,2,1,0
//  0,0,1,0,0,1,2,1
//  0,0,0,1,0,0,1,2]


// In CRS this matrix can be represented by the three arrays:
int rows[9] = {0,3,7,11,15,19,23,27,30};
int columns[30] =  {0,1,4,0,1,2,5,1,2,3,6,2,3,4,7,0,3,4,5,1,4,5,6,2,5,6,7,3,6,7};
double values[30] =    {2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2};


int main()
{
    // matrix multiply y = A*x
    double x[8] = {1,2,3,4,5,6,7,8};
    double y[8];
    for(unsigned int i=0; i<8; i++){
        y[i] = 0;
        for(unsigned int j=rows[i]; j<rows[i+1]; j++){
            y[i] += values[j]*x[columns[j]];
        }
    }

    // print arrays y 
    for(unsigned int i=0; i<8; i++){
        std::cout<<y[i]<<std::endl;
    }
}

Using the CRS format is efficient because it stores no non-zeros. Another nice things about this format is that by storing the sparse matrix with three one dimensional arrays means your memory is contiguous.