matrixMult.hh

Go to the documentation of this file.

 
#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