mtl.h

很好用的库

    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_lower()) { /* 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_upper())
        C(i,j) += *A_ki * B(k,j);
    }

    ++A_k;
  }
}


//: Specialization for triangular matrices
//!noindex:
template <class MatA, class MatB, class MatC>
inline void
matmat_mult(const MatA& A, const MatB& B, MatC& C, symmetric_tag)
{
  typedef typename matrix_traits<MatA>::orientation Orien;
  symm_simple_mult(A, B, C, Orien());
}

//: Specialization for triangular matrices
//!noindex
template <class MatA, class MatB, class MatC>
inline void
matmat_mult(const MatA& A, const MatB& B, MatC& C, triangle_tag)
{
  typedef typename matrix_traits<MatA>::size_type Int;
  typedef typename matrix_traits<MatA>::orientation Orien;
  if (A.is_unit()) {
    Int M = MTL_MIN(A.nrows(), A.ncols());
    Int N = B.ncols();
    for (Int i = 0; i < M; ++i)
      for (Int j = 0; j < N; ++j)
        C(i,j) += B(i,j);
  }

  simple_mult(A, B, C, mtl::dense_tag(), Orien());
}

//: Dispatch to row/column general and banded matrices
//!noindex:
template <class MatA, class MatB, class MatC>
inline void
matmat_mult(const MatA& A, const MatB& B, MatC& C, rectangle_tag)
{
  typedef typename matrix_traits<MatA>::sparsity Sparsity;
  typedef typename matrix_traits<MatA>::orientation Orien;
  simple_mult(A, B, C, Sparsity(), Orien());
}

template <class MatA, class MatB, class MatC>
inline void
matmat_mult(const MatA& A, const MatB& B, MatC& C, banded_tag)
{
  typedef typename matrix_traits<MatC>::sparsity Sparsity;
  typedef typename matrix_traits<MatA>::orientation Orien;
  simple_mult(A, B, C, Sparsity(), Orien());
}


//: Matrix multiplication  C <- C + A * B
//
//  The actual specialization of the algorithm used depends of the
//  types of matrices used. If all the matrices are dense and
//  rectangular the blocked algorithm is used (when --with-blais is
//  specified in the configure). Otherwise the traversal depends on
//  matrix A. Therefore if one is multiplying a sparse matrix by a
//  dense, one would want the sparse matrix as the A
//  argument. Typically, for performance reasons, one would not want
//  to use a sparse matrix for C.
//  <p>
//  Note: ignore the <tt>twod_tag</tt> argument and the underscores in
//  the name of this function.
//
//!precond: <tt>A.nrows() == C.nrows()</tt>
//!precond: <tt>A.ncols() == B.nrows()</tt>
//!precond: <tt>B.ncols() == C.ncols()</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!typereqs: the value types for each of the matrices must be compatible
//!typereqs: the multiplication operator must be defined for <tt>MatA::value_type</tt>
//!typereqs: the addition operator must be defined for <tt>MatA::value_type</tt>
template <class MatA, class MatB, class MatC>
inline void
mult_dim__(const MatA& A, const MatB& B, MatC& C, twod_tag)
{
  typedef typename MatA::shape Shape;
  matmat_mult(A, B, C, Shape());
}


//: Dispatch between matrix matrix and matrix vector mult.
//!noindex:
template <class LinalgA, class LinalgB, class LinalgC>
inline void
mult(const LinalgA& A, const LinalgB& B, MTL_OUT(LinalgC) C_)
{
  LinalgC& C = const_cast<LinalgC&>(C_);
  typedef typename linalg_traits<LinalgB>::dimension Dim;
  mult_dim__(A, B, C, Dim());
}

//: for column oriented
//!noindex:
template <class TriMatrix, class VecX>
inline void
tri_solve__(const TriMatrix& T, VecX& x, column_tag)
{
  typedef typename matrix_traits<TriMatrix>::size_type Int;
  typedef typename matrix_traits<TriMatrix>::value_type VT;
  typename VecX::value_type x_j; 

  if (T.is_upper()) {
    typename TriMatrix::const_reverse_iterator T_j; 
    typename TriMatrix::Column::const_reverse_iterator T_ji, T_jrend;

    for (T_j = T.rbegin(); T_j != T.rend(); ++T_j) {
      T_ji = (*T_j).rbegin();
      T_jrend = (*T_j).rend();
      //Int j = T_ji.column();
      Int j = T_j.index();
      
      //Paul C. Leopardi <[e-mail]> reported the fix 
      //for for a sparse matrix to have a completely empty row (or column)
      if ( (T_ji != T_jrend) && ! T.is_unit()) {
        x[j] /= *T_ji; /* the diagonal */
        ++T_ji;
      }
      x_j = x[j];

      while (T_ji != T_jrend) {
        Int i = T_ji.row();
        x[i] -= x_j * *T_ji;
        ++T_ji;
      }
    }
  } else {                      /* T is lower */
    typename TriMatrix::const_iterator T_j; 
    typename TriMatrix::Column::const_iterator T_ji, T_jend;

    for (T_j = T.begin(); T_j != T.end(); ++T_j) {
      T_ji = (*T_j).begin();
      T_jend = (*T_j).end();
      //Int j = T_ji.column(); //T_ji could be T_jend
      Int j = T_j.index();
      
      if ( (T_ji != T_jend) && ! T.is_unit()) {
        x[j] /= *T_ji; /* the diagonal */
        ++T_ji;
      }
      x_j = x[j];
      
      while (T_ji != T_jend) {
        Int i = T_ji.row();
        x[i] -= x_j * *T_ji;
        ++T_ji;
      }
    }
  }    
}

//: for row major
//!noindex:
template <class TriMatrix, class VecX>
inline void
tri_solve__(const TriMatrix& T, VecX& x, row_tag)
{
  typedef typename matrix_traits<TriMatrix>::value_type VT;
  typedef typename matrix_traits<TriMatrix>::size_type Int;

  if (T.is_upper()) {
    typename TriMatrix::const_reverse_iterator T_i, T_iend; 
    typename TriMatrix::Row::const_reverse_iterator T_ij;

    T_i = T.rbegin();
    T_iend = T.rend();

    if ( (T_i != T_iend) && ! T.is_unit()) {
      T_ij = (*T_i).rbegin();
      x[T_ij.row()] /= *T_ij;
      ++T_i;
    }

    while (T_i != T_iend) {
      T_ij = (*T_i).rbegin();
      //Int i = T_ij.row();
      Int i = T_i.index();
      VT t = x[i];

      typename TriMatrix::Row::const_reverse_iterator T_iend;
      T_iend = (*T_i).rend();
      if ( (T_ij != T_iend) && ! T.is_unit())
        --T_iend;

      Int j;
      while (T_ij != T_iend) {
        j = T_ij.column();
        t -= (*T_ij) * x[j];      
        ++T_ij;
      }
      if ( (*T_i).rbegin() != (*T_i).rend() && !T.is_unit()) //T_i is not empty
        t /= *T_ij;
        
      x[i] = t;

      ++T_i;
    }
  } else { /* T is lower */

    typename TriMatrix::const_iterator T_i; 
    typename TriMatrix::Row::const_iterator T_ij;

    T_i = T.begin();

    if (T_i != T.end() && ! T.is_unit()) {
      T_ij = (*T_i).begin();
      x[T_ij.row()] *= VT(1) / *T_ij;
      ++T_i;
    }

    while (T_i != T.end()) {
      T_ij = (*T_i).begin();
      //Int i = T_ij.row(); //T_ij could be bad
      Int i = T_i.index();
      VT t = x[i];

      typename TriMatrix::Row::const_iterator T_iend;
      T_iend = (*T_i).end();
      if ( ( T_ij != T_iend ) &&  ! T.is_unit())
        --T_iend;

      Int j;
      while (T_ij != T_iend) {
        j = T_ij.column();
        t -= (*T_ij) * x[j];
        ++T_ij;
      }
      if ( (*T_i).begin() !=(*T_i).end() &&  !T.is_unit())
        t /= *T_ij;

      x[i] = t;
      ++T_i;
    }
  }
}


//: Triangular Solve:  <tt>x <- T^{-1} * x</tt>
//  Use with trianguler matrixes only ie. use the <TT>triangle</TT>
//  adaptor class.
//
//  To use with a sparse matrix, the sparse matrix must be wrapped with
//  a triangle adaptor. You must specify "packed" in the triangle
//  adaptor. The sparse matrix must only have elements in the correct
//  side.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: tri_solve.cc
//!typereqs: <tt>Matrix::value_type</tt> and <tt>VecX::value_type</tt> must be the same type
//!typereqs: the multiplication operator must be defined for <tt>Matrix::value_type</tt>
//!typereqs: the division 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 TriMatrix, class VecX>
inline void
tri_solve(const TriMatrix& T, MTL_OUT(VecX) x_) MTL_THROW_ASSERTION
{
  VecX& x = const_cast<VecX&>(x_);
  MTL_ASSERT(T.nrows() <= x.size(), "mtl::tri_solve()");
  MTL_ASSERT(T.ncols() <= x.size(), "mtl::tri_solve()");
  MTL_ASSERT(T.ncols() == T.nrows(), "mtl::tri_solve()");
  typedef typename TriMatrix::orientation orien;
  tri_solve__(T, x, orien());
}





//: tri solve for left side
//!noindex:
template <class MatT, class MatB>
inline void
tri_solve__(const MatT& T, MatB& B, left_side)
{
  /*  const int M = B.nrows(); */
  const int N = B.ncols();

  /* unoptimized version */
  for (int j = 0; j < B.ncols(); ++j)
    mtl::tri_solve(T, columns(B)[j]);


  /* JGS need to do an optimized version of this
  if (T.is_upper()) {
    for (int k = M-1; k > 0; --k) {
      if (B(k,j) != 0) {
        if (! T.is_unit())
          B(k,j) /= T(k,k);
        for (int i = 0; i < k; ++i)
          B(i,j) -= B(k,j) * T(i,k);
      }
    }
  } else {
    for (int j = 0; j < N; ++j)
      for (int k = 0; k < M; ++k) {
        if (B(k,j) != 0) {
          if (! T.is_unit())
            B(k,j) /= T(k,k);
          for (int i = k; i < M; ++i)
            B(i,j) -= B(k,j) * T(i,k);
        }
      }
  }
  */
}


/* JGS untested!!! */

//: tri solve for right side
//!noindex:
template <class MatT, class MatB>
inline void
tri_solve__(const MatT& T, MatB& B, right_side)
{
  const int M = B.nrows();
  const int N = B.ncols();
  typedef typename MatT::PR PR;

  if (T.is_upper()) {
    for (int j = 0; j < N; ++j) {
      for (int k = 0; k < j; ++k)
        if (T(k,j) != PR(0))
          for (int i = 0; i < M; ++i)
            B(i,j) -=  T(k,j) * B(i,k);
      if (! T.is_unit()) {
        PR tmp = PR(1) / T(j,j);
        for (int i = 1; i < M; ++i)
          B(i,j) = tmp * B(i,j);
      }
    }
  } else { // T is lower
    for (int j = N - 1; j > 0; --j) {
      for (int k = j; k < N; ++k)
        if (T(k,j) != PR(0))
          for (int i = 0; i < M; ++i)
            B(i,j) -=  T(k,j) * B(i,k);
      if (! T.is_unit()) {
        PR tmp = PR(1) / T(j,j);
        for (int i = 1; i < M; ++i)
          B(i,j) = tmp * B(i,j);
      }
    }
  }
}

//: Triangular Solve: <tt>B <- A^{-1} * B  or  B <- B * A^{-1}</tt>
//
//  This solves the equation <tt>T*X = B</tt> or <tt>X*T = B</tt> where T
//  is an upper or lower triangular matrix, and B is a general
//  matrix. The resulting matrix X is written onto matrix B. The first
//  equation is solved if <tt>left_side</tt> is specified. The second
//  equation is solved if <tt>right_side</tt> is specified.
//
//  Currently only works with dense storage format.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n^3)
//!example: matmat_trisolve.cc
//!typereqs: <tt>MatT::value_type</tt> and <tt>MatB::value_type</tt> must be the same type
//!typereqs: the multiplication operator must be defined for <tt>MatT::value_type</tt>
//!typereqs: the division operator must be defined for <tt>MatT::value_type</tt>
//!typereqs: the addition operator must be defined for <tt>MatT::value_type</tt>
template <class MatT, class MatB, class Side>
inline void
tri_solve(const MatT& T, MTL_OUT(MatB) B, Side s)
{
  tri_solve__(T, const_cast<MatB&>(B), s);
}





//: Rank One Update:   <tt>A <- A  +  x * y^T</tt>
//
// Also known as the outer product of two vectors.
// <codeblock>
//       y = [ 1  2  3 ]
//
//     [ 1 ] [ 1  2  3 ]
// x = [ 2 ] [ 2  4  6 ] => A
//     [ 3 ] [ 3  6  9 ]
//     [ 4 ] [ 4  8 12 ]
// </codeblock>
// <p>
// When using this algorithm with a symmetric matrix, x and y
// must be the same vector, or at least have the same values.
// Otherwise the resulting matrix is not symmetric.
//
//!precond:  <TT>A.nrows() <= x.size()</TT>
//!precond:  <TT>A.ncols() <= y.size()</TT>
//!precond: A has rectangle shape and is dense
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: rank_one.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
rank_one_update(MTL_OUT(Matrix) A_, 
		const VecX& x, const VecY& y) MTL_THROW_ASSERTION
{
  Matrix& A = const_cast<Matrix&>(A_);
  MTL_ASSERT(A.nrows() <= x.size(), "mtl::rank_one_update()");
  MTL_ASSERT(A.ncols() <= y.size(), "mtl::rank_one_update()");
  typename Matrix::iterator i;
  typename Matrix::OneD::iterator j, jend;
  for (i = A.begin(); i != A.end(); ++i) {
    j = (*i).begin(); jend = (*i).end();
    for (; j != jend; ++j)
      *j += x[j.row()] * MTL_CONJ(y[j.column()]);
  }
}



/* 1. how will the scaling by alpha work into this
 * 2. is my placement of conj() ok with respect
 *    to both row and column oriented matrices
 * 3. Perhaps split this in two, have diff version for complex
 */

//: Rank Two Update:  <tt>A <- A  +  x * y^T  +  y * x^T</tt>
//
//
//!category: algorithms
//!component: function
//!precond:   <TT>A.nrows() == A.ncols()</TT>
//!precond:   <TT>A.nrows() == x.size()</TT>
//!precond:   <TT>x.size() == y.size()</TT>
//!precond: A has rectangle shape and is dense.
//!definition: mtl.h
//!example: rank_2_symm_sparse.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
rank_two_update(MTL_OUT(Matrix) A_,
		const VecX& x, const VecY& y) MTL_THROW_AS