matvec.h

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

/******************** matvec.h ***************************
 * This is a template library header file for manipulating and solving
 * matrices and vectors, whose entries can be int, float, double,
 * long double, and complex<double>, etc. This file contains

    ~ definitions for template classes 
        * AbsMtx<T>:     an abstract matrix, 
        * Vcr<T>:        a vector class, 
        * FullMtx<T>:    matrix in full storage, 
        * SparseMtx<T>:  sparse matrix in compressed sparse row storage, 
        * BandMtx<T>:    banded matrix, 
        * SSMtx<T>:      symmetric sparse matrix, 
    ~ declaration of complete instantiation Vcr< complex<double> >

  In particular, it contains

    ~ basic operations (eg add, multiply) on matrices and vectors
    ~ norms and dot products for vectors of real and complex entries
    ~ preconditioned conjugate gradient method CG()
    ~ preconditioned restarted GMRES(m) for general nonsingular matrices, 
    ~ Gauss elimination without pivoting GaussElim() for full & band matrices
    ~ Gauss elim with partial pivoting GaussElimPP() for full & band matrices
    ~ Gauss elim with complete pivoting GaussElimCP() for full matrices

 Functions CG(), GMRES(), GaussElim(), GaussElimPP(), GaussElimCP() 
 apply to linear systems of real and complex entries. 
 However, GaussElim() and CG() are guaranted to work only for 
 Hermitian positive definite matrices and GaussElim() 
 also for diagonally dominated matrices. 

 *                        written by Daoqi Yang
 *                        last modified September 28, 1999.
 **********************************************************/

#ifndef MATVEC_H
#define MATVEC_H

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

//using namespace std;

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

template<class T> class AbsMtx {               // base matrix, an abstract class
protected:
  int nrows;                                   // number of rows in the matrix 
  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
      // 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 pn);  
      // preconditioned Conjugate gradient method for Ax=b. A: sym pos def
      // 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.
      // pn: =0 if no preconditioner, =1 if diag precr, =2 if SSOR precr
      // it returns 0 for sucessful return and 1 for breakdowns
  int GMRES(Vcr<T>& x, const Vcr<T>& b, double& eps, int& iter, int pn, int m);
      // restarted GMRES for solve Ax=b, where A may be non-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 upon program termination
      // iter: on input, max number of iterations allowed; 
      //       on return, actual number of iterations taken.
      // pn: =0 if no preconditioner is used, =1 if diag precr, =2 if SSOR precr
      // m: number of iterations for restarting GMRES
      // it returns 0 for sucessful return and 1 for breakdowns
};

template<class T> class FullMtx: public AbsMtx<T> {        // full matrix
private: 
  int ncols;                                    // # of columns in the matrix 
  T** mx;                                       // entries of the matrix
  Vcr<T> preconding(const Vcr<T>&, int i = 0) const; 
      // i = 0: no preconditioner, i = 1: diag precr, i = 2: SSOR precr
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;
  }

  // implement + as a member fcn and  implement - as a friend
  FullMtx& operator=(const FullMtx&);           // overload =
  FullMtx& operator+=(const FullMtx&);          // overload +=
  FullMtx& operator-=(const FullMtx&);          // overload -=
  FullMtx& operator+();                         // unary +, 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]; }  // returns i-th row
  friend FullMtx operator-(const FullMtx&);     // unary -, m1 = -m2;
  friend FullMtx operator-(const FullMtx&,const FullMtx&);  // binary -, m=m1-m2
  friend ostream& operator<<(ostream&, const FullMtx&);     // overload << 
  void GaussElim(Vcr<T>&) const;        // Gauss elimination
  void GaussElimPP(Vcr<T>&) const;      // Gauss Elim with partial pivot
  void GaussElimPP2(Vcr<T>&) const;     // Gauss Elim with partial pivot
  void GaussElimCP(Vcr<T>&) const;      // Gauss Elim with complte pivot
};

template<class T>
class SSMtx: public AbsMtx<T> {                 // Symmetric Sparse matrix
private: 
  int lenth;        // # of diag entries and nonzero off-diag entries 
                    // in the lower triangular part of the matrix, 
  T* sra;           // store diag entries and non-zero off-diagonal entries 
  int *dgn;         // position of matrix diagonals as stored in sra 
  int *clm;         // columns in original matrix for entries in sra 
  Vcr<T> preconding(const Vcr<T> &, int = 0) const; 
      // solve preconditioning matrix P z = r with preconditioner P
      // int = 0: no preconditioning, = 1: diag precond, = 2: SSOR precond
public: 
  SSMtx(int n, int m, T* t, int* c, int* i);     // n = nrows, m = lenth
      // t: 1D array storing diag entries and nonzero off-diag entries
      // c: colunms (in the original matrix) of entries of the 1D array t 
      // i: index of diags in the 1D array.
  SSMtx(int n, int m);                           // initialize a zero SSMtx
      // n: number of rows of the Sym Sparse matrix = number of columns
      // m: # of diag entries & nonzero off-diag entries in lower triangular part
  SSMtx(const SSMtx &);                          // copy constructor
  ~SSMtx(){ delete[] sra; delete[] dgn; delete[] clm; }
	
  SSMtx& operator=(const SSMtx & );              // overload =
  T& operator[](int i) const { return sra[i]; }
  int& getdgn(int i) const { return dgn[i]; }
  int& getclm(int i) const { return clm[i]; }
  Vcr<T> operator*(const Vcr<T>&) const;
};

template<class T>                            // sparse matrix
class SparseMtx: public AbsMtx<T> {          // compressed sparse row format
private: 
  int lenth;             // # of nonzero entries of the original matrix 
  T *sra;                // array for storing the non-zero entries 
  int *clm;              // column indexes in matrix of the entries in sra 
  int *fnz;              // position in sra of first nonzero entires of each row
  Vcr<T> preconding(const Vcr<T>&, int i = 0) const; 
      // i = 0: no preconditioner, i = 1: diag precr, i = 2: SSOR precr
public: 
  SparseMtx(int n, int m, T* t, int* c, int* f);     
      // n: number of rows (and columns) of the original matrix
      // m: length of array sra for nonzero entries.
      // t: nonzero entries of the original matrix
      // c: colunm indexes (in the original matrix) of entries in sra 
      // f: index in sra of first nonzero entry in each row
  SparseMtx(int n, int m);                       // intialize all entris to zero
      // 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]; }  // subscripting
  int& getfnz(int i) const { return fnz[i]; }    // 1st nonzero entry ofeach row
  int& getclm(int i) const { return clm[i]; }    // column index
};

template<class T> 
class BandMtx: public AbsMtx<T> {                // banded matrix
private:
  int bwleft;                                    // left bandwidth 
  int bwrit;                                     // right bandwidth 
  T** bdmx;                                      // entries within the 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);         // initialize all entries to t
  BandMtx(int n, int p, int r, const FullMtx<T>& m);    // m: full matrix format
  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]; } // i-th row of bdmx
  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, column pivoting
  void GaussElimPP2(Vcr<T>& b) const;            // Gauss elim, column pivoting
//  void LUdecompn();                            // LU decomposition
//  Vcr<int> LUdecompnCP();                      // LU decompn, column pivoting
};

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

  int size() const { return lenth; }             // return 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
  friend Vcr operator+(const Vcr&);              // unary operator, v = + v2
  friend Vcr operator-(const Vcr&);              // unary operator, v = - v2
  friend Vcr operator+(const Vcr&, const Vcr&);  // v= v1 + v2
  friend Vcr operator-(const Vcr&, const Vcr&);  // v= v1 - v2
  friend T dot(const Vcr&, const Vcr&);          // dot product
  friend Vcr operator*(T, const Vcr&);           // vec-scalar multiply 
  friend Vcr operator*(const Vcr&, T);           // vec-scalar multiply 
  friend Vcr operator*(const Vcr&, const Vcr&);  // component-wise multiply
  friend Vcr operator/(const Vcr&, T);           // vec-scalar division 
  friend ostream& operator<<(ostream&, const Vcr&);    
};

template<class T> T dot(T* a, 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);                 // construct a vector from cp
  Vcr(int, complex<double> c = 0);      // constructor, all entries = c
  Vcr(const Vcr&);                               // copy constructor
  ~Vcr(){ delete[] vr;}                          // destructor

  int size() const { return lenth; }             // return 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&);    
  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+(const Vcr&);              // unary +, v = + v2
  friend Vcr operator-(const Vcr&);              // unary -, v = - v2
  friend Vcr operator+(const Vcr&, const Vcr&);  // v= v1+v2
  friend Vcr operator-(const Vcr&, const Vcr&);  // v= v1-v2
  friend complex<double> dot(const Vcr&, const Vcr&);      // dot product
  friend Vcr operator*(complex<double>, const Vcr&);       // vec-scalar x 
  friend Vcr operator*(const Vcr&, complex<double>);       // vec-scalar x
  friend Vcr operator*(const Vcr&, const Vcr&);            // component-wise x  
  friend Vcr operator/(const Vcr&, complex<double>);       // vec-scalar / 
  friend ostream& operator<<(ostream&, const Vcr&);    
};

complex<double> dot(complex<double>* a, complex<double>* b, int n);

//***** 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];
  }
  return *this;
}

template<class T>                     // usefull for testing small matrices 
ostream& operator<<(ostream& s, const FullMtx<T>& mat) {
  for (int i =0; i< mat.nrows; i++) {
    s << "\n" << i << "-th row:    \n";
    for (int j =0; j< mat.ncols; j++) {
      s << mat.mx[i][j] << "  ";
      if (j%10 ==9) s << "\n";
    }
    s << "\n";
  }
  return s;
}

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+() {                // usage: mat1 = + mat2;
  return *this;
}

template<class T> 
FullMtx<T> operator-(const FullMtx<T>&  mat) {        // usage: mat1 = - mat2;
  FullMtx<T> minus = mat;  
  for (int i = 0; i < mat.nrows; i++) {
    for (int j = 0; j < mat.ncols; j++) minus[i][j] = - mat[i][j];
  }
  return minus;
}

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

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 in FullMtx::operator*()");
  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> 
FullMtx<T> operator-(const FullMtx<T> & mat1, const FullMtx<T>& mat2) {
 if(mat1.nrows !=mat2.nrows || mat1.ncols !=mat2.ncols) 
   error("bad matrix sizes in FullMtx subtraction");
 FullMtx<T> sum = mat1;
 sum -= mat2;
 return sum;