mtl.h

很好用的库

      r_ = MTL_ABS(g);
    } else {
      C fs, gs, ss, t;
      T d, di, f1, f2, fa, g1, g2, ga;
      f1 = abs1(f);
      g1 = abs1(g);
      if (f1 >= g1) {
        gs = g / f1;
        g2 = abs_sq(gs);
        fs = f / f1;
        f2 = abs_sq(fs);
        d = sqrt(T(1) + g2 / f2);
        cs = T(1) / d;
        sn = MTL_CONJ(gs) * fs * (cs / f2);
        r_ = f * d;
      } else {
        fs = f / g1;
        f2 = abs_sq(fs);
        fa = sqrt(f2);
        gs = g / g1;
        g2 = abs_sq(gs);
        ga = sqrt(g2);
        d = sqrt(T(1) + f2 / g2);
        di = T(1) / d;
        cs = (fa / ga ) * di;
        ss = (MTL_CONJ(gs) * fs) / (fa * ga);
        sn = ss * di;
        r_ = g * ss * d;
      }
    }
#endif
  }
  //:  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;
    }
  }
  //: Apply plane rotation to two complex scalars.
  inline void scalar_apply(C& x, C& y) {
    complex<T> temp  =  MTL_CONJ(cs) * x + MTL_CONJ(sn) * y;
    y = cs * y - sn * x; 
    x = temp;
  }
  inline void set_cs(const T& cs_, const C& sn_) {
    cs = cs_; sn = sn_;
  }

  inline T c() { return cs; }
  inline C s() { return sn; }
#ifndef MTL_BLAS_GROT
  inline C r() { return r_; }
#endif

protected:
  T cs;
  C sn;
#ifndef MTL_BLAS_GROT
  C r_;
#endif
};

#else

//:  The specialization for complex numbers.
//!category: functors
//!component: type

class givens_rotation < std::complex<double> > {
  typedef double 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);
      r_ = MTL_ABS(g);
    } else {
      C fs, gs, ss, t;
      T d, di, f1, f2, fa, g1, g2, ga;
      f1 = abs1(f);
      g1 = abs1(g);
      if (f1 >= g1) {
        gs = g / f1;
        g2 = abs_sq(gs);
        fs = f / f1;
        f2 = abs_sq(fs);
        d = sqrt(T(1) + g2 / f2);
        cs = T(1) / d;
        sn = MTL_CONJ(gs) * fs * (cs / f2);
        r_ = f * d;
      } else {
        fs = f / g1;
        f2 = abs_sq(fs);
        fa = sqrt(f2);
        gs = g / g1;
        g2 = abs_sq(gs);
        ga = sqrt(g2);
        d = sqrt(T(1) + f2 / g2);
        di = T(1) / d;
        cs = (fa / ga ) * di;
        ss = (MTL_CONJ(gs) * fs) / (fa * ga);
        sn = ss * di;
        r_ = g * ss * d;
      }
    }
#endif
  }
  //:  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;
    }
  }
  //: Apply plane rotation to two complex scalars.
  inline void scalar_apply(C& x, C& y) {
    complex<T> temp  =  MTL_CONJ(cs) * x + MTL_CONJ(sn) * y;
    y = cs * y - sn * x; 
    x = temp;
  }
  T c() { return cs; }
  C s() { return sn; }
#ifndef MTL_BLAS_GROT
  inline C r() { return r_; }
#endif
protected:
  T cs;
  C sn;
#ifndef MTL_BLAS_GROT
  C r_;
#endif
};

//:  The specialization for complex numbers.
//!category: functors
//!component: type
class givens_rotation < std::complex<float> > {
  typedef float 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);
      r_ = MTL_ABS(g);
    } else {
      C fs, gs, ss, t;
      T d, di, f1, f2, fa, g1, g2, ga;
      f1 = abs1(f);
      g1 = abs1(g);
      if (f1 >= g1) {
        gs = g / f1;
        g2 = abs_sq(gs);
        fs = f / f1;
        f2 = abs_sq(fs);
        d = sqrt(T(1) + g2 / f2);
        cs = T(1) / d;
        sn = MTL_CONJ(gs) * fs * (cs / f2);
        r_ = f * d;
      } else {
        fs = f / g1;
        f2 = abs_sq(fs);
        fa = sqrt(f2);
        gs = g / g1;
        g2 = abs_sq(gs);
        ga = sqrt(g2);
        d = sqrt(T(1) + f2 / g2);
        di = T(1) / d;
        cs = (fa / ga ) * di;
        ss = (MTL_CONJ(gs) * fs) / (fa * ga);
        sn = ss * di;
        r_ = g * ss * d;
      }
    }
#endif
  }
  //:  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;
    }
  }
  //: Apply plane rotation to two complex scalars.
  inline void scalar_apply(C& x, C& y) {
    complex<T> temp  =  MTL_CONJ(cs) * x + MTL_CONJ(sn) * y;
    y = cs * y - sn * x; 
    x = temp;
  }
  T c() { return cs; }
  C s() { return sn; }
#ifndef MTL_BLAS_GROT
  inline C r() { return r_; }
#endif

protected:
  T cs;
  C sn;
#ifndef MTL_BLAS_GROT
  C r_;
#endif
};
#endif


#undef MTL_BLAS_GROT
//do allow internal macro escape out of the scope

//: Modified Givens Transformation
//!category: functors
//!component: type
//  
//  This class is under construction.  Like the givens rotation class,
//  there will be a real and complex class.
template <class T>
class modified_givens {

};


template <class T>
inline T two_norm3(const T& x, const T& y, const T& z) {
  return sqrt(x*x + y*y + z*z);
}

//: Generate Householder Transform
//
// Ok to alias x and v to the same vector.
// T can be real or complex.
// Equivalent to LAPACK's xLARFG
template <class T, class Vec>
inline void generate_householder(T& alpha, const Vec& x, 
                                 Vec& v, T& tau) MTL_THROW_ASSERTION
{
  MTL_ASSERT(x.size() == v.size(), "mtl::generate_householder");
  typedef typename number_traits<T>::magnitude_type Real;
  typename Vec::subrange_type subx = x(0, x.size() - 1);
  typename Vec::subrange_type subv = v(0, v.size() - 1);
  v[v.size() - 1] = x[x.size() - 1];
  Real xnorm = two_norm(x);
  Real alpha_r = real(alpha);
  Real alpha_i = imag(alpha);
  
  if (xnorm == Real(0) && alpha_i == Real(0))
    tau = T(0); // H = I
  else {
    Real beta = -xfer_sign(two_norm3(alpha_r, alpha_i, xnorm), alpha_r);
    Real safe_min = std::numeric_limits<Real>::min();
    Real r_safe_min = Real(1) / safe_min;
    
    int count = 0;
    while (MTL_ABS(beta) < safe_min) { // xnorm and beta may be inaccurate
      if (count == 0)                  // so scale x and recompute them
        copy(mtl::scaled(subx, r_safe_min), subv);
      else
        scale(subv, r_safe_min);
      beta *= r_safe_min;
      alpha *= r_safe_min;
      ++count;
    }
    if (count != 0) {
      alpha_r = real(alpha);
      alpha_i = imag(alpha);
      xnorm = two_norm(x);
      beta = -xfer_sign(two_norm3(alpha_r, alpha_i, xnorm), alpha_r);
    }
    tau = beta - (alpha / beta);
    alpha = T(1) / (alpha - beta);
    scale(subv, alpha);
    alpha = beta;
    for (int j = 0; j < count; ++j)
      alpha *= safe_min;
  }
}



//: Householder Transform
//
// Constructor does the generation, then call apply
//
template <class T>
class householder_transform {
  typedef typename number_traits<T>::magnitude_type Real;
  typedef dense1D<T> Vec;
public:
  template <class Vec>
  inline householder_transform(const T& alpha, Vec& x, const T& tau)
    : _alpha(alpha), _v(x.size()), _tau(tau) {
    generate_householder(_alpha, x, _v, _tau);
  }

#if 0
  // JGS conj ???
  // Equivalent to LAPACK xLARF
  template <class MatrixC>
  inline void apply(MatrixC& C, right_side) {
    if (_tau != T(0)) {
      Vec w(C.nrows());
      mult(C, _v, w);		       // w <- C * v
      rank_one_update(mtl::caled(w, -_tau),// C <- C - w * v'
		      conj(v), C);
    }
  }
  template <class MatrixC>
  inline void apply(MatrixC& C, left_side) {
    if (_tau != T(0)) {
      Vec w(C.nrows());
      mult(conj(C), _v, w);	       // w <- C' * v
      rank_one_update(mtl::scaled(w,-_tau), // C <- C - w * v'
		      conj(v), C);
    }
  }
#endif
protected:
  T _alpha;
  Vec _v;
  T _tau;
};



//: Transpose in Place:  <tt>A <- A^T</tt>
// Currently this algorithm only applies to square dense matrices
// Plan to include all rectangular dense matrices..
//!category: algorithms
//!component: function
//!definition: mtl.h
template <class Matrix>
inline void
transpose(MTL_OUT(Matrix) A_) MTL_THROW_ASSERTION
{
  Matrix& A = const_cast<Matrix&>(A_);
  MTL_ASSERT(A.nrows() == A.ncols(), "mat::transpose()");
  typedef typename matrix_traits<Matrix>::value_type T;
  typedef typename mtl::matrix_traits<Matrix>::size_type Int;
  for (Int i = 0; i < A.nrows(); ++i)
    for (Int j = i; j < A.ncols(); ++j) {
      T tmp = A(i, j);
      A(i, j) = A(j, i);
      A(j, i) = tmp;
    }
}


//: Transpose: <tt>B <- A^T</tt>
//!precond:  <tt> B(i,j) = 0 & B = A^T </tt>
//
//  When matrix B is banded, it is up to the user to ensure
//  that the bandwidth is sufficient to contain the elements
//  from A^T. If there are elements of A^T that do not
//  fall within the bandwidth, an exception will be thrown.
//  (exception not implemented yet).
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n^2)
template <class MatA, class MatB>
inline void
transpose(const MatA& A, MTL_OUT(MatB) B_) MTL_THROW_ASSERTION
{
  MatB& B = const_cast<MatB&>(B_);
  MTL_ASSERT(A.nrows() <= B.ncols(), "matmat::transpose()");
  MTL_ASSERT(A.ncols() <= B.nrows(), "matmat::transpose()");

  typename MatA::const_iterator i;
  typename MatA::OneD::const_iterator j, jend;

  for (i = A.begin(); i != A.end(); ++i) {
    j = (*i).begin(); jend = (*i).end();
    for (; j != jend; ++j)
      B(j.column(), j.row()) = *j;
  }
}


/*
  This version of the algorithm depends on the compiler
  hoisting the reference of z[j.row()] out of the inner loop
  (for the row major case)
  KCC doesn't do this, and niether does the underlying Sun C
  compiler.

  In order to hoist the reference by hand, I'll have to write
  specializations for column major matrix and for row major matrices.
  While I'm at it I'll the the unrolling stuff too.

*/