Main Question:- If the elements in the band formed by the three diagonals of a tridiagonal matrix X are represented column-wise in an array Y, with X[1,1] being stored in Y1, then write an algorithm to determine the value of X[i, j], 1 <= i, j <= n from array Y.
I did some calculation in rough and found the size of 1D array required and index to row and column index for main diagonal and upper diagonal but failing to do so for lower diagonal First image Identified the diagonals and non zero elements and calculated the size of array to store in column major order Second image Relation between matrix indices (i and j) and array index(k)
CodePudding user response：
Summary
You'll find that (indexing from 1), for an NxN tridiagonal matrix X linearized into an array Y, the entry in row r:
- on the main diagonal of X is at index
(3 * r - 3)of Y for r = 1..N; - on the upper-diagonal of X is at
(3 * r - 2)of Y (for r = 1..N-1); - on the lower-diagonal of X is at
(3 * r - 1)of Y (for r = 2..N).
For an element with row r, column c:
- if |r - c| > 1 then the value is zero;
- if r == c then it is on the leading diagonal;
- if r = c - 1 then it is on the upper diagonal;
- if r = c 1 then it is on the lower diagonal.
Of course, in C, arrays are indexed from zero, not one.
Cogitations
NxN tri-diagonal matrix X
Linearized into array Y with 3*N-2 entries
Indexing from 0
X c = 0 1 2 3 4 5
----------------------
r = 0 | 1 2 0 0 0 0
r = 1 | 3 4 5 0 0 0
r = 2 | 0 6 7 8 0 0
r = 3 | 0 0 9 10 11 0
r = 4 | 0 0 0 12 13 14
r = 5 | 0 0 0 0 15 16
Linearized (D = diagonal, L = lower band, U = upper band):
D U L D U L D U L D U L D U L D
Y = { 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 }
- Elements on main diagonal at indexes: 0, 3, 6, 9, 12, 15
- Elements on upper diagonal at indexes: 1, 4, 7, 10, 13
- Elements on lower diagonal at indexes: 2, 5, 8, 11, 14
Non-zero elements have c = r { -1, 0, 1 }
X[r,c] = 0 if |r - c| > 1X[r,c] = 3 * r 0 if r == c 0Main diagonalX[r,c] = 3 * r 1 if r == c - 1Upper diagonalX[r,c] = 3 * r 2 if r == c 1Lower diagonal
Always subject to 0 <= r < N; 0 <= c < N.
Indexing from 1
X c = 1 2 3 4 5 6
----------------------
r = 1 | 1 2 0 0 0 0
r = 2 | 3 4 5 0 0 0
r = 3 | 0 6 7 8 0 0
r = 4 | 0 0 9 10 11 0
r = 5 | 0 0 0 12 13 14
r = 6 | 0 0 0 0 15 16
Linearized (D = diagonal, L = lower diagonal, U = upper diagonal):
D U L D U L D U L D U L D U L D
Y = { 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 }
- Elements on main diagonal at indexes: 1, 4, 7, 10, 13, 16
- Elements on upper diagonal at indexes: 2, 5, 8, 11, 14
- Elements on lower diagonal at indexes: 3, 6, 9, 12, 15
Non-zero elements have c = r { -1, 0, 1 }
X[r,c] = 0 if |r - c| > 1X[r,c] = 3 * r - 2 if r == c 0Main diagonalX[r,c] = 3 * r - 1 if r == c - 1Upper diagonalX[r,c] = 3 * r - 0 if r == c 1Lower diagonal
Always subject to 1 <= r <= N; 1 <= c <= N.
Check:
- X[3,3] = Y[3*3-2] = Y[7] = 7 ✓
- X[3,4] = Y[3*3-1] = Y[8] = 8 ✓
- X[3,2] = Y[3*3-0] = y[6] = 6 ✓
Column-major instead of row-major order
You are correct. Here you are considering row-major order. I wanted to store the elements in the array as:
{ 1, 3, 2, 4, 6, 5, 7, 9, 8, 10, 12, 11, 13, 15, 14, 16 }as per your example
So, you can make a parallel argument to the one I gave. Annotate your column-major array Z similarly to the way I annotated the row-major array Y:
D L U D L U D L U D L U D L U D
{ 1, 3, 2, 4, 6, 5, 7, 9, 8, 10, 12, 11, 13, 15, 14, 16 }
As before, non-zero elements have c = r { -1, 0, 1 } or, equivalently, r = c { -1, 0, 1 }.
Assuming 1-based indexing (note the change to use c for column over r for row), arrange an appropriate permutation of the additions:
X[r,c] = 0 if |c - r| > 1X[r,c] = 3 * c - 2 if c == r 0Main diagonalX[r,c] = 3 * c - 3 if c == r 1Upper diagonalX[r,c] = 3 * c - 1 if c == r - 1Lower diagonal
Always subject to 1 <= r <= N; 1 <= c <= N.
Check:
- X[3,3] = Z[3*3-2] = Z[7] = 7 ✓
- X[3,4] = Z[3*4-3] = Z[9] = 8 ✓
- X[3,2] = Z[3*2-1] = Z[5] = 6 ✓
CodePudding user response：
Assuming d is the main diagonal, a is the "above" diagonal, and b is the "below" diagonal:
--- --- --- --- --- --- --- ---
| d | a | 0 | 0 | 0 | 0 | 0 | 0 |
--- --- --- --- --- --- --- ---
| b | d | a | 0 | 0 | 0 | 0 | 0 |
--- --- --- --- --- --- --- ---
| 0 | b | d | a | 0 | 0 | 0 | 0 |
--- --- --- --- --- --- --- ---
| 0 | 0 | b | d | a | 0 | 0 | 0 |
--- --- --- --- --- --- --- ---
| 0 | 0 | 0 | b | d | a | 0 | 0 |
--- --- --- --- --- --- --- ---
| 0 | 0 | 0 | 0 | b | d | a | 0 |
--- --- --- --- --- --- --- ---
| 0 | 0 | 0 | 0 | 0 | b | d | a |
--- --- --- --- --- --- --- ---
| 0 | 0 | 0 | 0 | 0 | 0 | b | d |
--- --- --- --- --- --- --- ---
Note that the problem assumes indexes are "origin 1" (e.g. 1 <= i <= n) but C arrays want "origin 0" (e.g. 0 <= i < n)
We can use an "inner" struct to represent the Y array (e.g. elem_t) and the basic functions are:
#include <stdio.h>
#include <stdlib.h>
// sparse matrix Y vector
typedef struct {
int a; // above diagonal
int d; // main diagonal
int b; // below diagonal
} elem_t;
// sparse matrix
typedef struct {
int n;
elem_t *base;
} sparse_t;
// sparse_create -- create sparse tridiagonal matrix
sparse_t *
sparse_create(int n)
{
sparse_t *mtx = malloc(sizeof(*mtx));
mtx->n = n;
mtx->base = calloc(n,sizeof(*mtx->base));
return mtx;
}
// sparse_pointer -- point to sparse matrix element
int *
sparse_pointer(sparse_t *mtx,int coli,int rowi)
{
int n = mtx->n;
elem_t *eptr;
int *iptr = NULL;
// convert to origin 0
coli -= 1;
rowi -= 1;
do {
// check column bounds
if (coli < 0)
break;
if (coli >= n)
break;
// check row bounds
if (rowi < 0)
break;
if (rowi >= n)
break;
eptr = &mtx->base[coli];
// main diagonal
if (rowi == coli) {
iptr = &eptr->d;
break;
}
// above diagnonal
if ((rowi == (coli - 1)) && (coli > 0)) {
iptr = &eptr->a;
break;
}
// below diagnonal
if ((rowi == (coli 1)) && (coli < (n - 1))) {
iptr = &eptr->b;
break;
}
} while (0);
return iptr;
}
// sparse_getval -- get matrix value
int
sparse_getval(sparse_t *mtx,int coli,int rowi)
{
int *iptr = sparse_pointer(mtx,coli,rowi);
int val;
if (iptr != NULL)
val = *iptr;
else
val = 0;
return val;
}
// sparse_setval -- set matrix value
void
sparse_setval(sparse_t *mtx,int coli,int rowi,int val)
{
int *iptr = sparse_pointer(mtx,coli,rowi);
if (iptr != NULL)
*iptr = val;
}
By contrast, a regular/full matrix would be:
#include <stdio.h>
#include <stdlib.h>
// regular/full matrix
typedef struct {
int n;
int *base;
} regmtx_t;
// regmtx_create -- create regular matrix
regmtx_t *
regmtx_create(int n)
{
regmtx_t *mtx = malloc(sizeof(*mtx));
mtx->n = n;
mtx->base = calloc(n * n,sizeof(*mtx->base));
return mtx;
}
// regmtx_pointer -- point to regular matrix element
int *
regmtx_pointer(regmtx_t *mtx,int coli,int rowi)
{
int *iptr = NULL;
// convert to origin 0
coli -= 1;
rowi -= 1;
do {
// check column bounds
if (coli < 0)
break;
if (coli >= mtx->n)
break;
// check row bounds
if (rowi < 0)
break;
if (rowi >= mtx->n)
break;
iptr = &mtx->base[(coli * mtx->n) rowi];
} while (0);
return iptr;
}