求解的方法：有限差分法（FDM）
求解第二类边界条件的Poisson方程，其方程如下：
上述的k为一常数。
为了保证解的唯一性，在此令k为一个小量或者微量。
上述方程的精确解为
利用c/C++来实现，暂时还没有用到线性方程组求解器，只是测试系数矩阵和右端向量是否正确。
代码在VS2008平台上来操作。
包含以下几个文件：GMRES.h matvec.cpp matvec.h solpoisson.h testposmain.cpp 共5个。
GMRES.h
View Code
#include "matvec.h" void gmres(SparseMtx<double> a,Vcr<double>& x, const Vcr<double>& b, double& eps, int& iter, int pn, int m) { int ret = a.GMRES(x, b,eps,iter,pn,m); // return ret; }
matvec.cpp
View Code
/******************** matvec.cc *************************** * This file contains definition of a complete specialization * Vcr< complex< double> >. It contains ~ add, subtract, multiply complex vectors ~ norm and dot product of complex vectors * written by Daoqi Yang * last modified September 28, 1998. **********************************************************/ #include "matvec.h" // *** definitions for members in class Vcr< complex<double> >. Vcr< complex<double> >::Vcr(int n, complex<double>* abd) { vr = new complex<double> [lenth =n]; for (int i = 0; i < lenth; i++) vr[i]= *(abd +i); } Vcr<complex<double> >::Vcr(int n, complex<double> a ) { vr = new complex<double> [lenth =n]; for (int i = 0; i < lenth; i++) vr[i] = a; } Vcr<complex<double> >::Vcr(const Vcr & vec) { vr = new complex<double> [lenth = vec.lenth]; for (int i = 0; i < lenth; i++) vr[i] = vec.vr[i]; } void Vcr<complex<double> >::reset(complex<double> a ) { for (int i = 0; i < lenth; i++) vr[i] = a; } Vcr<complex<double> >& Vcr<complex<double> >::operator=(const Vcr& vec) { if (this != &vec) { if (lenth != vec.lenth ) error("bad vector sizes in vector assignment"); for (int i = 0; i < lenth; i++) vr[i] = vec.vr[i]; } return *this; } Vcr<complex<double> > & Vcr<complex<double> >::operator+=(const Vcr& vec) { if (lenth != vec.lenth ) error("bad vector sizes in vector +="); for (int i = 0; i < lenth; i++) vr[i] += vec.vr[i]; return *this; } Vcr<complex<double> > & Vcr<complex<double> >::operator-=(const Vcr& vec) { if (lenth != vec.lenth ) error("bad vector sizes in vector -="); for (int i = 0; i < lenth; i++) vr[i] -= vec.vr[i]; return *this; } Vcr<complex<double> > operator+(const Vcr<complex<double> >& vec) { return vec; // usage: vec1 = + vec2; } Vcr<complex<double> > operator-(const Vcr<complex<double> >& vec) { Vcr<complex<double> > minus = vec; // usage: vec1 = - vec2; for (int i = 0; i < vec.lenth; i++) minus[i] = - vec[i]; return minus; } Vcr<complex<double> > operator+(const Vcr<complex<double> >& v1, const Vcr<complex<double> >& v2) { // v=v1+v2 if (v1.lenth != v2.lenth ) error("bad vector sizes in vector addition"); Vcr<complex<double> > sum = v1; sum += v2; return sum; } Vcr<complex<double> > operator-(const Vcr<complex<double> >& v1, const Vcr<complex<double> > & v2) { // v=v1-v2 if (v1.lenth != v2.lenth ) error("bad vector sizes in vector subtraction"); Vcr<complex<double> > sum = v1; sum -= v2; return sum; } ostream & operator<<(ostream & s, const Vcr<complex<double> > & vec ) { for (int i =0; i < vec.lenth; i++ ) { s << vec[i] << " "; if (i%10 == 9) s << "\n"; } return s; } double Vcr<complex<double> >::twonorm() const { double nm = abs(vr[0]); for (int i = 1; i < lenth; i++) { double vi = abs(vr[i]); if (nm < 100) nm = sqrt(nm*nm + vi*vi); else { // to avoid overflow for fn "sqrt" when nm is large double tm = vi/nm; nm *= sqrt(1.0 + tm*tm); } } return nm; } double Vcr<complex<double> >::maxnorm() const { double nm = abs(vr[0]); for (int i = 1; i < lenth; i++) { double vi = abs(vr[i]); if (nm < vi) nm = vi; } return nm; } complex<double> dot(const Vcr<complex<double> > & v1, const Vcr<complex<double> > & v2) { if (v1.lenth != v2.lenth ) error("bad vector sizes in vector dot product"); complex<double> tm = v1[0]*conj(v2[0]); for (int i = 1; i < v1.lenth; i++) tm += v1[i]*conj(v2[i]); return tm; } Vcr<complex<double> > operator*(complex<double> scalar, const Vcr<complex<double> > & vec) { Vcr<complex<double> > tm(vec.lenth); for (int i = 0; i < vec.lenth; i++) tm[i] = scalar*vec[i]; return tm; } Vcr<complex<double> > operator*(const Vcr<complex<double> >& vec, complex<double> scalar) { Vcr<complex<double> > tm(vec.lenth); for (int i = 0; i < vec.lenth; i++) tm[i] = scalar*vec[i]; return tm; } Vcr<complex<double> > operator*(const Vcr<complex<double> > & vec1, const Vcr<complex<double> > & vec2) { if (vec1.lenth != vec2.lenth ) error("bad vector sizes in vector mulitply"); Vcr<complex<double> > tm(vec1.lenth); for (int i = 0; i < vec1.lenth; i++) tm[i] = vec1[i]*vec2[i]; return tm; } Vcr<complex<double> > operator/(const Vcr<complex<double> > & vec, complex<double> scalar) { if (abs(scalar) == 0) error("divisor is zero in vector-scalar division"); Vcr< complex<double> > tm(vec.lenth); for (int i = 0; i < vec.lenth; i++) tm[i] = vec[i]/scalar; return tm; } complex<double> dot(complex<double>* a, complex<double>* b, int k) { complex<double> init = 0; for (int i = 0; i < k; i++) init += *a++ * conj(*b++); return init; }
matvec.h
View Code
/******************** 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 2001-08-13 **********************************************************/ #ifndef MATVEC_H #define MATVEC_H #include <cstdlib> #include <cmath> #include <iostream> #include <complex> #include <algorithm> using namespace std; template<class T> class Vcr; // forward declaration template<> class Vcr< complex<double> >; // forward declaration 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 template <class S> friend FullMtx<S> operator-(const FullMtx<S>&); // unary -, m1 = -m2; template <class S> friend FullMtx<S> operator-(const FullMtx<S>&, const FullMtx<S>&); // binary -, m=m1-m2 template <class S> friend ostream& operator<<(ostream&, const FullMtx<S>&); // 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 template <class S> friend Vcr<S> operator+(const Vcr<S>&); // unary operator, v = + v2 template <class S> friend Vcr<S> operator-(const Vcr<S>&); // unary operator, v = - v2 template <class S> friend Vcr<S> operator+(const Vcr<S>&, const Vcr<S>&); // v= v1 + v2 template <class S> friend Vcr<S> operator-(const Vcr<S>&, const Vcr<S>&); // v= v1 - v2 template <class S> friend S dot(const Vcr<S>&, const Vcr<S>&); // dot product template <class S> friend Vcr<S> operator*(S, const Vcr<S>&); // vec-scalar multiply template <class S> friend Vcr<S> operator*(const Vcr<S>&, S); // vec-scalar multiply template <class S> friend Vcr<S> operator*(const Vcr<S>&, const Vcr<S>&); // component-wise multiply template <class S> friend Vcr<S> operator/(const Vcr<S>&, S); // vec-scalar division template <class S> friend ostream& operator<<(ostream&, const Vcr<S>&); }; template<class T> T dot(T* a, T* b, int n); // a vector class template<> 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, ostream& s = cout) { //template<class T> void error(const T& t) { cout << t << ". program exited." << "\n"; exit(1); } const double SMALL = 1.0e-15; //const double SMALL = numeric_limits<double>::epsilon(); // *** 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; } 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.0) error("pivot is zero in FullMtx::GaussElim()"); for (int i = k + 1; i < nrows; i++) { if (tmpx[i][k] != 0.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.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.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; } template<class T> Vcr<T> operator*(const Vcr<T>& vec, T scalar) { Vcr<T> tm(vec.lenth); for (int i = 0; i < vec.lenth; i++) tm[i] = scalar*vec[i]; return tm; } template<class T> Vcr<T> operator*(const Vcr<T> & vec1, const Vcr<T> & vec2) { if (vec1.lenth != vec2.lenth ) error("bad vector sizes in vector multiply"); Vcr<T> tm(vec1.lenth); for (int i = 0; i < vec1.lenth; i++) tm[i] = vec1[i]*vec2[i]; return tm; } template<class T> Vcr<T> operator/(const Vcr<T> & vec, T scalar) { if (scalar == 0) error("divisor is zero in vector-scalar division"); Vcr<T> tm(vec.lenth); for (int i = 0; i < vec.lenth; i++) tm[i] = vec[i]/scalar; return tm; } template<class T> T dot(const Vcr<T> & v1, const Vcr<T> & v2) { if (v1.lenth != v2.lenth ) error("bad vector sizes in dot product"); T tm = v1[0]*v2[0]; for (int i = 1; i < v1.lenth; i++) tm += v1[i]*v2[i]; return tm; } template<class T> T dot(T* a, T* b, int n) { T init = 0; for (int i = 0; i < n; i++) init += *a++ * *b++; return init; } // *** definitions of iterative solvers for AbsMtx template<class T> int AbsMtx<T>::CG(Vcr<T>& x, const Vcr<T>& b, double& eps, int& iter, int prec) { // Conjugate gradient method. // x: on entry, initial guess; on retrun: approximate solution // b: right hand side vector as in A x = b; // eps: on entry, tolerance; on retrun: absolute residual in Euclid norm // iter: on entry, max number of iterations allowed; // on return, actual number of iterations used // prec= 0 if no preconditioning, 1 if diag prec, 2 if SSOR prec if (nrows != b.size()) error("matrix and vector sizes do not match"); const int maxiter = iter; Vcr<T> r = b - (*this)*x; // initial residual Vcr<T> z = preconding(r,prec); // solve the precond system Vcr<T> p = z; // p: search direction T zr = dot(z,r); // inner prod of z and r const double stp = eps*b.twonorm(); // stopping criterion if (r.maxnorm() == 0.0) { // if intial guess is true soln, eps = 0.0; // return. Otherwise division by iter = 0; // zero would occur. return 0; } for (iter = 0; iter < maxiter; iter++) { // main loop of CG method Vcr<T> mp = (*this)*p; // one matrix-vector multiply T pap = dot(mp,p); // one of two inner products T alpha = zr/pap; // pap is 0 only when r is 0 x += alpha*p; // update the iterative soln r -= alpha*mp; // update residual if (r.twonorm() <= stp) break; // stop if convergence achieved z = preconding(r,prec); // preconditioning T zrold = zr; zr = dot(z,r); // another of two inner products T beta= zr/zrold; // zrold = 0 only when r is 0 p = z + beta*p; // update search direction } eps = r.twonorm(); if (iter == maxiter) return 1; else return 0; } // end CG() template<class T> int AbsMtx<T>::GMRES(Vcr<T>& x, const Vcr<T>& b, double& eps, int& iter, int pcn, int m) { // Restarted GMRES for solving Ax = b, called GMRES(m) // Return 0 if sucessful and 1 if breakdown // It outputs the approximate solution, residual, and // number of iterations taken. // // 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 const int maxiter = iter; const double stp = (preconding(b, pcn)).twonorm() * eps; Vcr<T> r = preconding(b - (*this) * x, pcn); T beta = r.twonorm(); // T may be complex<C> bool conv = false; if (m > nrows) error("In GMRES(m), m is bigger than number of rows"); if (abs(beta) <= stp) { // return if initial guess eps = abs(beta); // is true solution iter = 0; return 0; } // orthonormal basis for Krylov space, // v[i] is a pointer to ith basis vector Vcr<T>** v = new Vcr<T>* [m+1]; for (int i = 0; i <= m; i++) v[i] = new Vcr<T>(nrows); // ith orthonormal basis // Hessenburg matrix h, (m+1) by m; // h[i] stores ith column of h that has i+2 entries. // Only non-zero part of h is stored T** h = new T* [m]; for (int i = 0; i < m; i++) h[i] = new T [i+2]; iter = 1; while (iter <= maxiter) { // iterations for gmres(m) *v[0] = r / beta; Vcr<T> g(m+1); // vector in least squares problem g[0] = beta; Vcr<T> cs(m), sn(m); // Givens rotations int k; for (k = 0; k < m && iter <= maxiter; k++, iter++) { // orthogonalization Vcr<T> w = preconding((*this) * (*v[k]), pcn); T nmw = w.twonorm(); for (int i = 0; i <= k; i++) { h[k][i] = dot(w, *v[i]); w -= h[k][i] * (*v[i]); } h[k][k+1] = w.twonorm(); // if h[k][k+1] is small, do reorthogonalization if (nmw + 1.0e-4*h[k][k+1] == nmw) { for (int i = 0; i <= k; i++) { T hri = dot(w, *v[i]); // re-orthogonalization h[k][i] += hri; w -= hri * (*v[i]); } h[k][k+1] = w.twonorm(); } if (h[k][k+1] == 0.0) error("unexpected zero-divisor in GMRES()"); *v[k+1] = w / h[k][k+1]; // apply Givens G_0, G_1, ...,G_{k-1} to column k of h for (int i = 0; i < k; i++) { // T tmp = conj(cs[i])*h[k][i] + conj(sn[i])*h[k][i+1]; T tv[2] = { cs[i], sn[i] }; T tmp = dot(&h[k][i], tv, 2); h[k][i+1] = - sn[i]*h[k][i] + cs[i]*h[k][i+1]; h[k][i] = tmp; } // generate Givens rotation G_k from kth column of h if (h[k][k+1] == 0.0) { cs[k] = 1; sn[k] = 0; } else { T tpm = sqrt(dot(&h[k][k], &h[k][k], 2)); cs[k] = h[k][k]/tpm; sn[k] = h[k][k+1]/tpm; } // apply Givens rotation G_k to kth column of h and to g T tv[2] = { cs[k], sn[k] }; h[k][k] = dot(&h[k][k], tv, 2); T tmp = dot(&g[k], tv, 2); g[k+1] = - sn[k]*g[k] + cs[k]*g[k+1]; g[k] = tmp; if (abs(g[k+1]) <= stp) { // stop if residual small k++; break; } } // back solve the upper triangular system for (int i = k-1; i >= 0; i--){ for (int j = i+1; j < k; j++) g[i] -= h[j][i]*g[j]; g[i] /= h[i][i]; } // 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 = 0) { 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; 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] = 0; } for (int i = 0; i< nrows; i++) { // fm may be non-symmetric 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 preconditioning const T omega = 1; // SSOR parameter for preconditioning 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 in BandMtx::preconding() not implemented"); } } // BandGauss template<class T> void BandMtx<T>::GaussElim(Vcr<T>& bb) const { // Gauss elim without pivoting 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.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.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() // BandGauss 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.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]; } } } pvt[nrowsmone] = nrowsmone; // 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() // BandGauss elimination with partial pivoting template<class T> void BandMtx<T>::GaussElimPP2(Vcr<T>& bb) const { if (nrows != bb.size()) error("matrix-vector sizes not match"); // Gauss eliminations with paritial pivoting // for band matrices. This version interchanges rows // in L and U and does not have any advantage since // it requires more storage than using full matrix // storage format. A straight forward programming without // thinking leads to this version. // This version should never be used. instead use // BandMtx<T>::GaussElimPP(), which requires more thinking // to understand. BandMtx<T> tx(nrows, nrows-1, 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 for (int k = 0; k < nrows; k++) pvt[k] = k; // LU decomposition with column partial pivoting for (int k = 0; k < tx.nrows - 1; k++) { int kbrow = min(nrows - 1 - k, tx.bwleft); int kbcol = min(nrows - 1 - 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"); // interchange pivot row and k-th row if (pc != k) { swap(pvt[k], pvt[pc]); for (int j = - min(k, tx.bwleft); j <= kbcol; j++) swap(tx[pc][j+k-pc], tx[k][j]); } // now eliminate column entries below main // diagonal entry at row k 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]; } } } Vcr<T> rhv(nrows); for (int k = 0; k < nrows; k++) rhv[k] = bb[pvt[k]]; bb = rhv; // Forward substitution LY = Pb for (int k = 1; k < tx.nrows; k++) { int kb = min(k, tx.bwleft); for (int j = 1; j <= kb; j++) bb[k] -= tx[k][-j]*bb[k-j]; } // Backward substitution for (int k = tx.nrows-1; k>= 0; k--) { int kb = min(tx.nrows -1-k, tx.bwrit); for (int j = 1; j <= kb; j++) bb[k] -= tx[k][j]*bb[k+j]; bb[k] /= tx[k][0]; } } // end GaussElimPP2() // *** definitions for members of class SSMtx template<class T> SSMtx<T>::SSMtx(int n, int m, T* ent, int* con, int* dia) { nrows = n; lenth = m; sra = new T [lenth]; clm = new int [lenth]; dgn = new int [nrows]; for (int i =0; i< lenth; i++) sra[i] = ent[i]; for (int i =0; i< nrows; i++) dgn[i] = dia[i]; for (int i =0; i< lenth; i++) clm[i] = con[i]; } template<class T> SSMtx<T>::SSMtx(int n, int m) { nrows = n; lenth = m; sra = new T [lenth]; clm = new int [lenth]; dgn = new int [nrows]; for (int i =0; i< lenth; i++) sra[i] = 0; for (int i =0; i< nrows; i++) dgn[i] = 0; for (int i =0; i< lenth; i++) clm[i] = 0; } template<class T> SSMtx<T>::SSMtx(const SSMtx & mat) { // copy constructor nrows = mat.nrows; lenth = mat.lenth; sra = new T [lenth]; clm = new int [lenth]; dgn = new int [nrows]; for (int i = 0; i< nrows; i++) dgn[i] = mat.dgn[i]; for (int i = 0; i < lenth; i++) { sra[i] = mat[i]; clm[i] = mat.clm[i]; } } template<class T> SSMtx<T>& SSMtx<T>::operator=(const SSMtx & ssm) { if(nrows !=ssm.nrows || lenth != ssm.lenth) error("bad matrix sizes"); for (int i = 0; i < lenth; i++) { sra[i] = ssm[i]; clm[i] = ssm.clm[i]; } for (int i = 0; i < nrows; i++) dgn[i] = ssm.dgn[i]; return *this; } template<class T> Vcr<T> SSMtx<T>::operator*(const Vcr<T> & vec) const { if (nrows != vec.size()) error("matrix and vector sizes do not match"); Vcr<T> tm(nrows); tm[0] += sra[0]*vec[0]; for (int j = 1; j < nrows; j++) { int m = dgn[j-1]+1; if (clm[m] == 0) tm[0] += sra[m]*vec[j]; } for (int i = 1; i < nrows; i++) { for (int j = dgn[i-1] + 1; j <= dgn[i]; j++) tm[i] += sra[j]*vec[clm[j]]; for (int j = i+1; j < nrows; j++) { for (int m = dgn[j-1]+1; m < dgn[j]; m++) if (clm[m] == i) tm[i] += sra[m]*vec[j]; } } return tm; } template<class T> Vcr<T> SSMtx<T>::preconding(const Vcr<T>& r, int prec) const { // solve preconditioning matrix: Pz = r; // prec: 0 if no preconditioning, 1 if diag prec, 2 if SSOR prec if (prec == 0) { // no preconditioning return r; } else if (prec == 1 ){ // diagonal preconditioning Vcr<T> z(nrows); for (int i = 0; i < sz; i++) z[i] = r[i]/sra[dgn[i]]; return z; } else if (prec == 2) { // symmetric SOR preconditioning const T omega = 1; // SSOR parameter for preconditioning Vcr<T> z(nrows); Vcr<T> tm(nrows); tm[0] = r[0]/sra[0]; for (int i = 0; i < sz; i++) { if (abs(sra[dgn[i]]) <= SMALL) { error("SSOR preconditioner did not work"); break; } T sum = 0; for (int j = dgn[i-1]+1; j < dgn[i]; j++) sum += omega*sra[j]*tm[clm[j]]; tm[i] = (r[i] - sum)/sra[dgn[i]]; } for (int i = sz -1; i >= 0; i--) { T sum = 0; for (int j = i+1; j < sz; j++) { for (int m = dgn[j-1]+1; m < dgn[j]; m++) if (clm[m] == i) sum += omega*sra[m]*z[j]; } z[i] = tm[i] - sum/sra[dgn[i]]; } return z; } } // *** 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> SparseMtx<T>::SparseMtx(const SparseMtx & mat) { // copy constructor 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 in SparseMtx::operator*()"); 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 Pz = r, return z. (P: preconditioner) Vcr<T> SparseMtx<T>::preconding(const Vcr<T>& r, int precn) const { if (precn == 0) { // no preconditioning is used return r; } else if (precn == 1) { // diagonal preconditioning Vcr<T> diag(nrows); // find diagonal entries for (int i = 0; i < nrows; i++) { // geometrical average, more stable for (int j = fnz[i]; j < fnz[i+1]; j++) diag[i] += sra[j]*sra[j]; 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 for 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++) { if (clm[j] == i) diag[i] += sra[j]*sra[j]; // averaging is used else diag[i] += omega*omega*sra[j]*sra[j]; // it is more stable } diag[i] = sqrt(diag[i]); } Vcr<T> z(nrows); for (int i = 0; i < nrows; i++) { // solve lower triangular system T sum = 0; 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--) { // solve upper triangular system T sum = 0; 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 in SparseMtx::preconding not implemented"); } } // end SparseMtx::preconding() #endif //****** end matvec.h ***************/
solpoisson.h
View Code
//求解压力步修改版 2012/11/18 // 验证文章 二维泊松方程的两步预估校正格式 中的基于第二类边界条件的Poisson方程 的算例 #include<cmath> #include <iomanip> #include <fstream> #include "matvec.h" //#include"GS.h" #include"GMRES.h" #define PI 3.1415926 //#define Temp(i,j) Temp[(i)*(N+1)*(M+1)+j] //使用宏定义将二维转换成一维 void Sol_P(int N,int M,double delta_x,double delta_y,double delta_t ,double **ustar ,double **vstar ,double **upie ,double **vpie,double **Pres) { /*cout<<sin(double(N))<<PI<<endl;*/ std::ofstream outPres; outPres.open("F:\\fluid\\test_poissonEqution\\Pres.txt"); //void (*p)(double * ,int , double [],double []); //定义一个函数指针变量 p cout<<"N="<<N<<" "<<"M="<<M<<endl; double **A=new double *[(N+1)*(M+1)]; for(int k=0;k<(N+1)*(M+1);k++) A[k]= new double [(N+1)*(M+1)]; for(int i=0;i<(N+1)*(M+1);i++) for(int j=0;j<(N+1)*(M+1);j++) A[i][j]=0; //cout<<endl; // 输出初始化的 A矩阵 // cout<<" 初始化A矩阵 "<<endl; /* for(int i=0;i<(N+1)*(M+1);i++)m { for(int j=0;j<(N+1)*(M+1);j++) { cout<<A[i][j]<<" "; } cout<<endl; }*/ double *B =new double [(N+1)*(M+1)]; for(int bk=0;bk<(N+1)*(M+1);bk++) { B[bk]=0; // cout<<B[bk]<<" "; } cout<<endl; double *sol_P =new double [(N+1)*(M+1)]; for(int bk=0;bk<(N+1)*(M+1);bk++) { sol_P[bk]=0; // cout<<sol_P[bk]<<" "; } cout<<endl; /*---------------------------------------------------------------------------------------------*/ /*点(0,0)处的离散方程*/ A[0*(M+1)+0][0*(M+1)+0]=(-2*(1/pow(delta_x,2)+1/pow(delta_y,2)))/2; A[0*(M+1)+0][0*(M+1)+1]=(2*(1/pow(delta_y,2)))/2; A[0*(M+1)+0][1*(M+1)+0]=(2*(1/pow(delta_x,2)))/2; B[0*(M+1)+0]=(-13*PI*PI*sin(3*PI*0+PI/4)*sin(2*PI*0+PI/4)+2/delta_x*(3*PI*cos(PI/4)*sin(2*PI*0+PI/4))+2/delta_y*(2*PI*cos(PI/4)*sin(3*PI*0+PI/4)))/2;//C+2C1/delta_x+2C2/delta_y // 向后差分 /*i=0上的点（除去(0,0),(0,M)）*/ for(int j=1;j<=M-1;j++) { int i=0;//第一行 A[i*(M+1)+j][i*(M+1)+j-1]=1/pow(delta_y,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=1/pow(delta_y,2); A[i*(M+1)+j][(i+1)*(M+1)+j]=2.0/pow(delta_x,2); B[i*(M+1)+j]=-13*PI*PI*sin(3*PI*i+PI/4)*sin(2*PI*(j*delta_y)+PI/4)+2/delta_x*(3*PI*cos(PI/4)*sin(2*PI*(j*delta_y)+PI/4)); // 向后差分和中心差分（能用中心差分就用中心差分） } /*点(0,M)处的离散方程*/ A[0*(M+1)+M][0*(M+1)+M-1]=(2*(1/pow(delta_y,2)))/2; A[0*(M+1)+M][0*(M+1)+M]=(-2*(1/pow(delta_x,2)+1/pow(delta_y,2)))/2; A[0*(M+1)+M][1*(M+1)+M]=(2*(1/pow(delta_x,2)))/2; B[0*(M+1)+M]=(-13*PI*PI*sin(3*PI*(0*delta_x)+PI/4)*sin(2*PI*(M*delta_y)+PI/4)+2/delta_x*(3*PI*cos(PI/4)*sin(2*PI*(M*delta_y)+PI/4))-2/delta_y*(2*PI*cos(PI/4)*sin(3*PI*(0*delta_x)+PI/4)))/2; /*i=1~N-1,j=1~M-1 之间的离散方程*/ for(int i=1;i<N;i++) { for(int j=1;j<M;j++) { A[i*(M+1)+j][i*(M+1)+j-1]=(1*(1/pow(delta_y,2)))*2; A[i*(M+1)+j][i*(M+1)+j]=(-2*(1/pow(delta_x,2)+1/pow(delta_y,2)))*2; A[i*(M+1)+j][i*(M+1)+j+1]=(1*(1/pow(delta_y,2)))*2; A[i*(M+1)+j][(i-1)*(M+1)+j]=(1*(1/pow(delta_x,2)))*2; A[i*(M+1)+j][(i+1)*(M+1)+j]=(1*(1/pow(delta_x,2)))*2; B[i*(M+1)+j]=(-13*PI*PI*sin(3*PI*(i*delta_x)+PI/4)*sin(2*PI*(j*delta_y)+PI/4))*2;//对一阶导数采用中心差分 } } /*------------------------------------------------------------------------------------------------------*/ /*点(N,0)处的离散方程*/ A[N*(M+1)+0][N*(M+1)+0]=(-2*(1/pow(delta_x,2)+1/pow(delta_y,2)))/2; A[N*(M+1)+0][N*(M+1)+1]=(2*(1/pow(delta_y,2)))/2; A[N*(M+1)+0][(N-1)*(M+1)+0]=(2*(1/pow(delta_x,2)))/2; B[N*(M+1)+0]=(-13*PI*PI*sin(3*PI*(N*delta_x)+PI/4)*sin(2*PI*(0*delta_y)+PI/4)-2/delta_x*(3*PI*cos(PI/4)*sin(2*PI*(0*delta_y)+PI/4))+2/delta_y*(2*PI*cos(PI/4)*sin(3*PI*(N*delta_x)+PI/4)))/2; /*i=N上的点（除去(0,0),(0,M)）*/ for(int j=1;j<=M-1;j++) { int i=N;//第N行 A[i*(M+1)+j][i*(M+1)+j-1]=1/pow(delta_y,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=1/pow(delta_y,2); A[i*(M+1)+j][(i-1)*(M+1)+j]=2.0/pow(delta_x,2); B[i*(M+1)+j]=-13*PI*PI*sin(3*PI*(N*delta_x)+PI/4)*sin(2*PI*(j*delta_y)+PI/4)-2/delta_x*(3*PI*cos(PI/4)*sin(2*PI*(j*delta_y)+PI/4)); } /*点(N,M)处的离散方程*/ A[N*(M+1)+M][N*(M+1)+M-1]=(2*(1/pow(delta_y,2)))/2; A[N*(M+1)+M][N*(M+1)+M]=(-2*(1/pow(delta_x,2)+1/pow(delta_y,2)))/2; A[N*(M+1)+M][(N-1)*(M+1)+M]=(2*(1/pow(delta_x,2)))/2; B[N*(M+1)+M]=(-13*PI*PI*sin(3*PI*(N*delta_x)+PI/4)*sin(2*PI*(M*delta_y)+PI/4)-2/delta_x*(3*PI*cos(PI/4)*sin(2*PI*(M*delta_y)+PI/4))-2/delta_y*(2*PI*cos(PI/4)*sin(3*PI*(N*delta_x)+PI/4)))/2; /*---------------------------------------------------------------------------------------------------------------------*/ /*j=M上的点（除去(0,M),(N,M)）*/ for(int i=1;i<=N-1;i++) { int j=M;//第N行 A[i*(M+1)+j][(i-1)*(M+1)+j]=1/pow(delta_y,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j-1]=2*(1/pow(delta_y,2)); A[i*(M+1)+j][(i+1)*(M+1)+j]=1.0/pow(delta_x,2); B[i*(M+1)+j]=-13*PI*PI*sin(3*PI*(i*delta_x)+PI/4)*sin(2*PI*(j*delta_y)+PI/4)-2/delta_y*(2*PI*cos(PI/4)*sin(3*PI*(i*delta_x)+PI/4)); } /*---------------------------------------------------------------------------------------------------------------------*/ /*j=0上的点（除去(0,0),(N,0)）*/ for(int i=1;i<=N-1;i++) { int j=0;//第0行 A[i*(M+1)+j][(i-1)*(M+1)+j]=1/pow(delta_x,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=2*(1/pow(delta_y,2)); A[i*(M+1)+j][(i+1)*(M+1)+j]=1.0/pow(delta_x,2); B[i*(M+1)+j]=-13*PI*PI*sin(3*PI*(i*delta_x)+PI/4)*sin(2*PI*(j*delta_y)+PI/4)-2/delta_y*(2*PI*cos(PI/4)*sin(3*PI*(i*delta_x)+PI/4)); } /*---------------------------------------------------------------------------------------------------------------------*/ //double *Temp=new double[((N+1)*(M+1))*((N+1)*(M+1))]; //动态分配内存空间用于存储临时矩阵 ; //将矩阵Temp全部元素初始化为 //注：row为行号， col为列号 ofstream outAMatrix; outAMatrix.open("F:\\fluid\\test_poissonEqution\\AMatrix.txt"); for(int row=0;row<(N+1)*(M+1);row++) { for(int col=0;col<(N+1)*(M+1);col++) { //Temp(row,col)=A[row][col]; outAMatrix<<setprecision(5)<<A[row][col]<<" "; } outAMatrix<<endl; } ofstream outBVector; outBVector.open("F:\\fluid\\test_poissonEqution\\BVector.txt"); for(int i=0;i<(N+1)*(M+1);i++) { outBVector<<setprecision(5)<<B[i]<<" "; outBVector<<endl; } //p=GS; //将函数指针指向函数GS; //(*p)(Temp,(N+1)*(M+1),B,sol_P); //调用GS函数 //GMRES int n_gmres = (N+1)*(M+1), m_gmres = (N+1)*(M+1); //FullMtx<double> mat1(n_gmres ,m_gmres ); cout << "******** Sparse Matrix ******** \n"; SparseMtx<double> sm1(n_gmres,n_gmres*n_gmres); for(int i=0;i<n_gmres;i++) { for(int j=0;j<n_gmres;j++) { //mat1[row][col]=A[row][col]; sm1[i*n_gmres+ j] = A[i][j]; sm1.getclm(i*n_gmres+ j) = j; } sm1.getfnz(i) = i*n_gmres; } sm1.getfnz(n_gmres) = n_gmres*n_gmres; Vcr<double> vectest(n_gmres ) ; //vectest[0] =14;vectest[1] = 18;vectest[2] = 20; for(int row=0;row<(N+1)*(M+1);row++) { vectest[row]=B[row]; } int prec = 1; //预条件算子 int iter = n_gmres ; double eps = 1.0e-14; Vcr<double> vec2(n_gmres ); vec2.reset(); //p=gmres; //将函数指针指向函数gmres; gmres(sm1,vec2, vectest,eps,iter,prec,n_gmres); //调用GMRES函数 cout << "GMRES solution = " << vec2 << "\n"; for(int row=0;row<(N+1)*(M+1);row++) { sol_P[row]=vec2[row]; } //for(int row=0;row<(N+1)*(M+1);row++) //{ // for(int col=0;col<(N+1)*(M+1);col++) // { // A[row][col]=Temp(row,col); // } //} for(int row=1;row<=N;row++) { for(int col=1;col<=M;col++) { Pres[row][col]=sol_P[(row-1)*M+col-1]; //Pres represents Pressure } } //delete Temp;// 删除分配的内存 cout<< "输出求解的矩阵A"<<endl; for(int i=0;i<N+1;i++) { for(int j=0;j<M+1;j++) { cout<<"A["<<i<<"]["<<j<<"]"<< " "; // out<<"A["<<i<<"]["<<j<<"]"<< " "; //out<<setprecision(5)<<A[i][j]<<" "; } } cout<<endl; //out<<setprecision(5)<<endl; for(int i=0;i<(N+1)*(M+1);i++) { for(int j=0;j<(N+1)*(M+1);j++) { cout<<A[i][j]<< " "; //out<<setprecision(5)<<A[i][j]<<" "; } cout<<endl; //out<<setprecision(5)<<endl; } cout<<endl; cout<< "输出求解的矩阵Pres"<<endl; for(int i=0;i<N+1;i++) { for(int j=0;j<M+1;j++) { cout<<Pres[i][j]<< " "; outPres<<setprecision(5)<<Pres[i][j]<<" "; } cout<<endl; outPres<<endl; } /*cout<<"Output the Matrix B!"<<endl; for(int bk=0;bk<(N+1)*(M+1);bk++) { //B[bk]=0; cout<<B[0]<<" "; } cout<<endl;*/ cout<<"Output the ustar!"<<endl; for(int i=N;i>=0;i--) { for(int j=0;j<=M;j++) { cout<<" "<<ustar[i][j]; // out<<setprecision(5)<<ustar(i,j)<<" "; } cout<<endl; //out<<setprecision(5)<<endl; } for(int t=0;t<(N+1)*(M+1);t++) delete []A[t]; delete []A; /*撤销u*空间*/ delete [] B; delete [] sol_P; outPres.close(); outAMatrix.close(); }
testposmain.cpp
View Code
/*the Sum and the DATE:2012-8-18*/ #include<iostream> #include<math.h> #include <iomanip> #include <fstream> //#include "solpos.h" #include "solpoisson.h" using namespace std; void (*Sol_Press)(int ,int ,double ,double ,double ,double ** ,double ** ,double ** ,double **,double **); //定义一个函数指针变量 Sol_Press int main() { std::ofstream out; out.open("F:\\fluid\\test_poissonEqution\\CV.txt"); //double ustar; int N; cout<<"Please the row"<<endl;//输入行 cin>>N; int M; cout<<"Please the lie"<<endl;//输入列 cin>>M; double Length=1.0; //求解区域的长 double Wide=1.0; //求解区域的宽 double delta_x=Length/N; double delta_y=Wide/M; double delta_t=0.01; double Re=1; double **Pres=new double *[N+1];//Pres代表压力 for(int k=0;k<N+1;k++) Pres[k]= new double [M+1]; for(int i=0;i<N+1;i++) for(int j=0;j<M+1;j++) Pres[i][j]=0; double **ustar=new double *[N+1]; for(int k=0;k<N+1;k++) ustar[k]=new double [M+1]; for(int i=0;i<N+1;i++) for(int j=0;j<M+1;j++) ustar[i][j]=0; cout<<endl; double **vstar=new double *[N+1]; for(int k=0;k<N+1;k++) vstar[k]=new double [M+1]; for(int i=0;i<N+1;i++) for(int j=0;j<M+1;j++) vstar[i][j]=0; cout<<endl; /*给u'开辟空间并赋初值*/ double **upie=new double *[N+1]; for(int k=0;k<N+1;k++) upie[k]= new double [M+1]; for(int i=0;i<N+1;i++) for(int j=0;j<M+1;j++) upie[i][j]=0; /*给v'开辟空间并赋初值*/ double **vpie=new double *[N+1]; for(int k=0;k<N+1;k++) vpie[k]= new double [M+1]; for(int i=0;i<N+1;i++) for(int j=0;j<M+1;j++) vpie[i][j]=0; for(int k=1;k<2;k++) { Sol_Press=Sol_P; //将函数指针指向函数压力; (*Sol_Press)(N,M,delta_x,delta_y,delta_t,(double **) ustar,(double **) vstar,(double **) upie,(double **) vpie,(double **) Pres); //调用CV函数 } for(int t=0;t<N+1;t++) delete [] Pres[t]; delete [] Pres; for(int t=0;t<N+1;t++) delete []ustar[t]; delete []ustar; for(int t=0;t<N+1;t++) delete []vstar[t]; delete []vstar; /*撤销u'空间*/ for(int t=0;t<N+1;t++) delete [] upie[t]; delete [] upie; /*撤销v'空间*/ for(int t=0;t<N+1;t++) delete [] vpie[t]; delete [] vpie; out.close(); return 0; }