Trouble with algorithm for multiplying matrix (c++)

3
votes

Hi I'm helping a friend of mine write a program that multiplies two upper triangular matrices without expanding the matrices. When I say without expanding the matrix, I mean without filling the lower part of the upper triangular matrix with zeros (the goal is to conserve space). I'm reading in the matrix from a file which only has the upper triangular values for the matrices with the lower part omitted.

What I'm trying to accomplish is to write a function which given an array and a pair of indices, would return an element that an expanded matrix would have at that location (0 for below-diagonal, and values above-diagonal). I'm thinking of then 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 come up with a way to transform double matrix indices (i,j) (i goes along rows) to single array indices and vice versa (Just as a reminder I'm using a single dimensional array to store the upper triangular matrix). Any help would be greatly appreciated!

3
c++arraysmatrix
For j == i, elements that would be on the main diagonal are at 0, n, n + (n-1), n + (n-1) + (n-2) and so on (figuring out the ith element of this sequence is left as an exercise for the reader). For j > i, just add j - i to the result of the previous calculation. – Igor Tandetnik

3 Answers

0
votes
// 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 if clause takes care of the lower triangle. The else clause computes the array index like this:

j + (i*n) - i*(i+1)/2

Which is just the regular rectangular-matrix index function minus an offset which is exactly the ith triangular number, because in any row i, the storage has omitted triangle(i) elements.

0
votes

You can divide your algorithm to 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);
      }
   }
}
0
votes

if p = &a[0][0] then a[i][j] is same as *(p+i*row_size+j) in c++. where p is pointer to the same type of data as the matrix elements. I hope this is what you want, and that you know pointers. all the best.