matvec.h

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

}

template<class T>                 // solve Pz = r, return z. (P: preconditioner)
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 preconditioning 
    const T omega = 1;            // SSOR parameter for preconditioning
    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 in FullMtx::preconding() not implemented");
  } 
} // end FullMtx::preconding()


/*
template<class T>          // LU decomposition with partial (column) pivoting
Vcr<int> FullMtx<T>::LUdecompnCP() {
  Vcr<int> pvt(nrows);                          // store pivoting info
  for (int k = 0; k < nrows; k++) pvt[k] = k;

  for (int k = 0; k < nrows - 1; k++) {
    // find the pivot in column k below main diagonal
    int pc  = k; 
    double aet = abs(mx[k][k]);
    for (int i = k + 1; i < nrows; i++) {
      if (abs(mx[i][k]) > aet) {
        aet = abs(mx[i][k]); 
        pc = i;
      }
    }
    if (aet == 0) error("pivot is zero in FullMtx::LUdecompn()"); 
    swap(pvt[k], pvt[pc]);

    // interchange pivot row and k-th row
    if (pc != k)
      for (int j = k; j < ncols; j++) swap(mx[pc][j], mx[k][j]);

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

template<class T>     // Gauss Elim without pivoting
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 
  int nrowsmone = nrows - 1;
  for (int k = 0; k < nrowsmone; 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 = nrowsmone; 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 row interchanges
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;

  int nrowsmone = nrows - 1;
  for (int k = 0; k < nrowsmone; 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 = nrowsmone; 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 row interchanges, which normally
// should be less efficient
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;

  const int nrowsmone = nrows - 1;
  for (int k = 0; k < nrowsmone; 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 = nrowsmone; 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;

  const int nrowsmone = nrows - 1;
  for (int k = 0; k < nrowsmone; 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 = nrowsmone; 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, T* 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.vr[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>& vec) {             // usage: vec1 = - vec2;
  Vcr<T> minus = vec;
  for (int i = 0; i < vec.lenth; i++) minus[i] = - vec[i];
  return minus;
}

template<class T> 
Vcr<T> operator+(const Vcr<T>& v1, const Vcr<T>& v2) {         // v=v1+v2
  if (v1.lenth != v2.lenth ) 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.lenth != v2.lenth ) 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.lenth; i++ ) {
    s << vec[i] << "  ";
    if (i%10 == 9) s << "\n";
  }
  return s;
}

template<class T> T Vcr<T>::twonorm() const {
  T nm = abs(vr[0]);
  for (int i = 1; i < lenth; i++) {
    T vi = abs(vr[i]);
    if (nm < 100) nm = sqrt(nm*nm + vi*vi);
    else {                  // to avoid overflow for fn "sqrt" when nm is large
      T tm = vi/nm;
      nm *= sqrt(1.0 + tm*tm);
    }
  } 
  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;