Matrix multiplication algorithm problem (C ++)
Hi, I am helping my friend write a program that multiplies two upper triangular matrices without matrix expansion. When I say no matrix expansion, I mean not filling the bottom of the top triangular matrix with zeros (the goal is to preserve space). I am reading in a matrix from a file that only has upper triangular values for matrices with the lower part omitted.
What I'm trying to accomplish is to write a function that has given an array and a pair of indices, will return the element that the expanded matrix will have at that location (0 for diagonals below and values above diagonal). I think about it using the usual matrix multiplication algorithm that uses this function whenever it needs to access an element. I've been working on this for several hours and I can't think of a way to convert the indices of a double matrix (i, j) (i goes along the rows) to indices of a single array and vice versa (as a reminder I'm using a one-dimensional array to store an upper triangular matrix) ... Any help would be greatly appreciated!
source to share
// mat is an array containing the upper triangle data for a square matrix of size n
// returns element at (i,j), or 0 for the lower triangle
int getFromTriangle(const int* mat, int n, int i, int j)
{
if (i > j) {
return 0; // lower triangle
} else {
return mat[j + (i*n) - i*(i+1)/2];
}
}
The proposal if
takes care of the bottom triangle. The clause else
calculates the array index as follows:
j + (i*n) - i*(i+1)/2
This is just the usual function of the index of a rectangular matrix minus the offset, which is exactly the i
th triangular number , since in any row i
, the store has omitted triangle(i)
elements.
source to share
You can split the algorithm into smaller, easily understandable chunks.
// Get the array index given the rank, the row, and the column.
int getArrayIndex(int rank, int row, int col)
{
return (row*rank + col - col*(col+1)/2);
}
// Get the term of a matrix, given the rank, the row, and the column.
int getMatrixTerm(int a[], int rank, int row, int col)
{
if ( col < row )
{
return 0;
}
else
{
return a[getArrayIndex(rank, row, col)];
}
}
// Get the term for a row and column resulting from mulitiplication.
int getMultipliedTerm(int a[], int b[], int rank, int row, int col)
{
int term = 0;
int k = 0;
for ( ; k < rank; ++k )
{
term += getMatrixTerm(a, rank, row, k)*getMatrixTerm(b, rank, k, col);
}
return term;
}
// Set the term in c given the rank, the row, and the column.
void setMultipliedTerm(int a[], int b[], int c[], int rank, int row, int col)
{
if ( j >= i )
{
c[getArrayIndex(rank, i, j)] = getMultipliedTerm(a, b, rank, i, j);
}
}
// High level function to multiply two upper triangular matrices
// The lower part of the matrix is not stored at all.
void multiply(int a[], int b[], int c[], int rank)
{
int i = 0;
int j = 0;
for ( i = 0; i < rank; ++i )
{
for ( j = 0; j < rank; ++j )
{
setMultipliedTerm(a, b, c, rank, i, j);
}
}
}
source to share