mv.h

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

/******************** mv.h ***************************
 *
 *                        written by Daoqi Yang
 *                        last modified September 28, 1998.
 **********************************************************/

#ifndef MATVEC_H
#define MATVEC_H

#include <cstdlib>
#include <cmath>
#include <iostream.h> 
#include <complex>
#include <algorithm>

//using namespace std;

template<class T> class Vcr;      // forward declaration
template<class T> class Vcr< complex<T> >;      // forward declaration

template<class T> class AbsMtx {  // base matrix
protected:
  int nrows;                      // number of rows
  virtual Vcr<T> preconding(const Vcr<T>& r,int i=0) const=0; 
    // Solve P z = r with preconditiner P and vector r. 
    // It returns z. Here i = 0 if no preconditioning, 
    // i = 1 if diag prec, i = 2 if SSOR  prec
public: 
  virtual Vcr<T> operator*(const Vcr<T> &) const = 0;  
                                  // matrix vector multiply
  int 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
  int GMRES(Vcr<T>& x, const Vcr<T>& b, double& eps, 
            int& iter, int pcn, int m);
    // Restarted GMRES for solving Ax = b 
    // Returns 0 for sucessful return and 1 for breakdowns
    // It outputs the approximate solution, residual, and 
    // number of iterations used. 
    //
    // 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
};

template<class T> class FullMtx: public AbsMtx<T> {
private: 
  int ncols;                  // number of columns in matrix
  T** mx;                     // entries of the matrix
  Vcr<T> preconding(const Vcr<T>&, int i = 0) const; 
    // i = 0: no preconditioner, i = 1: diag precondr, 
    // i = 2: SSOR precondr
public: 
  FullMtx(int n,int m,T**);   // n:# of rows, m:# of columns
  FullMtx(int n,int m,T t=0); // all entries are set to t
  FullMtx(const FullMtx &);   // copy constructor
  ~FullMtx(){                 // destructor
    for (int i = 0; i< nrows; i++) delete[]  mx[i];
    delete[] mx;
  }

  FullMtx& operator=(const FullMtx&);  // overload =
  FullMtx& operator+=(const FullMtx&); // overload +=
  FullMtx& operator-=(const FullMtx&); // overload -=

  FullMtx& operator+();                // unary +, m1 = +m2;
  FullMtx& operator-();                // unary -, m1 = -m2;
  FullMtx operator+(const FullMtx&);   // binary +, m =m1+m2
  FullMtx operator-(const FullMtx&);   // binary -, m =m1-m2
  Vcr<T> operator*(const Vcr<T>&) const;  
                                  // matrix-vector multiply
  T* operator[](int i) const { return mx[i]; }  // subscript
  void GaussElim(Vcr<T>&) const;  // Gauss elimination
  void GaussElimPP(Vcr<T>&) const; 
    // Gauss Elim with partial pivoting, no interchanging rows
  void GaussElimPP2(Vcr<T>&) const; 
    // Gauss Elim with partial pivoting, interchange rows
  void FullMtx<T>::GaussElimCP(Vcr<T>& bb) const;
    // Gauss Elim with complete pivoting
};

template<class T>        // sparse matrix
class SparseMtx: public AbsMtx<T> { 
private:         // compressed sparse row format
  int lenth;     // # of nonzero entries of original matrix
  T *sra;        // array for storing the non-zero entries
  int *clm;      // column indexes of entries stored in sra
  int *fnz;      // 1st nonzero entry of each row in sra
  Vcr<T> preconding(const Vcr<T>&, int i = 0) const;
    // i = 0: no preconditioner, i = 1: diag precondr, 
    // i = 2: SSOR precondr
public: 
  SparseMtx(int n, int m, T* t, int* c, int* f);     
    // n: number of rows (and columns) of original matrix
    // m: length of array sra for nonzero entries
    // t: nonzero entries of the original matrix
    // c: colunm indexes of entries stored in sra 
    // f: index in sra of first nonzero entry in each row
  SparseMtx(int n, int m); // initialize all entries to 0 
    // n: number of rows (and columns) 
    // m: number of nonzero entries. 
  SparseMtx(const SparseMtx &);        // copy constructor
  ~SparseMtx(){ delete[] sra; delete[] fnz; delete[] clm; }
	
  SparseMtx& operator=(const SparseMtx &);   // overload =
  Vcr<T> operator*(const Vcr<T>&) const;     
                         // matrix vector multiply
  T& operator[](int i) const { return sra[i]; } // subscript
  int& getfnz(int i) const { return fnz[i]; }   
                         // first nonzero entry of each row
  int& getclm(int i) const { return clm[i]; } 
                         // column index of entries in sra
};

template<class T> 
class BandMtx: public AbsMtx<T> {   // band matrix
private:
  int bwleft;                       // left bandwidth
  int bwrit;                        // right bandwidth
  T** bdmx;                         // entries within band
  Vcr<T> preconding(const Vcr<T> &, int = 0) const; 
public: 
  BandMtx(int n, int p, int r, T** t);           
    // n: # of rows = # of columns, p: left bandwidth, 
    // r: right bandwidth, t: entries within the band 
  BandMtx(int n, int p, int r, T t = 0);    // constructor 
  BandMtx(int n, int p, int r, const FullMtx<T>& m);
  BandMtx(const BandMtx &);            // copy constructor
  ~BandMtx(){ 
    for (int i = 0; i < nrows; i++) 
      delete[] (bdmx[i] -= bwleft);
    delete[] bdmx; 
  }

  BandMtx& operator=(const BandMtx&);  // overload =
  T* operator[](int i) const { return bdmx[i]; } // row i
  Vcr<T> operator*(const Vcr<T>&) const;  
                                // matrix vector multiply
  void GaussElim(Vcr<T>& b) const;   // Gauss elimination
  void GaussElimPP(Vcr<T>& b) const; 
                           // Gauss elim, partial pivoting
};

template<class T> class Vcr {       
  int lenth;      // number of entries 
  T* vr;          // entries of the vector
public: 
  Vcr(int n, const T* const t);             // constructor from t
  explicit Vcr(int n, T t = 0); // all entries =t
  Vcr(const Vcr&);              // copy constructor
  ~Vcr(){ delete[] vr; }        // destructor

  int size() const { return lenth; }   // size of vector
  void reset(T = 0);       // reset all enties to a value
  T& operator[](int i) const { return vr[i]; }  // v[i]=10;
  Vcr& operator=(const Vcr&);   // overload =, v = v2
  Vcr& operator+=(const Vcr&);  // v += v2
  Vcr& operator-=(const Vcr&);  // v -= v2
  T maxnorm() const;            // maximum norm
  T twonorm() const;            // L-2 norm
};

template<class T>          // unary operator, v = + v2
Vcr<T> operator+(const Vcr<T>&);    

template<class T>          // unary operator, v = - v2
Vcr<T> operator-(const Vcr<T>&);    

template<class T>          // binary operator, v = u + w 
Vcr<T> operator+(const Vcr<T>&, const Vcr<T>&);

template<class T>          // binary operator, v = u - w 
Vcr<T> operator-(const Vcr<T>&, const Vcr<T>&);

template<class T>          // dot product
T dot(const Vcr<T>&, const Vcr<T>&);

template<class T>          // vec-scalar multiply 
Vcr<T> operator*(T, const Vcr<T>&);

template<class T>          // vec-scalar multiply 
Vcr<T> operator*(const Vcr<T>&, T);

template<class T>          // component-wise multiply
Vcr<T> operator*(const Vcr<T>&, const Vcr<T>&);  

template<class T>          // vec-scalar division 
Vcr<T> operator/(const Vcr<T>&, T);

template<class T>          // output operator
ostream& operator<<(ostream&, const Vcr<T>&);    


template<class T> T dot(T* a, T* b, int n);

template<class T> class Vcr< complex<T> > {       
  int lenth;             // number of entries 
  complex<T>* vr;        // entries of the vector
public: 
  Vcr(int n, const complex<T>* const t);  // constructor
  explicit Vcr(int n, complex<T> t = 0);  // constructor
  Vcr(const Vcr&);       // copy constructor
  ~Vcr(){ delete[] vr; } // destructor

  int size() const { return lenth; }  // size of vector
  void reset(complex<T> = 0); // reset all entries to a value 
  complex<T>& operator[](int i) const { return vr[i]; } 
                                // subscript, v[i] = 10;
  Vcr& operator=(const Vcr&);   // overload =, v = v2
  Vcr& operator+=(const Vcr&);  // v += v2
  Vcr& operator-=(const Vcr&);  // v -= v2
  T maxnorm() const;            // maximum norm
  T twonorm() const;            // L-2 norm

};

template<class T>           // unary operator, v = + v2
Vcr< complex<T> > 
operator+(const Vcr< complex<T> >&);    

template<class T>          // unary operator, v = - v2
Vcr< complex<T> > 
operator-(const Vcr< complex<T> >&);    

template<class T>          // v= v1 + v2
Vcr< complex<T> > operator+
(const Vcr< complex<T> >&, const Vcr< complex<T> >&);
template<class T>          // v= v1 - v2

Vcr< complex<T> > operator-
(const Vcr< complex<T> >&, const Vcr< complex<T> >&);

template<class T>          // dot product
complex<T> dot(const Vcr< complex<T> >&, const Vcr< complex<T> >&);

template<class T>          // vec-scalar multiply 
Vcr< complex<T> > operator*
(complex<T>, const Vcr< complex<T> >&);

template<class T>          // vec-scalar multiply 
Vcr< complex<T> > operator*
(const Vcr< complex<T> >&, complex<T>); 

template<class T>          // component-wise multiply
Vcr< complex<T> > operator*
(const Vcr< complex<T> >&, const Vcr< complex<T> >&);

template<class T>          // vec-scalar division 
Vcr< complex<T> > operator/
(const Vcr< complex<T> >&, complex<T>);

template<class T>          // output operator
ostream& operator<<
  (ostream&, const Vcr< complex<T> >&);    

template<class T> 
complex<T> dot (complex<T>* a, complex<T>* b, int n);

/*

class Vcr< complex<double> > { // a complete specialization
  int lenth;                   // number of entries in vector
  complex<double>* vr;         // entries of the vector
public: 
  Vcr(int, complex<double> *cp);            // constructor
  explicit Vcr(int, complex<double> c = 0); // constructor
  Vcr(const Vcr&);             // copy constructor
  ~Vcr(){ delete[] vr;}        // destructor

  int size() const { return lenth; } // size of vector
  void reset(complex<double> t = 0); // reset all enties to t
  complex<double> & operator[](int i) const { return vr[i]; }
  Vcr& operator=(const Vcr&);        // overload =
  Vcr& operator+=(const Vcr&);       // v += v2
  Vcr& operator-=(const Vcr&);       // v -= v2
  double maxnorm() const;            // maximum norm
  double twonorm() const;            // L-2 norm

  friend Vcr operator+<complex<double> >(const Vcr&);
                                    // unary +, v = + v2
  friend Vcr operator-<complex<double> >(const Vcr&);
                                         // unary -, v = - v2
  friend Vcr operator+<complex<double> >
    (const Vcr&, const Vcr&);            // v= v1+v2
  friend Vcr operator-<complex<double> >
    (const Vcr&, const Vcr&);            // v= v1-v2
  friend complex<double> dot<complex<double> >
    (const Vcr&, const Vcr&);            // dot product
  friend Vcr operator*<complex<double> >
    (complex<double>, const Vcr&);       // vec-scalar x 
  friend Vcr operator*<complex<double> >
    (const Vcr&, complex<double>);       // vec-scalar x

  friend Vcr operator*<complex<double> >
    (const Vcr&, const Vcr&);            // component-wise x
  friend Vcr operator/<complex<double> >
    (const Vcr&, complex<double>);       // vec-scalar /
  friend ostream& operator<<<complex<double> >
    (ostream&, const Vcr&);    
};
*/

//***** some usefull small functions and parameters *******

//*** report error and exit with return value 1.
template<class T> void error(const T& t) {
//template<class T> void error(const T& t, ostream& s = cout) {
  cout << t << ". program exited." << "\n";
  exit(1);
}

//const double SMALL = numeric_limits<double>::epsilon(); 
const double SMALL = 1.0e-15;

// *** definitions for members of class FullMtx

template<class T> 
FullMtx<T>::FullMtx(int n, int m, T** dbp) {
  nrows = n;
  ncols = m;
  mx = new T* [nrows]; 
  for (int i =  0; i< nrows; i++) {
    mx[i] = new T [ncols];          
    for (int j = 0; j < ncols; j++) mx[i][j] = dbp[i][j];
  }
}

template<class T> 
FullMtx<T>::FullMtx(int n, int m, T a) {
  nrows = n;
  ncols = m;
  mx = new T* [nrows]; 
  for (int i =  0; i< nrows; i++) {
    mx[i] = new T [ncols];
    for (int j = 0; j < ncols; j++) mx[i][j] = a;
  }
}

template<class T> 
FullMtx<T>::FullMtx(const FullMtx& mat) { 
  nrows = mat.nrows;
  ncols = mat.ncols;
  mx = new T* [nrows]; 
  for (int i =  0; i< nrows; i++) {
    mx[i] = new T [ncols];
    for (int j = 0; j < ncols; j++) mx[i][j] = mat[i][j];
  }
}

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