Serial access to collections

Random access is rather expensive

But often the random access to collection's elements isn't even necessary

Example: matrix multiplication C = A*B

Mathematical algorithm

Cij = SUM{ Aik * Bkj }

Serial algorithm

j-th column of C = A * j-th column of B =
{i-th column of A * j-th column of B}

for(Matrix_col_stream cj(C), bj(B);
    !cj.total_eof();
    cj.next_column(), bj.next_column())
  for(Matrix_row_stream ai(A);
      !ai.total_eof();
      ai.next_row(), bj.rewind())
  {
    register double cij = 0;
    while( !ai.eof() )
      cij += ai.get() * bj.get();
    cj.put(cij);
  }

Actual implementation

(from the LinAlg package)

void Matrix::_AmultB(const Matrix& A, const Matrix& B)
{
  A.is_valid();
  B.is_valid();

  if( A.ncols != B.nrows || A.col_lwb != B.row_lwb )
    A.info(), B.info(),
    _error("matrices above cannot be multiplied");

  allocate(A.nrows,B.ncols,A.row_lwb,B.col_lwb);

  register REAL * arp;  // Pointer to the i-th row of A
           REAL * bcp = B.elements; // Pointer to the j-th col of B
  register REAL * cp = elements;    // C is to be traversed in the natural
  while( cp < elements + nelems )   // order, col-after-col
  {
    for(arp = A.elements; arp < A.elements + A.nrows; )
    {
      register double cij = 0;
      register REAL * bccp = bcp;               // To scan the jth col of B
      while( arp < A.elements + A.nelems )      // Scan the i-th row of A and
         cij += *bccp++ * *arp, arp += A.nrows; // the j-th col of B
      *cp++ = cij;
      arp -= A.nelems - 1;  // arp points to (i+1)-th row
    }
    bcp += B.nrows;  // We're done with j-th col of both
  }                  // B and C. Set bcp to the (j+1)-th col

  assert( cp == elements + nelems && bcp == B.elements + B.nelems );
}

Mathematical algorithm seems to require a direct access to matrices. But it's not necessary. The whole deal is just traversing matrices in some specific ways.

The snippet with Matrix_stream is pseudocode. Although with Matrix_row_stream/Matrix_col_stream implementations (which is a piece of cake) it becomes a full-blown working code.

Note, there is no direct access to any element of the matrix anywhere in the code; everything is done in terms of the 'current' (element, row, col) and the next (element, row, col). If you're saying it's all so corny, here's a less trivial example

Next | Table of contents