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.
new
but nodelete
. In modern code you seldom usenew
at all. You better usestd::vector
or similar stuff that does the memory management for you. $\endgroup$ – Tobias Mar 26 '16 at 22:09