1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You should read Stroustrup's book The C++ programming language. Read about the RAII principle (Resource Acquisition Is Initialization) for an example. This avoids memory leaks. You have some new but no delete. In modern code you seldom use new at all. You better use std::vector or similar stuff that does the memory management for you. $\endgroup$ – Tobias Mar 26 '16 at 22:09
  • $\begingroup$ Ok, thing is I did so much work for this already switching to std::vector may make my life harder. Unless that is a simpler way. $\endgroup$ – Wolfy Mar 26 '16 at 23:12
  • $\begingroup$ You allocate memory for each row of the matrix. But this is not necessary -- you can algorithmically determine how many nonzero entries each row has, and consequently you could just allocate one single array for all entries. $\endgroup$ – Wolfgang Bangerth Mar 28 '16 at 2:07
1
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ CRS is a bit of an overkill for banded matrices. You know what the indices are, so there is no need to spend two to three times the memory to explicitly store the pattern. $\endgroup$ – Christian Waluga Mar 28 '16 at 12:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.