matvec.h

C++ source code for book-C++ and Object Oriented Numeric computing for scientists and engineers

    }
  } 
  return nm;
}

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.lenth);
  for (int i = 0; i < vec.lenth; i++) tm[i] = scalar*vec[i]; 
  return tm;
}

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

template<class T> 
Vcr<T> operator*(const Vcr<T> & vec1, const Vcr<T> & vec2) {
  if (vec1.lenth != vec2.lenth ) error("bad vector sizes in vector multiply");
  Vcr<T> tm(vec1.lenth);
  for (int i = 0; i < vec1.lenth; 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.lenth);
  for (int i = 0; i < vec.lenth; i++) tm[i] = vec[i]/scalar; 
  return tm; 
}

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

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 of iterative solvers for AbsMtx

template<class T>
int AbsMtx<T>::CG(Vcr<T>& x, const Vcr<T>& b, double& eps,  
	                       int& iter, int prec) {
// Conjugate gradient method. 
// x: on entry, initial guess; on retrun: approximate solution 
// b: right hand side vector as in A x = b;
// eps:  on entry, tolerance; on retrun: absolute residual in Euclid norm 
// iter: on entry, max number of iterations allowed; 
//       on return, actual number of iterations used 
// prec= 0 if no preconditioning, 1 if diag prec, 2 if SSOR prec

  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,prec);              // solve the precond 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 intial guess is true soln, 
    eps = 0.0;                                // return. Otherwise division by 
    iter = 0;                                 // zero would occur.
    return 0; 
  }

  for (iter = 0; iter < maxiter; iter++) {     // main loop of CG method
    Vcr<T> mp = (*this)*p;                     // one matrix-vector multiply
    T pap = dot(mp,p);                         // one of two inner products
    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 convergence achieved
    z = preconding(r,prec);                    // preconditioning
    T zrold = zr;
    zr = dot(z,r);                             // another of two inner products
    T beta= zr/zrold;                          // zrold = 0 only when r is 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.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.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];
    }

    // update the solution x = x0 + sum v[i]*y[i]
    for (int i = 0; i < k; i++) x += (*v[i])*g[i];

    // calculate residual and check for convergence
    r = preconding(b - (*this) * x, pcn);
    beta = r.twonorm();
    if (abs(beta) <= stp ) {    // stop if residual small
      conv = true;
      break;
    }
  }

  eps =  (b - (*this) * x).twonorm();  // get final residual

  // deallocate space for v and h
  for (int i = 0; i <= m; i++ ) delete v[i];
  delete[] v;
  for (int i = 0; i < m; i++) delete[] h[i];
  delete[] h;

  if (conv) return 0;
  else return 1;
}   // end of gmres(m)


// *** definitions for members of class BandMtx

template<class T>
BandMtx<T>::BandMtx(int n, int p, int r, T** t) {
  nrows = n;                                 // number of rows
  bwleft = p;                                // left bandwidth
  bwrit = r;                                 // right bandwidth
  bdmx = new T* [nrows]; 
  for (int i = 0; i < nrows; i++) {
   bdmx[i] = new T [bwleft + bwrit + 1];
   bdmx[i] += bwleft;                              
  }
  for (int i = 0; i < nrows; i++) 
    for (int j = - bwleft; j <= bwrit; j++) bdmx[i][j] = t[i][j];
}

template<class T>
BandMtx<T>::BandMtx(int n, int p, int r, T t = 0) {
  nrows = n;
  bwleft = p;
  bwrit = r;
  bdmx = new T* [nrows]; 
  for (int i = 0; i < nrows; i++) {
   bdmx[i] = new T [bwleft + bwrit + 1];
   bdmx[i] += bwleft;                              
  }
  for (int i = 0; i < nrows; i++) 
    for (int j = - bwleft; j <= bwrit; j++) bdmx[i][j] = t;
}

template<class T>
BandMtx<T>::BandMtx(const BandMtx & bd) {
  nrows = bd.nrows;
  bwleft = bd.bwleft;
  bwrit = bd.bwrit;
  bdmx = new T* [nrows]; 
  for (int i = 0; i < nrows; i++) {
   bdmx[i] = new T [bwleft + bwrit + 1];
   bdmx[i] += bwleft;                              
  }
  for (int i = 0; i < nrows; i++) 
    for (int j = - bwleft; j <= bwrit; j++) bdmx[i][j] = bd[i][j];
}


template<class T>
BandMtx<T>::BandMtx(int n, int p, int r, const FullMtx<T>& fm) {
  nrows = n;
  bwleft = p;
  bwrit = r;
  bdmx = new T* [nrows]; 
  for (int i = 0; i < nrows; i++) {
   bdmx[i] = new T [bwleft + bwrit + 1];
   bdmx[i] += bwleft;                              
  }
  for (int i = 0; i< nrows; i++) {
    for (int j = - bwleft; j <= bwrit; j++) 
    bdmx[i][j] = 0;
  }
  for (int i = 0; i< nrows; i++) {                   // fm may be non-symmetric
    int ip = max(i-bwleft, 0); 
    int ir = min(i+bwrit, nrows -1); 
    for (int j = ip; j <= ir; j++) bdmx[i][j-i] = fm[i][j];
  }
}

template<class T>
BandMtx<T>& BandMtx<T>::operator=(const BandMtx & bd) {
  if (nrows != bd.nrows || bwleft != bd.bwleft || bwrit != bd.bwrit) 
    error("bad matrix size in BandMtx assignment");
  for (int j = 0; j < nrows; j++) 
    for (int k = -bwleft; k <= bwrit; k++) bdmx[j][k] = bd[j][k];
  return *this;
}

template<class T>                                     // matrix vector multiply 
Vcr<T> BandMtx<T>::operator*(const Vcr<T> & vec) const {    
  if (nrows != vec.size()) error("matrix and vector sizes do not match");
  Vcr<T> tm(nrows);
  for (int i = 0; i < nrows; i++) {
    int ip = max(i-bwleft, 0); 
    int ir = min(i+bwrit, nrows - 1);
    for (int j = ip; j <= ir; j++) tm[i] += bdmx[i][j-i]*vec[j]; 
  }
  return tm;
}

template<class T>                      // solve Pz = r, return z
Vcr<T> BandMtx<T>::preconding(const Vcr<T> & r, int precn) const {
  if (precn == 0) {                    // no preconditioning
    return r;
  } else if (precn == 1) {             // diagonal preconditioning
    Vcr<T> z(nrows);
    for (int i = 0; i < nrows; i++) z[i] = r[i]/bdmx[i][0];
    return z;
  } else if (precn == 2) {             // symmetric SOR preconditioning
    const T omega = 1;                 // SSOR parameter for preconditioning
    Vcr<T> z(nrows);
    for (int i = 0; i < nrows; i++) {
      T sum = 0;
      int ip = max(i-bwleft, 0);
      for (int j = ip; j < i; j++)  sum += omega*bdmx[i][j-i]*z[j];
      z[i] = (r[i] - sum)/bdmx[i][0];
    }
    for (int i = nrows -1; i >= 0; i--) {
      T sum = 0;
      int ir = min(i+bwrit, nrows -1); 
      for (int j = i+1; j <= ir; j++) sum += omega*bdmx[i][j-i]*z[j];
      z[i] -= sum/bdmx[i][0];
    }
    return z;
  } else {
    error("specified preconditioner in BandMtx::preconding() not implemented");
  }
}

// BandGauss
template<class T>
void BandMtx<T>::GaussElim(Vcr<T>& bb) const {
                           // Gauss elim without pivoting
  if (nrows != bb.size()) 
    error("matrix and vector sizes do not match");
  BandMtx<T> tmx = *this;

  // banded LU decomposition without pivoting
  const int nrowsmone = nrows - 1;
  for (int k = 0; k < nrowsmone; k++) {
    if (tmx[k][0] == 0.0) 
      error("pivot is zero in BandMtx::GaussElim()");
    int kbf = min(nrowsmone - k, bwleft);
    int kbr = min(nrowsmone - k, bwrit);
    for (int i = 1; i <= kbf; i++) {
      int kpi = k+i;
      if (tmx[kpi][-i] != 0.0) {   
        T dmul = tmx[kpi][-i]/tmx[k][0];
        tmx[kpi][-i] = dmul;  // tmx[k+i][-i] = a[k+i][k]
        for (int j = 1; j <= kbr; j++) 
          tmx[kpi][j-i] -= dmul*tmx[k][j];   
      }                       // a[k+i][k+j] = tmx[k+i][j-i]
    }
  }

  // Forward substitution
  for (int k = 1; k < nrows; k++) {
    int kf = min(k, bwleft);
    for (int j = 1; j <= kf; j++) 
      bb[k] -= tmx[k][-j]*bb[k-j];
  }                           // tmx[k][-j] = a[k][k-j]

  // Backward substitution
  for (int k = nrowsmone; k >= 0; k--) {