mv.h

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

  }
  return *this;
}

template<class T> 
FullMtx<T> & FullMtx<T>::operator+=(const FullMtx &mat) {
  if (nrows !=mat.nrows || ncols !=mat.ncols) 
    error("bad matrix sizes in FullMtrx::operator+=()");
  for (int i = 0; i < nrows; i++)
    for (int j = 0; j < ncols; j++) mx[i][j] += mat[i][j];
  return *this;
}

template<class T> 
FullMtx<T> & FullMtx<T>::operator-=(const FullMtx &mat) {
  if (nrows !=mat.nrows || ncols !=mat.ncols) 
    error("bad matrix sizes in FullMtrx::operator-=()");
  for (int i = 0; i < nrows; i++)
    for (int j = 0; j < ncols; j++) mx[i][j] -= mat[i][j];
  return *this;
}

template<class T> 
FullMtx<T>& FullMtx<T>::operator+() {  // m1 = + m2
  return *this;
}

template<class T> 
FullMtx<T>& FullMtx<T>::operator-() {  // m1 = - m2
  FullMtx<T> zero(nrows, ncols);
  return zero - *this;
}

template<class T> 
FullMtx<T> FullMtx<T>::operator+(const FullMtx& mat) {
 FullMtx sum = *this;            // m = m1 + m2
 sum += mat;
 return sum;
}

template<class T> 
FullMtx<T> FullMtx<T>::operator-(const FullMtx& mat) {
 FullMtx minus = *this;            // m = m1 - m2
 minus -= mat;
 return minus;
}

template<class T> 
Vcr<T> FullMtx<T>::operator*(const Vcr<T> & vec) const {
  if (ncols != vec.size()) 
    error("matrix and vector sizes do not match");
  Vcr<T> tm(nrows);
  for (int i = 0; i < nrows; i++) 
    for (int j=0; j<ncols; j++) tm[i] += mx[i][j]*vec[j]; 
  return tm;
}

template<class T>   // solve Pz = r, return z
Vcr<T> FullMtx<T>::preconding(const Vcr<T>& r, 
                              int precn) const {
  if (precn == 0) {           // no preconditioner is used
    return r;
  } else if (precn == 1) {    // diagonal preconditioning
    Vcr<T> z(nrows);
    for (int i = 0; i < nrows; i++) z[i] = r[i]/mx[i][i];
    return z;
  } else if (precn == 2) {    // symmetric SOR precondition
    const T omega = 1;        // SSOR parameter for precond
    Vcr<T> z(nrows);
    for (int i = 0; i < nrows; i++) {
      T sum = 0;
      for (int j=0; j<i; j++) sum += omega*mx[i][j]*z[j];
      z[i] = (r[i] - sum)/mx[i][i];
    } 
    for (int i = nrows -1; i >= 0; i--) {
      T sum = 0;
      for (int j = i +1; j < nrows; j++) 
        sum += omega*mx[i][j]*z[j];
      z[i] -= sum/mx[i][i];
    }
    return z;
  } else {
    error("specified preconditioner not implemented");
  } 
} // end FullMtx::preconding()

template<class T> 
void FullMtx<T>::GaussElim(Vcr<T>& bb) const{
  if (nrows != bb.size() || nrows != ncols) 
    error("matrix or vector sizes do not match");
  FullMtx<T> tmpx = *this;
  
  // LU decomposition without pivoting 
  for (int k = 0; k < nrows - 1; k++) {
    if (tmpx[k][k] == 0) 
      error("pivot is zero in FullMtx::GaussElim()");
    for (int i = k + 1; i < nrows; i++) {
      if (tmpx[i][k] != 0) {   // tmpx[i][k] can be complex
        T mult = tmpx[i][k]/tmpx[k][k];
        tmpx[i][k] = mult;
        for (int j = k + 1; j < ncols; j++) 
          tmpx[i][j] -= mult*tmpx[k][j];
      }
    }
  }

  // forwad substitution for L y = b. y still stored in bb
  for (int i = 1; i < nrows; i++) 
    for (int j = 0; j < i; j++) bb[i] -= tmpx[i][j]*bb[j];

  // back substitution for U x = y. x still stored in bb
  for (int i = nrows - 1; i >= 0; i--) {
    for (int j=i+1; j<nrows; j++) bb[i] -= tmpx[i][j]*bb[j];
    bb[i] /= tmpx[i][i];
  }
} // end GaussElim()

// Gauss Elim with partial pivoting
// This version does not do explicit row interchange
template<class T>     
void FullMtx<T>::GaussElimPP(Vcr<T>& bb) const {
  if (nrows != bb.size() || nrows != ncols) 
    error("matrix or vector sizes do not match");
  FullMtx<T> tmpx = *this;

  Vcr<int> pvt(nrows);       // pivoting vector 
  for (int k = 0; k < nrows; k++) pvt[k] = k;

  for (int k = 0; k < nrows - 1; k++) {  // main loop

    // find the pivot in column k in rows pvt[k], 
    // pvt[k+1], ..., pvt[n-1]
    int pc  = k; 
    double aet = abs(tmpx[pvt[k]][k]);
    for (int i = k + 1; i < nrows; i++) {
      if (abs(tmpx[pvt[i]][k]) > aet) {
        aet = abs(tmpx[pvt[i]][k]); 
        pc = i;
      }
    }
    if (!aet) 
      error("pivot is zero in FullMtx::GaussElimPtP()"); 
    if (pc != k) swap(pvt[k], pvt[pc]);
    int pvtk = pvt[k];                  // pivot row
    T pivot = tmpx[pvtk][k];            // pivot

    // now eliminate column entries logically 
    // below tmpx[pvt[k]][k]
    for (int i = k + 1; i < nrows; i++) {
      int pvti = pvt[i];
      if (tmpx[pvti][k] != 0) {
        T mult = tmpx[pvti][k]/pivot;
        tmpx[pvti][k] = mult;
        for (int j = k + 1; j < ncols; j++) 
          tmpx[pvti][j] -= mult*tmpx[pvtk][j];
      }
    }
  }

  // forwad substitution for L y = Pb.
  for (int i = 1; i < nrows; i++)  
    for (int j = 0; j < i; j++) 
      bb[pvt[i]] -= tmpx[pvt[i]][j]*bb[pvt[j]];

  // back substitution for Ux = y
  Vcr<T> xx(nrows);  // xx stores solution in correct order
  for (int i = nrows - 1; i >= 0; i--) {
    for (int j = i+1; j < ncols; j++) 
      bb[pvt[i]] -= tmpx[pvt[i]][j]*xx[j];
    xx[i] = bb[pvt[i]] / tmpx[pvt[i]][i];
  }

  bb = xx;             // put solution in bb
} // end GaussElimPP()


// Gauss Elim with partial (column) pivoting
// This version does explicit row interchange
template<class T>   
void FullMtx<T>::GaussElimPP2(Vcr<T>& bb) const {

// this function interchanges row k and pivot row

  if (nrows != bb.size() || nrows != ncols) 
    error("matrix or vector sizes do not match in FullMtx::GaussElimPtP()");
  FullMtx<T> tmpx = *this;

  Vcr<int> pvt(nrows);                          // pivoting vector 
  for (int k = 0; k < nrows; k++) pvt[k] = k;

  for (int k = 0; k < nrows - 1; k++) {         // main loop

    // find the pivot in column k below tmpx[k][k]
    int pc  = k; 
    double aet = abs(tmpx[k][k]);
    for (int i = k + 1; i < nrows; i++) {
      if (abs(tmpx[i][k]) > aet) {
        aet = abs(tmpx[i][k]); 
        pc = i;
      }
    }
    if (!aet) error("pivot is zero in FullMtx::GaussElimCP()"); 
    if (pc != k) {
      swap(pvt[k], pvt[pc]);
      for (int j = 0; j < ncols; j++) swap(tmpx[pc][j], tmpx[k][j]);
    }

    // now eliminate the column entries below tmpx[k][k]
    for (int i = k + 1; i < nrows; i++) {
      if (tmpx[i][k] != 0) {
        T mult = tmpx[i][k]/tmpx[k][k];
        tmpx[i][k] = mult;
        for (int j = k + 1; j < ncols; j++) 
          tmpx[i][j] -= mult*tmpx[k][j];
      }
    }
  }

  Vcr<T> rhv(nrows);
  for (int k = 0; k < nrows; k++) rhv[k] = bb[pvt[k]];
  bb = rhv;

  // forwad substitution for L y = Pb.
  for (int i = 1; i < nrows; i++)  
    for (int j = 0; j < i; j++) bb[i] -= tmpx[i][j]*bb[j];

  // back substitution  for Ux = y
  for (int i = nrows - 1; i >= 0; i--) {
    for (int j = i+1; j < ncols; j++) bb[i] -= tmpx[i][j]*bb[j];
    bb[i] /= tmpx[i][i];
  }

} // end GaussElimPP2()

template<class T>    // Gauss Elim with complete pivoting
void FullMtx<T>::GaussElimCP(Vcr<T>& bb) const {
  if (nrows != bb.size() || nrows != ncols) 
    error("matrix vector sizes no match in GaussElimCP()");
  FullMtx<T> tmpx = *this;

  Vcr<int> px(nrows);             // row pivoting vector 
  Vcr<int> qy(nrows);             // column pivoting vector 
  for (int k = 0; k < nrows; k++) px[k] = qy[k] = k;

  for (int k = 0; k < nrows - 1; k++) {  // main loop

    // find pivot entry in columns qy[k], qy[k+1], ..., 
    // qy[n-1] and in rows px[k], px[k+1], ..., px[n-1]
    int pct = k, qdt = k; 
    double aet = 0;
    for (int i = k; i < nrows; i++) {
      for (int j = k; j < nrows; j++) {
        double tmp = abs(tmpx[px[i]][qy[j]]); 
        if (tmp > aet) { aet = tmp; pct = i; qdt = j; }
      }
    }
    if (!aet) error("pivot is zero in GaussElimCP()"); 
    swap(px[k], px[pct]);    // swap px[k] and px[pct]
    swap(qy[k], qy[qdt]);

    // eliminate column entries logically below and right
    // to the entry mx[px[k]][qy[k]]
    for (int i = k + 1; i < nrows; i++) {
      if (tmpx[px[i]][qy[k]] != 0) {
        T mult = tmpx[px[i]][qy[k]]/tmpx[px[k]][qy[k]]; 
        tmpx[px[i]][qy[k]] = mult;
        for (int j = k + 1; j < nrows; j++) 
          tmpx[px[i]][qy[j]] -= mult*tmpx[px[k]][qy[j]];
      }
    }
  }

  // forwad substitution for L y = Pb. Store y in b
  for (int i = 1; i < nrows; i++)  
    for (int j = 0; j < i; j++) 
      bb[px[i]] -= tmpx[px[i]][qy[j]]*bb[px[j]];

  // back substitution for Uz = y and x = Q z
  Vcr<T> xx(nrows);    // xx stores solution      
  for (int i = nrows - 1; i >= 0; i--) {
    for (int j = i+1; j < nrows; j++) 
      bb[px[i]] -= tmpx[px[i]][qy[j]]*xx[qy[j]];
    xx[qy[i]] = bb[px[i]] / tmpx[px[i]][qy[i]];
  }

  bb = xx;
} // end GaussElimCP()

// *** definitions for members in class Vcr.

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

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

template<class T> 
Vcr<T>::Vcr(const Vcr & vec) {
  vr = new T [lenth = vec.lenth]; 
  for (int i = 0; i < lenth; i++)  vr[i] = vec[i]; 
}

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

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

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

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

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

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

template<class T> 
Vcr<T> operator-(const Vcr<T>& v1, const Vcr<T> & v2) {         // v=v1-v2
  if (v1.size()!= v2.size()) error("bad vector sizes in vector subtraction");
  Vcr<T> sum = v1;           // It would cause problem without copy constructor
  sum -= v2;
  return sum;
}

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

template<class T> T Vcr<T>::twonorm() const {
  T nm = vr[0]*vr[0];
  for (int i = 1; i < lenth; i++) nm += vr[i]*vr[i];
  return sqrt(nm);