Serial access to collections

Random access is rather expensive

see detailed discussion in the paper


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

Serial access would suffice. It's faster, it's safer, and shows off the algorithm better


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