SolveCG.h

//===========================================================================
// GoTools - SINTEF Geometry Tools version 1.1
//
// GoTools module: CORE
//
// Copyright (C) 2000-2007 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.        
//===========================================================================
#ifndef _GOSOLVECG_H_
#define _GOSOLVECG_H_


//   -----------------------------------------------------------------------
//      Interface file for class SolveCG
//   -----------------------------------------------------------------------
//
//       Solve the equation system Ax=b where A is a symmetric
//       positive definite matrix using Conjugate Gradient Method.
//
//   -----------------------------------------------------------------------
//    Written by: Vibeke Skytt                            09.99
//    Based on  : PrCG.h written by Mike Floater
//   -----------------------------------------------------------------------

#include <vector>

namespace Go
{


class SolveCG
{
public:

    SolveCG();

    ~SolveCG();

    void attachMatrix(double *gmat, int nn);

    void precondRILU(double relaxfac);

    int solve(double *ex, double *eb, int nn);

    void setTolerance(double tolerance = 1.0e-6)
    {tolerance_ = tolerance;}

    void setMaxIterations(int max_iterations)
    {max_iterations_ = max_iterations;}


private:

    std::vector<double> A_;   // Sparse matrix containing the left side
                              // of the equation system.
    int nn_;           // Size of equation system, i.e. number of unknowns.
    int np_;           // Number of non-zero entries in the equation system.
    std::vector<int> irow_;  // The indexes in A_ and jcol_ of the
                          // first non-zeros of the nn_ rows.
    std::vector<int> jcol_;  // The np_ indexes j of the non-zero elements

    double  tolerance_; // The numerical tolerance deciding if we have reached a solution.
    int     max_iterations_; // The maximal number of iterations to be used by solver.

    // Parameters used in RILU preconditioning.

    std::vector<double> M_;  // Preconditioning matrix.
    double omega_;        // Relaxation parameter.

    std::vector<int> diagonal_;  // Index of diagonal elements in the jcol
    int diagset_; // Whether the index of the diagonal elements has been set.

    template <typename RandomIterator1, typename RandomIterator2>
    void matrixProduct(RandomIterator1 sx, RandomIterator2 sy)
    {
        int kj, ki;
        for(kj=0; kj<nn_; kj++) {
            sy[kj] = 0.0;
            for(ki=irow_[kj]; ki<irow_[kj+1]; ki++) {
                sy[kj] += A_[ki] * sx[jcol_[ki]];
            }
        }
    }

    int getIndex(int ki, int kj);

    void forwBack(double *r, double *s);

    int solveRILU(double *ex, double *eb, int nn);

    // Solve the equation system by conjugate gradient method
    int solveStd(double *ex, double *eb, int nn);

    void printPrecond();        // Print LU factorised preconditioning matrix

};

} // namespacew Go

#endif

---