// -*- c++ -*-//// Copyright 1997, 1998, 1999 University of Notre Dame.// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee//// This file is part of the Matrix Template Library//// You should have received a copy of the License Agreement for the// Matrix Template Library along with the software;  see the// file LICENSE.  If not, contact Office of Research, University of Notre// Dame, Notre Dame, IN  46556.//// Permission to modify the code and to distribute modified code is// granted, provided the text of this NOTICE is retained, a notice that// the code was modified is included with the above COPYRIGHT NOTICE and// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE// file is distributed with the modified code.//// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.// By way of example, but not limitation, Licensor MAKES NO// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS// OR OTHER RIGHTS.////===========================================================================#ifndef _MTL_MTL_H_#define _MTL_MTL_H_#include <functional>#include <iostream>#include "mtl/mtl_limits.h"#include "mtl/mtl_complex.h"#include "mtl/fast.h"#include "mtl/dense1D.h"#include "mtl/mtl_exception.h"#include "mtl/matrix_traits.h"#include "mtl/transform_iterator.h"#include "mtl/scaled1D.h"#include "mtl/abs.h"#ifdef USE_DOUBLE_DOUBLE#include "contrib/double_double/double_double.h"#endif#if USE_BLAIS#include "mtl/blais.h"#endif#include "mtl/matrix.h"/*  This is a nasty hack necessitated by several things:  1. C++ does not allow temporaries to be passed    into a reference argument.  2. Many MTL expressions result in temporaries  3. Some MTL matrix classes (static matrix) can    not be passed by value for the output argument    since they are not handles. */#define MTL_OUT(X) const X&namespace mtl {template <class T>inline T sign(const T& x) { return (x < 0) ? T(-1) : T(1); }template <class T>inline T xfer_sign(const T& x, const T& y){  return (y < 0) ? -MTL_ABS(x) : MTL_ABS(x);}//: for tri_solve and others//!noindex:class right_side { };//: for tri_solve and others//!noindex:class left_side { };#include "mtl/dim_calc.h"template <class Vector> inlinetypename linalg_traits<Vector>::value_typesum__(const Vector& x, fast::count<0>){  typedef typename linalg_traits<Vector>::value_type vt;  return mtl_algo::accumulate(x.begin(), x.end(), vt());}#if USE_BLAIStemplate <class Vector, int N> inlinetypename linalg_traits<Vector>::value_typesum__(const Vector& x, fast::count<N>){  typedef typename linalg_traits<Vector>::value_type vt;  return fast::accumulate(x.begin(), fast::count<N>(), vt());}#endif//: Sum:  <tt>s <- sum_i(x(i))</tt>//!category: algorithms//!component: function//!definition: mtl.h//!example: vec_sum.cc//!complexity: linear//!typereqs: The addition operator must be defined for <TT>Vector::value_type</TT>.// The sum of all of the elements in the container.template <class Vector> inlinetypename linalg_traits<Vector>::value_typesum(const Vector& x){  return sum__(x, dim_n<Vector>::RET());}#include "mtl/mtl_set.h"template <class S, class T, class R>struct mtl_multiplies : std::binary_function<S, T, R> {  typedef S first_argument_type;  typedef T second_argument_type;  typedef R result_type;  R operator () (const S& x, const T& y) const { return x * y; }};template <class Vector, class T> inlinevoidoned_scale(Vector& x, const T& alpha, fast::count<0>){  typedef typename Vector::value_type VT;  mtl_algo::transform(x.begin(), x.end(), x.begin(),                      std::bind1st(mtl_multiplies<T,VT,VT>(), alpha));}#if USE_BLAIStemplate <class Vector, class T, int N> inlinevoidoned_scale(Vector& x, const T& alpha, fast::count<N>){  typedef typename Vector::value_type VT;  fast::transform(x.begin(), fast::count<N>(), x.begin(),                  std::bind1st(mtl_multiplies<T,VT,VT>(), alpha));}#endiftemplate <class Vector, class T> inlinevoidscale_dim(Vector& x, const T& alpha, oned_tag){  oned_scale(x, alpha, dim_n<Vector>::RET());}template <class Matrix, class T>inline voidscale_dim(Matrix& A, const T& alpha, twod_tag){  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 *= alpha;  }}//: Scale:  <tt>A <- alpha*A or x <- alpha x</tt>//// Multiply all the elements in <tt>A</tt> (or <tt>x</tt>) by// <tt>alpha</tt>.// //!category: algorithms//!component: function//!example: vec_scale_algo.cc//!complexity: O(n)//!definition: mtl.h//!typereqs: <TT>Vector</TT> must be mutable//!typereqs: <TT>T</TT> is convertible to <TT>Vector</TT>'s <TT>value_type</TT>//!typereqs: The multiplication operator must be defined for <TT>Vector::value_type</TT> and <tt>T</tt>template <class LinalgObj, class T>inline voidscale(MTL_OUT(LinalgObj) A, const T& alpha){  typedef typename linalg_traits<LinalgObj>::dimension Dim;  scale_dim(const_cast<LinalgObj&>(A), alpha, Dim());  }//: Set Diagonal:  <tt>A(i,i) <- alpha</tt>//// Set the value of the elements on the main diagonal of A to alpha.////!category: algorithms//!component: function//!definition: mtl.h//!example: tri_pack_sol.cc//!typereqs: <tt>T</tt> must be convertible to <tt>Matrix::value_type</tt>.//!complexity: O(min(m,n)) for dense matrices, O(nnz) for sparse matrices (except envelope, which is O(m))template <class Matrix, class T>inline voidset_diagonal(MTL_OUT(Matrix) A_, const T& alpha){  Matrix& A = const_cast<Matrix&>(A_);  typedef typename mtl::matrix_traits<Matrix>::size_type Int;  if (! A.is_unit())    for (Int i = 0; i < A.nrows() && i < A.ncols(); ++i)      A(i,i) = alpha;}//: add absolute value//!noindex:struct abs_add {  template <class T, class U>  T operator()(const T& a, const U& b) {    return a + MTL_ABS(b);  }};template <class Vector>inline typename linalg_traits<Vector>::magnitude_typeoned_one_norm(const Vector& x, fast::count<0>){  typedef typename linalg_traits<Vector>::magnitude_type T;  return mtl_algo::accumulate(x.begin(), x.end(), T(), abs_add());}#if USE_BLAIStemplate <class Vector, int N>inline typename linalg_traits<Vector>::magnitude_typeoned_one_norm(const Vector& x, fast::count<N>){  typedef typename     number_traits<typename Vector::value_type>::magnitude_type T;  return fast::accumulate(x.begin(), fast::count<N>(), T(), abs_add());}#endiftemplate <class Vector>inline typename linalg_traits<Vector>::magnitude_typeone_norm(const Vector& x, oned_tag){  return oned_one_norm(x, dim_n<Vector>::RET());}//: add square//!noindex:struct sqr_add {   template <class T, class U>  T operator()(const T& a, const U& b) {    return a + MTL_ABS(b * b);  }};template <class Vector>inline typename linalg_traits<Vector>::magnitude_typeoned_two_norm(const Vector& x, fast::count<0>){  typedef typename Vector::value_type T;  typedef typename number_traits<T>::magnitude_type M;  using std::sqrt;  return ::sqrt(mtl_algo::accumulate(x.begin(), x.end(), M(), sqr_add()));}#if USE_BLAIStemplate <class Vector, int N>inline typename linalg_traits<Vector>::magnitude_typeoned_two_norm(const Vector& x, fast::count<N>){  typedef typename Vector::value_type T;  typedef typename number_traits<T>::magnitude_type M;  using std::sqrt;  return ::sqrt(fast::accumulate(x.begin(), fast::count<N>(), M(), sqr_add()));}#endif//: Two Norm: <tt>s <- sqrt(sum_i(|x(i)^2|))</tt>////  The square root of the sum of the squares of the elements of the container.////!category: algorithms//!component: function//!definition: mtl.h//!example: vec_two_norm.cc//!complexity: O(n)//!typereqs: <tt>Vector</tt> must have an associated magnitude_type that is the type of the absolute value of <tt>Vector::value_type</tt>.//!typereqs: There must be <tt>abs()</tt> defined for <tt>Vector::value_type</tt>.//!typereqs: The addition must be defined for magnitude_type.//!typereqs: <tt>sqrt()</tt> must be defined for magnitude_type.template <class Vector>inline typename linalg_traits<Vector>::magnitude_typetwo_norm(const Vector& x){  return oned_two_norm(x, dim_n<Vector>::RET());}//: add square//!noindex:struct sqr_ {   template <class T, class U>  T operator()(const T& a, const U& b) {    return a + MTL_ABS(b * b);  }};//: Sum of the Squares////!category: algorithms//!component: function//!definition: mtl.h//!complexity: O(n)template <class Vector>inline typename linalg_traits<Vector>::value_typesum_squares(const Vector& x){  typedef typename linalg_traits<Vector>::value_type T;  return mtl_algo::accumulate(x.begin(), x.end(), T(), sqr_add());}//: compare absolute values//!noindex:struct abs_cmp { template <class T>bool operator()(const T& a, const T& b) {  return MTL_ABS(a) < MTL_ABS(b);}};template <class Vec>inline typename linalg_traits<Vec>::magnitude_typeinfinity_norm(const Vec& x, oned_tag){  return MTL_ABS(*mtl_algo::max_element(x.begin(), x.end(), abs_cmp()));}//: use by one and inf norm//!noindex:template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typemajor_norm__(const Matrix& A){  typedef typename linalg_traits<Matrix>::magnitude_type T;  typedef typename matrix_traits<Matrix>::size_type Int;  T norm = 0;  T sum = 0;  typename Matrix::const_iterator i;  typename Matrix::OneD::const_iterator j;  i = A.begin();  /* get the first sum */  if (i != A.end()) {    j = (*i).begin();    sum = T(0);    for (; j != (*i).end(); ++j)      sum = sum + MTL_ABS(*j);    norm = sum;    ++i;  }  for (; i != A.end(); ++i) {    j = (*i).begin();    if (A.is_unit() && Int(i.index()) < MTL_MIN(A.nrows(), A.ncols()))      sum = T(1);    else sum = T(0);    for (; j != (*i).end(); ++j)      sum = sum + MTL_ABS(*j);    norm = MTL_MAX(MTL_ABS(norm), MTL_ABS(sum));  }  return norm;}//: used by one and inf norm//!noindex:template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typeminor_norm__(const Matrix& A){  typedef typename linalg_traits<Matrix>::magnitude_type T;  typedef typename matrix_traits<Matrix>::size_type Int;  typename Matrix::const_iterator i;  typename Matrix::OneD::const_iterator j, jend;  dense1D<T> sums(A.minor(), T());  if (A.is_unit()) {    for (Int x = 0; x < MTL_MIN(A.nrows(), A.ncols()); ++x)      sums[x] = T(1);  }  for (i = A.begin(); i != A.end(); ++i) {    j = (*i).begin(); jend = (*i).end();    for (; j != jend; ++j)      sums[j.index()] += MTL_ABS(*j);  }  return infinity_norm(sums, oned_tag());}/* this handles both the major and minor norm for symmetric matrices */template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typesymmetric_norm(const Matrix& A, row_tag){  typedef typename linalg_traits<Matrix>::magnitude_type T;  typename Matrix::const_iterator i;  typename Matrix::OneD::const_iterator j, jend;    dense1D<T> sums(A.minor(), T(0));  for (i = A.begin(); i != A.end(); ++i) {    j = (*i).begin();    jend = (*i).end();    if (A.is_upper()) { /* handle the diagonal elements */      sums[j.row()] += MTL_ABS(*j);      ++j;    } else      --jend;    for (; j != jend; ++j) {      sums[j.row()] += MTL_ABS(*j);      sums[j.column()] += MTL_ABS(*j);    }    if (A.is_lower())      sums[j.row()] += MTL_ABS(*j);  }  return infinity_norm(sums, oned_tag());}template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typesymmetric_norm(const Matrix& A, column_tag){  typedef typename linalg_traits<Matrix>::magnitude_type T;  typename Matrix::const_iterator i;  typename Matrix::OneD::const_iterator j, jend;    dense1D<T> sums(A.minor(), T(0));  for (i = A.begin(); i != A.end(); ++i) {    j = (*i).begin();    jend = (*i).end();    if (A.is_lower()) { /* handle the diagonal elements */      sums[j.row()] += MTL_ABS(*j);      ++j;    } else      --jend;    for (; j != jend; ++j) {      sums[j.row()] += MTL_ABS(*j);      sums[j.column()] += MTL_ABS(*j);    }    if (A.is_upper())      sums[j.row()] += MTL_ABS(*j);  }  return infinity_norm(sums, oned_tag());}template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typesymmetric_norm(const Matrix& A){  typedef typename matrix_traits<Matrix>::orientation Orien;  return symmetric_norm(A, Orien());}template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typediagonal_one_norm(const Matrix& A){  typedef typename linalg_traits<Matrix>::magnitude_type T;  typename Matrix::const_iterator i;  typename Matrix::OneD::const_iterator j, jend;