/* this is generic */template <class Matrix, class VecX, class VecZ>inline voidmult_generic__(const Matrix& A, const VecX& xx, VecZ& zz) MTL_THROW_ASSERTION{  MTL_ASSERT(A.nrows() <= zz.size(), "mtl::mult()");  MTL_ASSERT(A.ncols() <= xx.size(), "mtl::mult()");  typedef typename matrix_traits<Matrix>::value_type T;  typename Matrix::const_iterator i;  typename Matrix::OneD::const_iterator j, jend;  typename VecX::const_iterator x = xx.begin();  typename VecZ::iterator z = zz.begin();  for (i = A.begin(); i != A.end(); ++i) {    j = (*i).begin(); jend = (*i).end();    for (; j != jend; ++j)      z[j.row()] += *j * x[j.column()];  }}template <class Matrix, class VecX, class VecZ>inline voidmult_shape__(const Matrix& A, const VecX& x, VecZ& z,             banded_tag) MTL_THROW_ASSERTION{  mult_generic__(A, x, z);}/* this is fast */template <class Matrix, class VecX, class VecZ>inline voidrect_mult(const Matrix& A, const VecX& xx, VecZ& zz,           row_tag, dense_tag) MTL_THROW_ASSERTION{  MTL_ASSERT(A.nrows() <= zz.size(), "mtl::mult()");  MTL_ASSERT(A.ncols() <= xx.size(), "mtl::mult()");  typedef typename matrix_traits<Matrix>::value_type T;  typename Matrix::const_iterator i, iend;  typename Matrix::OneD::const_iterator j, jend;  typename VecX::const_iterator x = xx.begin();  typename VecZ::iterator z = zz.begin();  i = A.begin();  iend = A.end();  for (; i != iend; ++i) {    j = (*i).begin(); jend = (*i).end();    T tmp = z[j.row()];    for (; j != jend; ++j)      tmp += *j * x[j.column()];    z[j.row()] = tmp;  }}/*  This is slow */template <class Matrix, class VecX, class VecZ>inline voidrect_mult(const Matrix& A, const VecX& xx, VecZ& zz,          column_tag, dense_tag) MTL_THROW_ASSERTION{  MTL_ASSERT(A.nrows() <= zz.size(), "mtl::mult()");  MTL_ASSERT(A.ncols() <= xx.size(), "mtl::mult()");  typedef typename matrix_traits<Matrix>::value_type T;  typedef typename matrix_traits<Matrix>::size_type Int;  typename VecX::const_iterator x = xx.begin();  typename VecZ::iterator z = zz.begin();  typename Matrix::const_iterator i, iend;  typename Matrix::OneD::const_iterator j, jend;  i = A.begin();  iend = A.end();  for (; i != iend; ++i) {    j = (*i).begin(); jend = (*i).end();    for (; j != jend; ++j)      z[j.row()] += *j * x[j.column()];  }}// x is sparse// A is column orientedtemplate <class Matrix, class VecX, class VecY>void rect_mult(const Matrix& A, const VecX& x, VecY& y, 	  column_tag, sparse_tag) {  typename VecX::const_iterator xi = x.begin();  for (; xi != x.end(); ++xi) {    mtl::add(mtl::scaled(A[xi.index()], *xi), y);  }}  // x is sparse// A is row orientedtemplate <class Matrix, class VecX, class VecY>void rect_mult(const Matrix& A, const VecX& x, VecY& y, 	  row_tag, sparse_tag){  typename Matrix::const_iterator Ai;  for (Ai = A.begin(); Ai != A.end(); ++Ai) {    y[Ai.index()] = mtl::dot(*Ai, x); // this is a sparse dot  }}template <class Matrix, class VecX, class VecZ>inline voidmult_shape__(const Matrix& A, const VecX& x, VecZ z,             rectangle_tag) MTL_THROW_ASSERTION{  typedef typename matrix_traits<Matrix>::orientation Orien;  typedef typename linalg_traits<VecX>::sparsity SparseX;  rect_mult(A, x, z, Orien(), SparseX());}template <class Matrix, class VecX, class VecZ>inline voidmult_shape__(const Matrix& A, const VecX& x, VecZ& z,              triangle_tag){  mult_shape__(A, x, z, rectangle_tag());  if (A.is_unit()) {    /* actually, this still isn't quite right,       should do       add_n(x, z, z, MTL_MIN(A.nrows(), A.ncols()));       instead       */    if (z.size() <= x.size())      mtl::add(z, x, z);    else      mtl::add(x, z, z);        }}template <class Matrix, class VecX, class VecZ>inline voidmult_symm__(const Matrix& A, const VecX& x, VecZ& z, row_tag){  typedef typename matrix_traits<Matrix>::value_type T;  typename Matrix::const_iterator i;  typename Matrix::OneD::const_iterator j, jend;  for (i = A.begin(); i != A.end(); ++i) {    T tmp = z[i.index()];    j = (*i).begin();    jend = (*i).end();    if (A.is_upper()) {      tmp += *j * x[j.column()];      ++j;    } else      --jend;    for (; j != jend; ++j) {      /* normal side */      tmp += *j * x[j.column()];      /* symmetric side */      z[j.column()] += *j * x[j.row()];    }    if (A.is_lower())      tmp += *j * x[j.column()];    z[i.index()] = tmp;  }}template <class Matrix, class VecX, class VecZ>inline voidmult_symm__(const Matrix& A, const VecX& x, VecZ& z, column_tag){  typedef typename matrix_traits<Matrix>::value_type T;  typename Matrix::const_iterator i;  typename Matrix::OneD::const_iterator j, jend;  for (i = A.begin(); i != A.end(); ++i) {    T tmp = T(0);    j = (*i).begin();    jend = (*i).end();    if (A.is_lower()) {      z[j.column()] += *j * x[j.column()];      ++j;    } else      --jend;    for (; j != jend; ++j) {      /* normal side */      z[j.row()] += *j * x[j.column()];      /* symmetric side */      tmp += *j * x[j.row()];    }    if (A.is_upper())      tmp += *j * x[j.row()];          z[i.index()] += tmp;  }}template <class Matrix, class VecX, class VecZ>inline voidmult_shape__(const Matrix& A, const VecX& x, VecZ& z,              symmetric_tag){  typedef typename matrix_traits<Matrix>::orientation Orien;  mult_symm__(A, x, z, Orien());}//: Multiplication:  <tt>z <- A x + y</tt>//!category: algorithms//!component: function //!definition: mtl.h//!precond:  <TT>A.nrows() <= y.size()</TT>//!precond:  <TT>A.nrows() <= z.size()</TT>//!precond:  <TT>A.ncols() <= x.size()</TT>//!precond:  no aliasing in the arguments//!example: symm_sparse_vec_prod.cc//!typereqs: <tt>Matrix::value_type</tt>, <tt>VecX::value_type</tt>, <tt>VecY::value_type</tt>, and <tt>VecZ::value_type</tt> must be the same type//!typereqs: the multiplication operator must be defined for <tt>Matrix::value_type</tt>//!typereqs: the addition operator must be defined for <tt>Matrix::value_type</tt>template <class Matrix, class VecX, class VecY, class VecZ>inline voidmult(const Matrix& A, const VecX& x, const VecY& y, MTL_OUT(VecZ) z_)  MTL_THROW_ASSERTION{  VecZ& z = const_cast<VecZ&>(z_);  mtl::copy(y, z);  typedef typename matrix_traits<Matrix>::shape Shape;  mult_shape__(A, x, z, Shape());}//: Matrix Vector Multiplication:  <tt>y <- A x</tt>//// Multiplies matrix A times vector x and stores the result in vector y.// <p>// Note: ignore the <tt>oned_tag</tt> parameter and the underscores in// the name of this function.////!category: algorithms//!component: function//!definition: mtl.h//!example: general_matvec_mult.cc, banded_matvec_mult.cc, symm_matvec_mult.cc//!precond:  <TT>A.nrows() <= y.size()</TT>//!precond:  <TT>A.ncols() <= x.size()</TT>//!precond:  x and y not same vector//!example: symm_matvec_mult.cc//!typereqs: <tt>Matrix::value_type</tt>, <tt>VecX::value_type</tt>, and <tt>VecY::value_type</tt> must be the same type//!typereqs: the multiplication operator must be defined for <tt>Matrix::value_type</tt>//!typereqs: the addition operator must be defined for <tt>Matrix::value_type</tt>template <class Matrix, class VecX, class VecY>inline voidmult_dim__(const Matrix& A, const VecX& x, VecY& y, oned_tag) MTL_THROW_ASSERTION{  mtl::mult(A, x, mtl::scaled(y, 0), y);#if 0  typedef typename matrix_traits<Matrix>::shape Shape;  mult_shape__(A, x, y, Shape());#endif}template <class Matrix, class VecX, class VecY>inline voidmult_add(const Matrix& A, const VecX& x, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION{  VecY& y = const_cast<VecY&>(y_);  typedef typename matrix_traits<Matrix>::shape Shape;  mult_shape__(A, x, y, Shape());}//: simple 3 loop version of matmat mult//!noindex:template <class MatA, class MatB, class MatC, class Orien>inline voidsimple_mult(const MatA& A, const MatB& B, MatC& C, dense_tag, Orien){  typedef typename matrix_traits<MatA>::size_type Int;  typename MatA::const_iterator A_k;  typename MatA::OneD::const_iterator A_ki;  A_k = A.begin();  while (not_at(A_k, A.end())) {    for (Int j = 0; j < B.ncols(); ++j) {      A_ki = (*A_k).begin();      while (not_at(A_ki, (*A_k).end())) {        Int k = A_ki.column();        Int i = A_ki.row();        C(i,j) += *A_ki * B(k,j);        ++A_ki;      }    }    ++A_k;  }}/* Assumes A and B are also row oriented */template <class MatrixA, class MatrixB, class MatrixC>inline voidsimple_mult(const MatrixA& A, const MatrixB& B, MatrixC& C,             sparse_tag, row_tag){  typedef typename matrix_traits<MatrixA>::value_type T;  typedef typename matrix_traits<MatrixA>::size_type Int;  T scal;  Int len = 0;  Int jj, k;  Int nzmax = C.capacity();  Int M = A.nrows();  Int N = B.ncols();  dense1D<Int> ic(M + 1, 0);  dense1D<Int> jc(nzmax);  dense1D<T> c(nzmax);    typedef typename dense1D<Int>::iterator di_iter;  typedef typename dense1D<T>::iterator dt_iter;    compressed1D<T> tmp1(N), tmp2(N), tmp3(N);  tmp1.reserve(N);  tmp2.reserve(N);    typedef typename compressed1D<T>::iterator tmpiter;    typename MatrixA::const_iterator Ai;  typename MatrixA::Row::const_iterator Aij;    for (Ai = A.begin(); Ai != A.end(); ++Ai) {    copy(C[Ai.index()], tmp1);    for (Aij = (*Ai).begin(); Aij != (*Ai).end(); ++Aij) {      scal = *Aij;      jj = Aij.column();      // add B[jj] and tmp1 into tmp2      add(mtl::scaled(B[jj], scal), tmp1, tmp2);      tmp1.clear();      // swap tmp1 and tmp2      tmp3 = tmp1; tmp1 = tmp2; tmp2 = tmp3;    }    // copy tmp1 into C[ii]    k = len;    if (k + tmp1.nnz() > nzmax) {      std::cerr << "Not enough work space, increase capacity of C" << std::endl;      return;    }    for (tmpiter t = tmp1.begin(); t != tmp1.end(); ++t, ++k) {      c[k] = *t;      jc[k] = t.index();    }        len += tmp1.nnz();    ic[Ai.index() + 1] = len;  }  typedef typename matrix<T, rectangle<>,     compressed<Int, external>,     row_major>::type  SpMat;  SpMat CC(M, N, len, c.data(), ic.data(), jc.data());  copy(CC, C);}/* Assumes A and B are also column oriented */template <class MatrixA, class MatrixB, class MatrixC>inline voidsimple_mult(const MatrixA& A, const MatrixB& B, MatrixC& C,             sparse_tag, column_tag){  typedef typename matrix_traits<MatrixA>::value_type T;  typedef typename matrix_traits<MatrixA>::size_type Int;  T scal;  Int len = 0;  Int kk, k;  Int nzmax = C.capacity();  Int M = A.nrows();  Int N = B.ncols();  dense1D<Int> ic(N + 1, 0);  dense1D<Int> jc(nzmax);  dense1D<T> c(nzmax);    typedef typename dense1D<Int>::iterator di_iter;  typedef typename dense1D<T>::iterator dt_iter;    compressed1D<T> tmp1(M), tmp2(M), tmp3(M);  tmp1.reserve(M);  tmp2.reserve(M);    typedef typename compressed1D<T>::iterator tmpiter;    typename MatrixB::const_iterator Bj;  typename MatrixB::Column::const_iterator Bjk;    for (Bj = B.begin(); Bj != B.end(); ++Bj) {    copy(C[Bj.index()], tmp1);    for (Bjk = (*Bj).begin(); Bjk != (*Bj).end(); ++Bjk) {      scal = *Bjk;      kk = Bjk.row();      // add A[kk] and tmp1 into tmp2      add(mtl::scaled(A[kk], scal), tmp1, tmp2);      tmp1.clear();      // swap tmp1 and tmp2      tmp3 = tmp1; tmp1 = tmp2; tmp2 = tmp3;    }    // copy tmp1 into C[ii]    k = len;    if (k + tmp1.nnz() > nzmax) {      std::cerr << "Not enough work space, increase capacity of C" << std::endl;      return;    }    for (tmpiter t = tmp1.begin(); t != tmp1.end(); ++t, ++k) {      c[k] = *t;      jc[k] = t.index();    }        len += tmp1.nnz();    ic[Bj.index() + 1] = len;  }    typedef typename matrix<T, rectangle<>,                  compressed<Int, external>,                  column_major>::type  SpMat;  SpMat CC(M, N, len, c.data(), ic.data(), jc.data());  copy(CC, C);}//: Symmetric version, row-major//!noindex:template <class MatA, class MatB, class MatC>inline voidsymm_simple_mult(const MatA& A, const MatB& B, MatC& C, row_tag){  typedef typename matrix_traits<MatA>::size_type Int;  typename MatA::const_iterator A_k;  typename MatA::OneD::const_iterator A_ki, A_kiend;  A_k = A.begin();  while (not_at(A_k, A.end())) {    for (Int j = 0; j < B.ncols(); ++j) {      A_ki = (*A_k).begin();      A_kiend = (*A_k).end();      Int k = A_ki.column();      Int i = A_ki.row();      if (A.is_upper()) { /* handle the diagonal elements */        C(i,j) += *A_ki * B(k,j);        ++A_ki;      } else        --A_kiend;      while (not_at(A_ki, A_kiend)) {        k = A_ki.column();        i = A_ki.row();        C(i,j) += *A_ki * B(k,j);        C(k,j) += *A_ki * B(i,j);        ++A_ki;      }      k = A_ki.column();      i = A_ki.row();      if (A.is_lower())        C(i,j) += *A_ki * B(k,j);    }    ++A_k;  }}//: Symmetric version, column-major//!noindex:template <class MatA, class MatB, class MatC>inline voidsymm_simple_mult(const MatA& A, const MatB& B, MatC& C, column_tag){  typedef typename matrix_traits<MatA>::size_type Int;  typename MatA::const_iterator A_k;  typename MatA::OneD::const_iterator A_ki, A_kiend;  A_k = A.begin();  while (not_at(A_k, A.end())) {