operator/matrixMult.hh

Go to the documentation of this file.

/* matrix-matrix-multiplication function
 */

#ifndef matrixmult_hh
#define matrixmult_hh

#include "matrixIterator.hh"
#include "basics/exceptions.hh"
#include "basics/debug.hh"

#define MatrixMultRowSort_D 0

namespace concepts {

  // ********************************************** matrixMultiplyRowSorting **

  template<class F, class G, class H>
  void matrixMultiplyRowSorting(const F& factL, const G& factR,
                                Matrix<H>& dest) {
    conceptsAssert(factL.nofCols() == factR.nofRows(), Assertion());
    conceptsAssert(factL.nofRows() == dest.nofRows(), Assertion());
    conceptsAssert(factR.nofCols() == dest.nofCols(), Assertion());
    DEBUGL(MatrixMultRowSort_D, "factL = " << factL << ", factR = " << factR);

    // current column and row in the left matrix
    uint col, row;
    // current value in left matrix
    typename F::type value;

    typename F::const_iterator factL_end = factL.end();
    typename G::const_iterator factR_end = factR.end();
    // iterator over one row in factR
    for(typename F::const_iterator i = factL.begin(); i != factL_end; ++i) {
      DEBUGL(MatrixMultRowSort_D, "i = " << i);
      col = i.col();
      row = i.row();
      value = *i;
      for(typename G::const_iterator j = factR.begin(col);
          j != factR_end && j.row() == col; ++j) {
        DEBUGL(MatrixMultRowSort_D, "j = " << j);
        dest(row, j.col()) += value * *j;
      }
    }
    DEBUGL(MatrixMultRowSort_D, "done");
  }

} // namespace concepts

#endif // matrixmult_hh