mtl.h

MTL C++ Numeric Library

  typename Matrix::OneD::const_iterator j, jend;

  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_type
diagonal_infinity_norm(const Matrix& A)
{
  typedef 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_type
one_norm__(const Matrix& A, column_tag)
{
  return major_norm__(A);
}


//: dispatch function
//!noindex:
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
one_norm__(const Matrix& A, row_tag)
{
  return minor_norm__(A);
}


template <class Matrix, class Shape>
inline typename linalg_traits<Matrix>::magnitude_type
twod_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_type
twod_one_norm(const Matrix& A, symmetric_tag)
{
  return symmetric_norm(A);
}

template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
twod_one_norm(const Matrix& A, diagonal_tag)
{
  return diagonal_one_norm(A);
}


template <class Linalg>
inline typename linalg_traits<Linalg>::magnitude_type
one_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_type
one_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_type
infinity_norm__(const Matrix& A, row_tag)
{
  return major_norm__(A);
}

//: dispatch function
//!noindex:
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
infinity_norm__(const Matrix& A, column_tag)
{
  return minor_norm__(A);
}

template <class Matrix, class Shape>
inline typename linalg_traits<Matrix>::magnitude_type
twod_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_type
twod_infinity_norm(const Matrix& A, symmetric_tag)
{
  return symmetric_norm(A);
}

template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
twod_infinity_norm(const Matrix& A, diagonal_tag)
{
  return diagonal_infinity_norm(A);
}


template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
infinity_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_type
infinity_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_type
max_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_type
max_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_type
min_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_type
min_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 value
template <class VectorT>
inline typename VectorT::value_type
max(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_type
min(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_; }
#endif
protected:
#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: type
template <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)) {