mv.h

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

    // 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) {
  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;                      // fm may be non-symmetric
  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++) {   
    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 preconditioner
    const T omega = 1;       // SSOR parameter for precondit
    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 not implemented");
  }
}

// Gauss elim without pivoting
template<class T>
void BandMtx<T>::GaussElim(Vcr<T>& bb) const {
  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) 
      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) {   
        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--) {
    int kbr = min(nrowsmone - k, bwrit);
    for ( int j = 1; j <= kbr; j++) 
      bb[k] -= tmx[k][j]*bb[k+j];
    bb[k] /= tmx[k][0];
  }
} // end GaussElim()

// Band Gauss elimination with partial pivoting
template<class T>
void BandMtx<T>::GaussElimPP(Vcr<T>& bb) const {
  if (nrows != bb.size()) 
    error("matrix-vector sizes not match");

  BandMtx<T> tx(nrows, bwleft,
                       min(nrows-1, bwleft+bwrit));
  for (int i = 0; i < nrows; i++) 
    for (int j = - bwleft; j <= bwrit; j++) 
      tx[i][j] = bdmx[i][j];

  Vcr<int> pvt(nrows);   // store pivoting info

  // LU decomposition with column partial pivoting
  const int nrowsmone = tx.nrows - 1;
  for (int k = 0; k < nrowsmone; k++) {
    int kbrow = min(nrowsmone - k, tx.bwleft); 
    int kbcol = min(nrowsmone - k, tx.bwrit); 

    // find the pivot in the k-th column
    int pc = k;
    double aet = abs(tx[k][0]);
    for (int i = 1; i <= kbrow; i++) {
      if (abs(tx[k+i][-i]) > aet) {
        aet = abs(tx[k+i][-i]);
        pc = k + i;
      }
    }
    if (!aet) error("pivot is zero in GaussElimPP");
    pvt[k] = pc;

    // interchange pivot row and k-th row in U, not in L
    if (pc != k) {
      for (int j = 0; j <= kbcol; j++) 
        swap(tx[pc][k+j-pc], tx[k][j]);
    }

    // now eliminate column entries 
    for (int i = 1; i <= kbrow; i++) { 
      int kpi = k + i;
      if (tx[kpi][-i] != 0) {
        T dmul = tx[kpi][-i]/tx[k][0];
        tx[kpi][-i] = dmul;
        for (int j = 1; j <= kbcol; j++) 
          tx[kpi][j-i] -= dmul*tx[k][j];
      } 
    }
  }

  // Forward substitution LY = b
  for (int k = 0; k < nrowsmone; k++) {
    int kbrow = min(nrowsmone - k, tx.bwleft); 
    int pvtk = pvt[k];
    T sb = bb[pvtk];
    if (k != pvtk) swap(bb[k], bb[pvtk]);
    for (int j = 1; j <= kbrow; j++) 
      bb[k+j] -= tx[k+j][-j]*sb;
  }

  // Backward substitution U x = y
  for (int k = nrowsmone; k>= 0; k--) {
    int kb = min(nrowsmone -k, tx.bwrit); 
    for (int j = 1; j <= kb; j++) bb[k] -= tx[k][j]*bb[k+j];
    bb[k] /= tx[k][0];
  }
} // end GaussElimPP()

// *** definitions for members of class SparseMtx

template<class T>
SparseMtx<T>::SparseMtx(int n,int m,T* et,int* cn,int* da) {
  nrows = n;
  lenth = m;
  sra = new T [lenth];
  clm = new int [lenth];
  fnz = new int [nrows +1];

  for (int i =0; i< lenth; i++) {
    sra[i] = et[i];
    clm[i] = cn[i];
  }
  for (int i = 0; i <= nrows; i++)  fnz[i] = da[i];
}

template<class T>
SparseMtx<T>::SparseMtx(int n, int m) {
  nrows = n;
  lenth = m;
  sra = new T [lenth];
  clm = new int [lenth];
  fnz = new int [nrows +1];

  for (int i =0; i< lenth; i++)  {
    sra[i] = 0;
    clm[i] = 0;
  }
  for (int i =0; i <= nrows; i++)  fnz[i] = 0;
}

template<class T>           // copy constructor
SparseMtx<T>::SparseMtx(const SparseMtx& mat) {
  nrows = mat.nrows;
  lenth = mat.lenth;
  sra = new T [lenth]; 
  clm = new int [lenth]; 
  fnz = new int [nrows +1]; 
  for (int i = 0; i < lenth; i++) {
    sra[i] = mat[i]; 
    clm[i] = mat.clm[i]; 
  }
  for (int i =  0; i <= nrows; i++) fnz[i] = mat.fnz[i];
}

template<class T>
SparseMtx<T>& SparseMtx<T>::operator=(const SparseMtx& ssm) {
  if(nrows != ssm.nrows || lenth != ssm.lenth) 
    error("bad matrix sizes in SparseMtx::operator=()");
  for (int i = 0; i < lenth; i++) {
    sra[i]  = ssm[i];
    clm[i]  = ssm.clm[i];
  }
  for (int i = 0; i <= nrows; i++) fnz[i] = ssm.fnz[i]; 
  return *this;
}

template<class T>
Vcr<T> SparseMtx<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++)  
    for (int j = fnz[i]; j < fnz[i +1]; j++) 
      tm[i] += sra[j]*vec[clm[j]];
  return tm;
}

template<class T>    // solve preconditioned Pz = r, return z
Vcr<T> SparseMtx<T>::preconding(const Vcr<T>& r, 
                                int precn) const {
  if (precn == 0) {             // no preconditioning
    return r;
  } else if (precn == 1) {      // diagonal preconditioning
    Vcr<T> diag(nrows);         // find diagonal entries
    for (int i = 0; i < nrows; i++) { 
      for (int j = fnz[i]; j < fnz[i+1]; j++) 
        diag[i] += sra[j]*sra[j];  // geom. average, stable
      diag[i] = sqrt(diag[i]); 
    }

    Vcr<T> z(nrows);
    for (int i = 0; i < nrows; i++) z[i] = r[i]/diag[i];
    return z;
  } else if (precn == 2) {      // SSOR preconditioning
    const T omega = 1;          // SOR parameter

    Vcr<T> diag(nrows);         // find diagonal entries
    for (int i = 0; i < nrows; i++) {        
      for (int j = fnz[i]; j < fnz[i+1]; j++) {
        if (clm[j] == i) diag[i] += sra[j]*sra[j];   
        else diag[i] += omega*omega*sra[j]*sra[j];  
      }          // averaging is used, it is more stable
      diag[i] = sqrt(diag[i]);
    }

    Vcr<T> z(nrows);
    for (int i = 0; i < nrows; i++) {  
      T sum = 0;        // solve lower triangular system
      for (int j = fnz[i]; j < fnz[i+1]; j++) {
        if (clm[j] < i) sum += omega*sra[j]*z[clm[j]];
      }
      z[i] = (r[i] - sum)/diag[i];
    }
    for (int i = nrows -1; i >= 0; i--) {
      T sum = 0;        // solve upper triangular system
      for (int j = fnz[i]; j < fnz[i+1]; j++) {
        if (clm[j] > i) sum += omega*sra[j]*z[clm[j]];
      }
      z[i] -= sum/diag[i];
    }
    return z;
  } else {
    error("specified preconditioner not implemented");
  } 
} // end SparseMtx::preconding()

#endif

//****** end mv.h  ***************/