//=========================================================================== // GoTools Core - SINTEF Geometry Tools Core library, version 2.0.1 // // Copyright (C) 2000-2007, 2010 SINTEF ICT, Applied Mathematics, Norway. // // This program is free software; you can redistribute it and/or // modify it under the terms of the GNU General Public License // as published by the Free Software Foundation version 2 of the License. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program; if not, write to the Free Software // Foundation, Inc., // 59 Temple Place - Suite 330, // Boston, MA 02111-1307, USA. // // Contact information: E-mail: tor.dokken@sintef.no // SINTEF ICT, Department of Applied Mathematics, // P.O. Box 124 Blindern, // 0314 Oslo, Norway. // // Other licenses are also available for this software, notably licenses // for: // - Building commercial software. // - Building software whose source code you wish to keep private. //===========================================================================
00001 #ifndef _SOLVECG_H_ 00002 #define _SOLVECG_H_ 00003 00004 00005 // ----------------------------------------------------------------------- 00006 // Interface file for class SolveCG 00007 // ----------------------------------------------------------------------- 00008 // 00009 // Solve the equation system Ax=b where A is a symmetric 00010 // positive definite matrix using Conjugate Gradient Method. 00011 // 00012 // ----------------------------------------------------------------------- 00013 // Written by: Vibeke Skytt 09.99 00014 // Based on : PrCG.h written by Mike Floater 00015 // ----------------------------------------------------------------------- 00016 00017 #include <vector> 00018 00019 namespace Go 00020 { 00021 00024 class SolveCG 00025 { 00026 public: 00027 00029 SolveCG(); 00030 00032 ~SolveCG(); 00033 00041 void attachMatrix(double *gmat, int nn); 00042 00045 virtual void precondRILU(double relaxfac); 00046 00053 int solve(double *ex, double *eb, int nn); 00054 00057 void setTolerance(double tolerance = 1.0e-6) 00058 {tolerance_ = tolerance;} 00059 00062 void setMaxIterations(int max_iterations) 00063 {max_iterations_ = max_iterations;} 00064 00065 00066 protected: 00067 00068 std::vector<double> A_; // Sparse matrix containing the left side 00069 // of the equation system. 00070 int nn_; // Size of equation system, i.e. number of unknowns. 00071 int np_; // Number of non-zero entries in the equation system. 00072 std::vector<int> irow_; // The indexes in A_ and jcol_ of the 00073 // first non-zeros of the nn_ rows. 00074 std::vector<int> jcol_; // The np_ indexes j of the non-zero elements 00075 00076 double tolerance_; // The numerical tolerance deciding if we have reached a solution. 00077 int max_iterations_; // The maximal number of iterations to be used by solver. 00078 00079 // Parameters used in RILU preconditioning. 00080 00081 std::vector<double> M_; // Preconditioning matrix. 00082 double omega_; // Relaxation parameter. 00083 00084 std::vector<int> diagonal_; // Index of diagonal elements in the jcol 00085 int diagset_; // Whether the index of the diagonal elements has been set. 00086 00090 template <typename RandomIterator1, typename RandomIterator2> 00091 void matrixProduct(RandomIterator1 sx, RandomIterator2 sy) 00092 { 00093 int kj, ki; 00094 for(kj=0; kj<nn_; kj++) { 00095 sy[kj] = 0.0; 00096 for(ki=irow_[kj]; ki<irow_[kj+1]; ki++) { 00097 sy[kj] += A_[ki] * sx[jcol_[ki]]; 00098 } 00099 } 00100 } 00101 00103 int getIndex(int ki, int kj); 00104 00109 void forwBack(double *r, double *s); 00110 00111 // Compute sy = A_^T * sx. 00112 void transposedMatrixProduct(double *sx, double *sy); 00113 00121 int solveRILU(double *ex, double *eb, int nn); 00122 00123 // Solve the equation system by conjugate gradient method 00129 int solveStd(double *ex, double *eb, int nn); 00130 00132 void printPrecond(); // Print LU factorised preconditioning matrix 00133 00134 }; 00135 00136 } // namespacew Go 00137 00138 #endif