dense1D<T> sums(A.ncols(), T(0));  for (i = A.begin(); i != A.end(); ++i) {    j = (*i).begin(); jend = (*i).end();    for (; j != jend; ++j)      sums[j.column()] += MTL_ABS(*j);  }  return infinity_norm(sums);}template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typediagonal_infinity_norm(const Matrix& A){  typedef typename linalg_traits<Matrix>::magnitude_type T;  typename Matrix::const_iterator i;  typename Matrix::OneD::const_iterator j, jend;  dense1D<T> sums(A.nrows(), T(0));  for (i = A.begin(); i != A.end(); ++i) {    j = (*i).begin(); jend = (*i).end();    for (; j != jend; ++j)      sums[j.row()] += MTL_ABS(*j);  }  return infinity_norm(sums);}//: dispatch function//!noindex:template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typeone_norm__(const Matrix& A, column_tag){  return major_norm__(A);}//: dispatch function//!noindex:template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typeone_norm__(const Matrix& A, row_tag){  return minor_norm__(A);}template <class Matrix, class Shape>inline typename linalg_traits<Matrix>::magnitude_typetwod_one_norm(const Matrix& A, Shape){  typedef typename Matrix::orientation Orien;  return one_norm__(A, Orien());}template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typetwod_one_norm(const Matrix& A, symmetric_tag){  return symmetric_norm(A);}template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typetwod_one_norm(const Matrix& A, diagonal_tag){  return diagonal_one_norm(A);}template <class Linalg>inline typename linalg_traits<Linalg>::magnitude_typeone_norm(const Linalg& A, twod_tag){  typedef typename matrix_traits<Linalg>::shape Shape;  return twod_one_norm(A, Shape());}//: One Norm:  <tt>s <- sum(|x_i|) or s <- max_i(sum_j(|A(i,j)|))</tt>//// For vectors, the sum of the absolute values of the elements.// For matrices, the maximum of the column sums.// Note: not implemented yet for unit triangle matrices.////!category: algorithms//!component: function//!definition: mtl.h//!example: vec_one_norm.cc//!complexity: O(n)//!typereqs: The vector or matrix must have an associated magnitude_type that//   is the type of the absolute value of its <tt>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.template <class LinalgObj>inline typename linalg_traits<LinalgObj>::magnitude_typeone_norm(const LinalgObj& A){  typedef typename linalg_traits<LinalgObj>::dimension Dim;  return one_norm(A, Dim());}//: dispatch function//!noindex:template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typeinfinity_norm__(const Matrix& A, row_tag){  return major_norm__(A);}//: dispatch function//!noindex:template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typeinfinity_norm__(const Matrix& A, column_tag){  return minor_norm__(A);}template <class Matrix, class Shape>inline typename linalg_traits<Matrix>::magnitude_typetwod_infinity_norm(const Matrix& A, Shape){  typedef typename Matrix::orientation Orien;  return infinity_norm__(A, Orien());}template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typetwod_infinity_norm(const Matrix& A, symmetric_tag){  return symmetric_norm(A);}template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typetwod_infinity_norm(const Matrix& A, diagonal_tag){  return diagonal_infinity_norm(A);}template <class Matrix>inline typename linalg_traits<Matrix>::magnitude_typeinfinity_norm(const Matrix& A, twod_tag){  typedef typename matrix_traits<Matrix>::shape Shape;  return twod_infinity_norm(A, Shape());}//: Infinity Norm: <tt>s <- max_j(sum_i(|A(i,j)|)) or s <- max_i(|x(i)|)</tt>//// For matrices, the maximum of the row sums.// For vectors, the maximum absolute value of any of its element.////!category: algorithms//!component: function//!definition: mtl.h//!complexity: O(n) for vectors, O(m*n) for dense matrices, O(nnz) for sparse//!example: vec_inf_norm.cc//!typereqs: The vector or matrix must have an associated magnitude_type that is the type of the absolute value of its <tt>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.template <class LinalgObj>inline typename linalg_traits<LinalgObj>::magnitude_typeinfinity_norm(const LinalgObj& A){  typedef typename linalg_traits<LinalgObj>::dimension Dim;  return infinity_norm(A, Dim());}//: Max Index:  <tt>i <- index of max(|x(i)|)</tt>//!category: algorithms//!component: function//!definition: mtl.h//!complexity: O(n)// The location (index) of the element with the maximum absolute value.//!example: max_index.cc//!typereqs: <tt>Vec::value_type</tt> must be LessThanComparible.template <class Vec>inline typename Vec::size_typemax_index(const Vec& x){  typename Vec::const_iterator maxi =    mtl_algo::max_element(x.begin(), x.end(), abs_cmp());  return maxi.index();}//: Maximum Absolute Index:  <tt>i <- index of max(|x(i)|)</tt>//!category: algorithms//!component: function//!definition: mtl.h//!complexity: O(n)// The location (index) of the element with the maximum absolute value.//!example: max_abs_index.cc//!typereqs: The vector or matrix must have an associated magnitude_type that//   is the type of the absolute value of its <tt>value_type</tt>.//!typereqs: There must be <tt>abs()</tt> defined for <tt>Vector::value_type</tt>.//!typereqs: The magnitude type must be LessThanComparible.template <class Vec>inline typename Vec::size_typemax_abs_index(const Vec& x){  typename Vec::const_iterator maxi =    mtl_algo::max_element(x.begin(), x.end(), abs_cmp());  return maxi.index();}//: Minimum Index:  <tt>i <- index of min(x(i))</tt>//!category: algorithms//!component: function//!definition: mtl.h//!complexity: O(n)// The location (index) of the element with the minimum value.//!example: min_abs_index.cc//!typereqs: <tt>Vec::value_type</tt> must be LessThanComparible.template<class Vec>inline typename Vec::size_typemin_index(const Vec& x) {  typename Vec::const_iterator mini =     mtl_algo::min_element(x.begin(), x.end());     return mini.index(); }                    //: Minimum Absolute Index:  <tt>i <- index of min(|x(i)|)</tt>//!category: algorithms//!component: function//!definition: mtl.h//!complexity: O(n)// The location (index) of the element with the minimum absolute value.//!example: max_index.cc//!typereqs: The vector or matrix must have an associated magnitude_type that//   is the type of the absolute value of its <tt>value_type</tt>.//!typereqs: There must be <tt>abs()</tt> defined for <tt>Vector::value_type</tt>.//!typereqs: The magnitude type must be LessThanComparible.template<class Vec>inline typename Vec::size_typemin_abs_index(const Vec& x) {  typename Vec::const_iterator mini =     mtl_algo::min_element(x.begin(), x.end(), abs_cmp());     return mini.index(); }                    //: Max Value:  <tt>s <- max(x(i))</tt>//!category: algorithms//!component: function//!definition: mtl.h//!example: vec_max.cc//!complexity: O(n)//!typereqs: <tt>Vec::value_type</tt> must be LessThanComparible.// Returns the value of the element with the maximum valuetemplate <class VectorT>inline typename VectorT::value_typemax(const VectorT& x){  return *mtl_algo::max_element(x.begin(), x.end());}//: Min Value:  <tt>s <- min(x_i)</tt>//!category: algorithms//!component: function//!complexity: O(n)//!definition: mtl.h//!typereqs: <tt>Vec::value_type</tt> must be LessThanComparible.template <class VectorT>inline typename VectorT::value_typemin(const VectorT& x){  return *mtl_algo::min_element(x.begin(), x.end());}#define MTL_BLAS_GROT//use blas version always, since there is a bug in the lapack verions // of givens_rotation according to Andy's email//: Givens Plane Rotation//!category: functors//!component: type//!definition: mtl.h//!example: apply_givens.cc//// Input a and b to the constructor to create a givens plane rotation// object. Then apply the rotation to two vectors. There is a// specialization of the givens rotation for complex numbers.//// <codeblock>// [  c  s ] [ a ] = [ r ]// [ -s  c ] [ b ]   [ 0 ]// </codeblock>////!typereqs: the addition operator must be defined for <tt>T</tt>//!typereqs: the multiplication operator must be defined for <tt>T</tt>//!typereqs: the division operator must be defined for <tt>T</tt>//!typereqs: the abs() function must be defined for <tt>T</tt>template <class T>class givens_rotation {public:  //: Default constructor  inline givens_rotation()     : #ifdef MTL_BLAS_GROT    a_(0), b_(0),#endif    c_(0), s_(0)#ifndef MTL_BLAS_GROT    , r_(0) #endif  { }  //: Givens Plane Rotation Constructor  inline givens_rotation(T a_in, T b_in) {#ifdef MTL_BLAS_GROT // old BLAS version    T roe;    if (MTL_ABS(a_in) > MTL_ABS(b_in))      roe = a_in;    else      roe = b_in;        T scal = MTL_ABS(a_in) + MTL_ABS(b_in);    T r, z;    if (scal != T(0)) {      T a_scl = a_in / scal;      T b_scl = b_in / scal;      r = scal * sqrt(a_scl * a_scl + b_scl * b_scl);      if (roe < T(0)) r *= -1;      c_ = a_in / r;      s_ = b_in / r;      z = 1;      if (MTL_ABS(a_in) > MTL_ABS(b_in))        z = s_;      else if (MTL_ABS(b_in) >= MTL_ABS(a_in) && c_ != T(0))        z = T(1) / c_;    } else {      c_ = 1; s_ = 0; r = 0; z = 0;          }    a_ = r;    b_ = z;#else // similar LAPACK slartg version, modified to the NEW BLAS proposal    T a = a_in, b = b_in;    if (b == T(0)) {      c_ = T(1);      s_ = T(0);      r_ = a;    } else if (a == T(0)) {      c_ = T(0);      s_ = sign(b);      r_ = b;    } else {      // cs = |a| / sqrt(|a|^2 + |b|^2)      // sn = sign(a) * b / sqrt(|a|^2 + |b|^2)      T abs_a = MTL_ABS(a);      T abs_b = MTL_ABS(b);      if (abs_a > abs_b) {        // 1/cs = sqrt( 1 + |b|^2 / |a|^2 )        T t = abs_b / abs_a;        T tt = sqrt(T(1) + t * t);        c_ = T(1) / tt;        s_ = t * c_;        r_ = a * tt;      } else {        // 1/sn = sign(a) * sqrt( 1 + |a|^2/|b|^2 )        T t = abs_a / abs_b;        T tt = sqrt(T(1) + t * t);        s_ = sign(a) / tt;        c_ = t * s_;        r_ = b * tt;      }    }#endif  }  inline void set_cs(T cin, T sin) { c_ = cin; s_ = sin; }  //: Apply plane rotation to two real scalars. (name change a VC++ workaround)  inline void scalar_apply(T& x, T& y) {    T tmp = c_ * x + s_ * y;    y = c_ * y - s_ * x;    x = tmp;  }  //: Apply plane rotation to two vectors.  template <class VecX, class VecY>  inline void apply(MTL_OUT(VecX) x_, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION {    VecX& x = const_cast<VecX&>(x_);    VecY& y = const_cast<VecY&>(y_);    MTL_ASSERT(x.size() <= y.size(), "mtl::givens_rotation::apply()");    typename VecX::iterator xi = x.begin();    typename VecX::iterator xend = x.end();    typename VecY::iterator yi = y.begin();    while (mtl::not_at(xi, xend)) {      scalar_apply(*xi, *yi);      ++xi; ++yi;    }  }#ifdef MTL_BLAS_GROT  inline T a() { return a_; }  inline T b() { return b_; }#endif  inline T c() { return c_; }  inline T s() { return s_; }#ifndef MTL_BLAS_GROT  inline T r() { return r_; }#endifprotected:#ifdef MTL_BLAS_GROT  T a_, b_;#endif  T c_, s_;#ifndef MTL_BLAS_GROT  T r_;#endif};using std::real;using std::imag;#if MTL_PARTIAL_SPEC//:  The specialization for complex numbers.//!category: functors//!component: typetemplate <class T>class givens_rotation < std::complex<T> > {  typedef std::complex<T> C;public:  //:  inline givens_rotation() : cs(0), sn(0)#ifndef MTL_BLAS_GROT    , r_(0)#endif  { }    inline T abs_sq(C t) { return real(t) * real(t) + imag(t) * imag(t); }  inline T abs1(C t) { return MTL_ABS(real(t)) + MTL_ABS(imag(t)); }  //:  inline givens_rotation(C a_in, C b_in) {#ifdef MTL_BLAS_GROT    T a = std::abs(a_in), b = std::abs(b_in);    if ( a == T(0) ) {      cs = T(0);      sn = C(1.);      //in zrotg there is an assignment for ca, what is that for?     } else {      T scale = a + b;      T norm = std::sqrt(abs_sq(a_in/scale)+abs_sq(b_in/scale)) * scale;          cs = a / norm;      sn = a_in/a * std::conj(b_in)/norm;      //in zrotg there is an assignment for ca, what is that for?     }#else // LAPACK version, clartg    C f(a_in), g(b_in);    if (g == C(0)) {      cs = T(1);      sn = C(0);      r_ = f;    } else if (f == C(0)) {      cs = T(0);      sn = MTL_CONJ(g) / MTL_ABS(g);