mv.h

来自「C++ source code for book-C++ and Object 」· C头文件 代码 · 共 1,506 行 · 第 1/4 页

}

template<class T> T Vcr<T>::maxnorm() const {
  T nm = abs(vr[0]);
  for (int i = 1; i < lenth; i++) {
    T vi = abs(vr[i]);
    if (nm < vi) nm = vi;
  }
  return nm;
}

template<class T> 
Vcr<T> operator*(T scalar, const Vcr<T>& vec) {
  Vcr<T> tm(vec.size());
  for (int i = 0; i < vec.size(); i++) tm[i] = scalar*vec[i]; 
  return tm;
}

template<class T> 
Vcr<T> operator*(const Vcr<T>& vec, T scalar) {
  Vcr<T> tm(vec.size());
  for (int i = 0; i < vec.size(); i++) tm[i] = scalar*vec[i]; 
  return tm;
}

template<class T> 
Vcr<T> operator*(const Vcr<T> & vec1, const Vcr<T> & vec2) {
  int size = vec1.size();
  if (size != vec2.size()) error("bad vector sizes in vector multiply");
  Vcr<T> tm(size);
  for (int i = 0; i < size; i++) tm[i] = vec1[i]*vec2[i]; 
  return tm;
}

template<class T> 
Vcr<T> operator/(const Vcr<T> & vec, T scalar) {
  if (scalar == 0) error("divisor is zero in vector-scalar division");
  Vcr<T> tm(vec.size());
  for (int i = 0; i < vec.size(); i++) tm[i] = vec[i]/scalar; 
  return tm; 
}

template<class T> 
T dot(const Vcr<T> & u, const Vcr<T> & v) {
  int size = u.size();
  if (size != v.size()) error("bad vector sizes in dot product");
  T tm = u[0]*v[0];
  for (int i = 1; i < size; i++) tm += u[i]*v[i];
  return tm;
}

// function for dot product for real numbers
template<class T> T dot(T* a, T* b, int n) {
  T init = 0;
  for (int i = 0; i < n; i++) init += *a++ * *b++;
  return init;
}

// *** definitions for members in class Vcr< complex<T> >.

template<class T>    // constructor
Vcr< complex<T> >::Vcr(int n, const complex<T>* const abd) {
  vr = new complex<T> [lenth =n]; 
  for (int i = 0; i < lenth; i++)  vr[i]= *(abd +i);
}

template<class T>
Vcr< complex<T> >::Vcr(int n, complex<T> abd) { // constructor
  vr = new complex<T> [lenth =n]; 
  for (int i = 0; i < lenth; i++)  vr[i]= abd;
}

template<class T>
Vcr< complex<T> >::Vcr(const Vcr & vec) {  //copy constructor
  vr = new complex<T> [lenth = vec.lenth]; 
  for (int i = 0; i < lenth; i++)  vr[i] = vec[i]; 
}
template<class T>          
T Vcr< complex<T> >::maxnorm() const {                  // maximum norm
  T nm = abs(vr[0]);
  for (int i = 1; i < lenth; i++) 
    nm = max(nm, abs(vr[i]));                           // #include <algorithm>
  return nm;
}

template<class T>          
T Vcr< complex<T> >::twonorm() const {                  // two norm
  T nm = pow(abs(vr[0]), 2);
  for (int i = 1; i < lenth; i++) nm += pow(abs(vr[i]),2);
  return sqrt(nm);
}

template<class T>
void Vcr< complex<T> >::reset(complex<T> a = 0) {
 for (int i = 0; i < lenth; i++)  vr[i] = a;
}

template<class T> Vcr<complex<T> >& 
Vcr<complex<T> >::operator=(const Vcr<complex<T> >& v) {
  if (this != &v) {
    if (lenth != v.lenth) error("bad vector sizes in assigment");
    for (int i = 0; i < lenth; i++) vr[i] = v[i];
  }
  return *this;
}

template<class T> Vcr<complex<T> >& 
Vcr<complex<T> >::operator+=(const Vcr<complex<T> >& v) {
  if (lenth != v.lenth) error("bad vector sizes in addition");
  for (int i = 0; i < lenth; i++) vr[i] += v[i];
  return *this;
}

template<class T> Vcr<complex<T> >& 
Vcr<complex<T> >::operator-=(const Vcr<complex<T> >& v) {
  if (lenth != v.lenth) error("bad vector sizes in subtract");
  for (int i = 0; i < lenth; i++) vr[i] -= v[i];
  return *this;
}

template<class T> 
Vcr<complex<T> > operator+(const Vcr<complex<T> >& v) {
  return v;                 // usage: u = + v;
}

template<class T>           // usage: vec1 = - vec2;
Vcr<complex<T> > operator-(const Vcr<complex<T> >& v) {
  Vcr< complex<T> > minus = v;
  for (int i = 0; i < v.size(); i++) minus[i] = - v[i];
  return minus;
}

template<class T> Vcr<complex<T> >       // w=u+v
operator+(const Vcr<complex<T> >& u, const Vcr<complex<T> >& v) {
  if (u.size()!= v.size()) error("bad vector sizes in addition");
  Vcr< complex<T> > sum = u;
  sum += v;
  return sum;
}

template<class T> Vcr<complex<T> >       // w=u-v
operator-(const Vcr<complex<T> >& u, const Vcr<complex<T> >& v) {
  if (u.size() != v.size()) error("bad vector sizes in addition");
  Vcr< complex<T> > minus = u;
  minus -= v;
  return minus;
}

template<class T> Vcr<complex<T> >       // w=a*v
operator*(complex<T> scalar, const Vcr<complex<T> >& v) {
  Vcr< complex<T> > tm(v.size());
  for (int i = 0; i < v.size(); i++) tm[i] = scalar*v[i];
  return tm;
}

template<class T> Vcr<complex<T> >       // w=a*v
operator*(const Vcr<complex<T> >& v, complex<T> scalar) {
  Vcr< complex<T> > tm(v.size());
  for (int i = 0; i < v.size(); i++) tm[i] = scalar*v[i];
  return tm;
}

template<class T> Vcr<complex<T> >       // w=a*v
operator/(const Vcr<complex<T> >& v, complex<T> scalar) {
  Vcr< complex<T> > tm(v.size());
  for (int i = 0; i < v.size(); i++) tm[i] = v[i]/scalar;
  return tm;
}


template<class T> Vcr<complex<T> >       // w=a*v
operator*(const Vcr<complex<T> >& u, const Vcr<complex<T> >& v) {
  int size = u.size();
  if (size != v.size()) error("bad vector sizes in multiply");
  Vcr< complex<T> > tm(size);
  for (int i = 0; i < size; i++) tm[i] = u[i]*v[i];
  return tm;
}


template<class T>                                       // dot product
complex<T> dot(const Vcr<complex<T> >& u, const Vcr<complex<T> >& v) {
  int size = u.size();
  if (size != v.size()) cout<< "bad vector sizes";
  complex<T> tm = u[0]*conj(v[0]);       // conjugation is needed
  for (int i = 1; i < size; i++) tm += u[i]*conj(v[i]);
  return tm;
}

template<class T>      
ostream& operator<<(ostream& s, const Vcr< complex<T> >&v) {
  for (int i =0; i < v.size(); i++) {
    s << v[i] << "  ";
    if (i%10 == 9) s << '\n';
  }
  return s;
}

                           // output operator

// function for dot product for complex numbers
template<class T> 
complex<T> dot(complex<T>* a, complex<T>* b, int k) {
  complex<T> init = 0;
  for (int i = 0; i < k; i++) init += *a++ * conj(*b++);
  return init;
}


// *** definitions of iterative solvers for AbsMtx

template<class T>
int AbsMtx<T>::CG(Vcr<T>& x, const Vcr<T>& b, double& eps,
                  int& iter, int pcn) {
    // Preconditioned conjugate gradient method for Ax = b
    // It returns 0 for sucessful return and 1 if breakdown
    // A: Hermitian positive definite coefficient matrix
    // x:  on entry:  initial guess;
    //     on return: approximate solution
    // b:  right side vector 
    // eps: on entry:  stopping criterion, epsilon
    //      on return: absolute residual in two-norm for 
    //                 approximate solution
    // iter: on entry:  max number of iterations allowed
    //       on return: actual number of iterations taken
    // pcn: = 0 if no preconditioner, = 1 if diag precond
    //      = 2 if SSOR preconditioner

  if (nrows != b.size()) 
    error("matrix and vector sizes do not match");
  const int maxiter = iter;
  Vcr<T> r = b - (*this)*x;        // initial residual
  Vcr<T> z = preconding(r, pcn);   // solve precondit system
  Vcr<T> p = z;                    // p: search direction
  T zr = dot(z,r);                 // inner prod of z and r
  const double stp = eps*b.twonorm();  // stopping criterion

  if (r.maxnorm() == 0.0) {   // if initial guess is true
    eps = 0.0;                // solution, return. Otherwise
    iter = 0;                 // division by zero would occur
    return 0; 
  }

  for (iter = 0; iter < maxiter; iter++) {     // main loop
    Vcr<T> mp = (*this)*p;    // one matrix-vector multiply
    T pap = dot(mp,p);        // one inner product
    T alpha = zr/pap;         // pap is 0 only when r is 0
    x += alpha*p;             // update the iterative soln
    r -= alpha*mp;            // update residual
    if (r.twonorm() <= stp) break;  // stop if converged
    z = preconding(r,pcn);   // preconditioning
    T zrold = zr;
    zr = dot(z,r);            // the other inner product
    T beta= zr/zrold;         // zrold = 0 only when r = 0
    p = z + beta*p;           // update search direction
  } 

  eps = r.twonorm();
  if (iter == maxiter) return 1;
  else return 0;
} // end CG()


template<class T>
int AbsMtx<T>::GMRES(Vcr<T>& x, const Vcr<T>& b, double& eps,
                     int& iter, int pcn, int m) {  
    // Restarted GMRES for solving Ax = b, called GMRES(m) 
    // Return 0 if sucessful and 1 if breakdown
    // It outputs the approximate solution, residual, and 
    // number of iterations taken. 
    //
    // A: any nonsingular matrix, may not be symmetric
    // x: on entry:  initial guess
    //    on return: approximate solution
    // b: right side vector
    // eps: on input, stopping criterion,
    //      on return, absolute two-norm residual
    // iter: on input, max number of iterations allowed 
    //       on return, actual number of iterations taken.
    // pcn: = 0 if no preconditioner, = 1 if diag precond, 
    //      = 2 if SSOR precond
    // m: number of iterations for restarting GMRES

  const int maxiter = iter;
  const double stp = (preconding(b, pcn)).twonorm() * eps;
  Vcr<T> r = preconding(b - (*this) * x, pcn);
  T beta = r.twonorm();             // T may be complex<C>
  bool conv = false;
  if (m > nrows) 
    error("In GMRES(m), m is bigger than number of rows");
  if (abs(beta) <= stp) {    // return if initial guess 
    eps = abs(beta);         // is true solution
    iter = 0;
    return 0;
  }

  // orthonormal basis for Krylov space,
  // v[i] is a pointer to ith basis vector
  Vcr<T>** v = new Vcr<T>* [m+1];    
  for (int i = 0; i <= m; i++) 
    v[i] = new Vcr<T>(nrows);  // ith orthonormal basis

  // Hessenburg matrix h, (m+1) by m;
  // h[i] stores ith column of h that has i+2 entries.
  // Only non-zero part of h is stored
  T** h = new T* [m];     
  for (int i = 0; i < m; i++) h[i] = new T [i+2];              

  iter = 1;
  while (iter  <= maxiter) {  // iterations for gmres(m)
    *v[0] = r / beta;  
    Vcr<T> g(m+1);    // vector in least squares problem
    g[0] = beta;

    Vcr<T> cs(m), sn(m);      // Givens rotations
    int k;  
    for (k = 0; k < m && iter <= maxiter; k++, iter++) {

      // orthogonalization
      Vcr<T> w = preconding((*this) * (*v[k]), pcn);
      T nmw = w.twonorm();
      for (int i = 0; i <= k; i++) {     
        h[k][i] = dot(w, *v[i]);
        w -= h[k][i] * (*v[i]);
      }
      h[k][k+1] = w.twonorm();
      // if h[k][k+1] is small, do reorthogonalization
      if (nmw + 1.0e-4*h[k][k+1] == nmw) {
        for (int i = 0; i <= k; i++) { 
          T hri = dot(w, *v[i]);   // re-orthogonalization
          h[k][i] += hri;
          w -= hri * (*v[i]);
        }
        h[k][k+1] = w.twonorm();
      }
      if (h[k][k+1] == 0) 
        error("unexpected zero-divisor in GMRES()");
      *v[k+1] = w / h[k][k+1]; 

      // apply Givens G_0, G_1, ...,G_{k-1} to column k of h
      for (int i = 0; i < k; i++) { 
//      T tmp = conj(cs[i])*h[k][i] + conj(sn[i])*h[k][i+1];
        T tv[2] = { cs[i], sn[i] };
        T tmp = dot(&h[k][i], tv, 2);
        h[k][i+1] = - sn[i]*h[k][i] + cs[i]*h[k][i+1];
        h[k][i] = tmp;
      }

      // generate Givens rotation G_k from kth column of h
      if (h[k][k+1] == 0) {
        cs[k] = 1;
        sn[k] = 0;
      } else {
        T tpm = sqrt(dot(&h[k][k], &h[k][k], 2));
        cs[k] = h[k][k]/tpm;
        sn[k] = h[k][k+1]/tpm;
      }

      // apply Givens rotation G_k to kth column of h and to g
      T tv[2] = { cs[k], sn[k] };
      h[k][k] = dot(&h[k][k], tv, 2);
      T tmp = dot(&g[k], tv, 2);
      g[k+1] = - sn[k]*g[k] + cs[k]*g[k+1];
      g[k] = tmp;

      if (abs(g[k+1]) <= stp) { // stop if residual small
        k++;
        break;
      }
    }

    // back solve the upper triangular system
    for (int i = k-1; i >= 0; i--){
      for (int j = i+1; j < k; j++) g[i] -= h[j][i]*g[j];
      g[i] /= h[i][i];
    }