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 voidmatmat_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//!noindextemplate <class MatA, class MatB, class MatC>inline voidmatmat_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 voidmatmat_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 voidmatmat_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 voidmult_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 voidmult(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 voidtri_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 <leopardi@bigpond.net.au> 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 voidtri_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 voidtri_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 voidtri_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 voidtri_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 voidtri_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 voidrank_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 voidrank_two_update(MTL_OUT(Matrix) A_,		const VecX& x, const VecY& y) MTL_THROW_AS