mtl.h

很好用的库

/* this is generic
 */
template <class Matrix, class VecX, class VecZ>
inline void
mult_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 void
mult_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 void
rect_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 void
rect_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 oriented
template <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 oriented
template <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 void
mult_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 void
mult_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 void
mult_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 void
mult_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 void
mult_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 void
mult(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 void
mult_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 void
mult_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 void
simple_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 void
simple_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 void
simple_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 void
symm_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 void
symm_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())) {