escript  Revision_
Preconditioner.h
Go to the documentation of this file.
00001 
00002 /*****************************************************************************
00003 *
00004 * Copyright (c) 2003-2014 by University of Queensland
00005 * http://www.uq.edu.au
00006 *
00007 * Primary Business: Queensland, Australia
00008 * Licensed under the Open Software License version 3.0
00009 * http://www.opensource.org/licenses/osl-3.0.php
00010 *
00011 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
00012 * Development 2012-2013 by School of Earth Sciences
00013 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
00014 *
00015 *****************************************************************************/
00016 
00017 #ifndef __PASO_PRECONDITIONER_H__
00018 #define __PASO_PRECONDITIONER_H__
00019 
00020 #include "SystemMatrix.h"
00021 #include "BOOMERAMG.h"
00022 
00023 #define PRECONDITIONER_NO_ERROR 0
00024 #define PRECONDITIONER_MAXITER_REACHED 1
00025 #define PRECONDITIONER_INPUT_ERROR -1
00026 #define PRECONDITIONER_MEMORY_ERROR -9
00027 #define PRECONDITIONER_BREAKDOWN -10
00028 #define PRECONDITIONER_NEGATIVE_NORM_ERROR -11
00029 #define PRECONDITIONER_DIVERGENCE -12
00030 
00031 namespace paso {
00032 
00033 struct MergedSolver;
00034 struct Preconditioner;
00035 typedef boost::shared_ptr<Preconditioner> Preconditioner_ptr;
00036 typedef boost::shared_ptr<const Preconditioner> const_Preconditioner_ptr;
00037 
00038 struct Preconditioner_Smoother;
00039 struct Preconditioner_AMG_Root;
00040 struct Solver_ILU;
00041 struct Solver_RILU;
00042 
00043 // general preconditioner interface
00044 struct Preconditioner
00045 {
00046     dim_t type;
00047     dim_t sweeps;
00049     Preconditioner_Smoother* jacobi;
00051     Preconditioner_Smoother* gs;
00053     Preconditioner_AMG_Root* amg;
00055     Solver_ILU* ilu;
00057     Solver_RILU* rilu;
00058 };
00059 
00060 void Preconditioner_free(Preconditioner*);
00061 Preconditioner* Preconditioner_alloc(SystemMatrix_ptr A, Options* options);
00062 void Preconditioner_solve(Preconditioner* prec, SystemMatrix_ptr A, double*, double*);
00063 
00064 
00065 // GAUSS SEIDEL & Jacobi
00066 struct Preconditioner_LocalSmoother
00067 {
00068     bool Jacobi;
00069     double* diag;
00070     double* buffer;
00071     index_t* pivot;
00072 };
00073 
00074 struct Preconditioner_Smoother
00075 {
00076     Preconditioner_LocalSmoother* localSmoother;
00077     bool is_local;
00078 };
00079 
00080 void Preconditioner_Smoother_free(Preconditioner_Smoother * in);
00081 void Preconditioner_LocalSmoother_free(Preconditioner_LocalSmoother * in);
00082 
00083 Preconditioner_Smoother* Preconditioner_Smoother_alloc(
00084         SystemMatrix_ptr A, bool jacobi, bool is_local, bool verbose);
00085 
00086 Preconditioner_LocalSmoother* Preconditioner_LocalSmoother_alloc(
00087         SparseMatrix_ptr A, bool jacobi, bool verbose);
00088 
00089 void Preconditioner_Smoother_solve(SystemMatrix_ptr A,
00090         Preconditioner_Smoother* gs, double* x, const double* b,
00091         dim_t sweeps, bool x_is_initial);
00092 
00093 void Preconditioner_LocalSmoother_solve(SparseMatrix_ptr A,
00094         Preconditioner_LocalSmoother* gs, double* x, const double* b,
00095         dim_t sweeps, bool x_is_initial);
00096 
00097 err_t Preconditioner_Smoother_solve_byTolerance(SystemMatrix_ptr A,
00098         Preconditioner_Smoother* gs, double* x, const double* b,
00099         double atol, dim_t* sweeps, bool x_is_initial);
00100 
00101 void Preconditioner_LocalSmoother_Sweep(SparseMatrix_ptr A,
00102         Preconditioner_LocalSmoother* gs, double* x);
00103 
00104 void Preconditioner_LocalSmoother_Sweep_sequential(
00105         SparseMatrix_ptr A, Preconditioner_LocalSmoother* gs,
00106         double* x);
00107 
00108 void Preconditioner_LocalSmoother_Sweep_tiled(SparseMatrix_ptr A,
00109         Preconditioner_LocalSmoother* gs, double* x);
00110 
00111 void Preconditioner_LocalSmoother_Sweep_colored(SparseMatrix_ptr A,
00112         Preconditioner_LocalSmoother* gs, double* x);
00113 
00114 
00115 typedef enum
00116 {
00117     PASO_AMG_UNDECIDED=-1,
00118     PASO_AMG_IN_F=0,
00119     PASO_AMG_IN_C=1
00120 } AMGBlockSelect;
00121 
00123 struct Preconditioner_AMG
00124 {
00125     dim_t level;
00127     SystemMatrix_ptr A_C;
00129     SystemMatrix_ptr P;
00131     SystemMatrix_ptr R;
00132 
00133     Preconditioner_Smoother* Smoother;
00134     dim_t post_sweeps;
00135     dim_t pre_sweeps;
00137     dim_t options_smoother;
00139     bool verbose;
00141     index_t reordering;
00143     dim_t refinements;
00145     double* r;
00147     double* x_C;
00149     double* b_C;
00151     MergedSolver* merged_solver;
00152     Preconditioner_AMG* AMG_C;
00153 };
00154 
00155 void Preconditioner_AMG_free(Preconditioner_AMG* in);
00156 Preconditioner_AMG* Preconditioner_AMG_alloc(SystemMatrix_ptr A, dim_t level, Options* options);
00157 void Preconditioner_AMG_solve(SystemMatrix_ptr A, Preconditioner_AMG* amg, double* x, double* b);
00158 void Preconditioner_AMG_setStrongConnections(SystemMatrix_ptr A,  dim_t* degree_S, index_t* offset_S, index_t* S, double theta, double tau);
00159 void Preconditioner_AMG_setStrongConnections_Block(SystemMatrix_ptr A, dim_t* degree_S, index_t* offset_S, index_t* S, double theta, double tau);
00160 SystemMatrix_ptr Preconditioner_AMG_getProlongation(SystemMatrix_ptr A_p, const index_t* offset_S, const dim_t* degree_S, const index_t* S, dim_t n_C, index_t* counter_C, index_t interpolation_method);
00161 void Preconditioner_AMG_setClassicProlongation(SystemMatrix_ptr P, SystemMatrix_ptr A, const index_t* offset_S, const dim_t* degree_S, const index_t* S,const index_t *counter_C);
00162 void Preconditioner_AMG_setClassicProlongation_Block(SystemMatrix_ptr P, SystemMatrix_ptr A, const index_t* offset_S, const dim_t* degree_S, const index_t* S, const index_t* counter_C);
00163 void Preconditioner_AMG_setDirectProlongation(SystemMatrix_ptr P, SystemMatrix_ptr A, const index_t* offset_S, const dim_t* degree_S, const index_t* S, const index_t* counter_C);
00164 void Preconditioner_AMG_setDirectProlongation_Block(SystemMatrix_ptr P, SystemMatrix_ptr A, const index_t* offset_S, const dim_t* degree_S, const index_t* S, const index_t* counter_C);
00165 double Preconditioner_AMG_getCoarseLevelSparsity(const Preconditioner_AMG* in);
00166 dim_t Preconditioner_AMG_getNumCoarseUnknowns(const Preconditioner_AMG* in);
00167 index_t Preconditioner_AMG_getMaxLevel(const Preconditioner_AMG* in);
00168 void Preconditioner_AMG_transposeStrongConnections(dim_t n, const dim_t* degree_S, const index_t* offset_S, const index_t* S, dim_t nT, dim_t* degree_ST, index_t* offset_ST, index_t* ST);
00169 void Preconditioner_AMG_CIJPCoarsening(dim_t n, dim_t my_n,
00170         AMGBlockSelect* split_marker, const dim_t* degree_S,
00171         const index_t* offset_S, const index_t* S, const dim_t* degree_ST,
00172         const index_t* offset_ST, const index_t* ST,
00173         Connector_ptr col_connector, const_Distribution_ptr col_dist);
00174 
00175 SystemMatrix_ptr Preconditioner_AMG_getRestriction(SystemMatrix_ptr P);
00176 
00177 SystemMatrix_ptr Preconditioner_AMG_buildInterpolationOperator(
00178         SystemMatrix_ptr A, SystemMatrix_ptr P, SystemMatrix_ptr R);
00179 
00180 SystemMatrix_ptr Preconditioner_AMG_buildInterpolationOperatorBlock(
00181         SystemMatrix_ptr A, SystemMatrix_ptr P, SystemMatrix_ptr R);
00182 
00183 SparseMatrix_ptr Preconditioner_AMG_mergeSystemMatrix(SystemMatrix_ptr A);
00184 
00185 void Preconditioner_AMG_mergeSolve(Preconditioner_AMG* amg);
00186 
00188 struct Preconditioner_LocalAMG
00189 {
00190    dim_t level;
00191    SparseMatrix_ptr A_C;  /* coarse level matrix */
00192    SparseMatrix_ptr P;    /* prolongation n x n_C*/
00193    SparseMatrix_ptr R;    /* restriction  n_C x n */
00194 
00195    Preconditioner_LocalSmoother* Smoother;
00196    dim_t post_sweeps;
00197    dim_t pre_sweeps;
00198    index_t reordering;  /* applied reordering in direct solver */
00199    dim_t refinements;  /* number of refinements in direct solver (typically =0) */
00200    double* r;         /* buffer for residual */
00201    double* x_C;       /* solution of coarse level system */
00202    double* b_C;       /* right hand side of coarse level system */
00203    struct Preconditioner_LocalAMG * AMG_C;
00204 };
00205 
00206 void Preconditioner_LocalAMG_free(Preconditioner_LocalAMG * in);
00207 Preconditioner_LocalAMG* Preconditioner_LocalAMG_alloc(SparseMatrix_ptr A, dim_t level, Options* options);
00208 void Preconditioner_LocalAMG_solve(SparseMatrix_ptr A, Preconditioner_LocalAMG * amg, double * x, double * b);
00209 
00210 void Preconditioner_LocalAMG_RungeStuebenSearch(const dim_t n, const index_t* offset, const dim_t* degree, const index_t* S, AMGBlockSelect*split_marker, const bool usePanel);
00211 void Preconditioner_LocalAMG_setStrongConnections_Block(SparseMatrix_ptr A, dim_t *degree, index_t *S, const double theta, const double tau);
00212 void Preconditioner_LocalAMG_setStrongConnections(SparseMatrix_ptr A, dim_t *degree, index_t *S, const double theta, const double tau);
00213 
00214 SparseMatrix_ptr Preconditioner_LocalAMG_getProlongation(
00215         SparseMatrix_ptr A_p, const index_t* offset_S,
00216         const dim_t* degree_S, const index_t* S, dim_t n_C,
00217         const index_t* counter_C, index_t interpolation_method);
00218 
00219 void Preconditioner_LocalAMG_setDirectProlongation_Block(SparseMatrix_ptr P_p, const_SparseMatrix_ptr A_p, const index_t *counter_C);
00220 
00221 void Preconditioner_LocalAMG_setDirectProlongation(SparseMatrix_ptr P_p, const_SparseMatrix_ptr A_p, const index_t *counter_C);
00222 void Preconditioner_LocalAMG_setClassicProlongation(SparseMatrix_ptr P_p, SparseMatrix_ptr A_p, const index_t* offset_S, const dim_t* degree_S, const index_t* S, const index_t *counter_C);
00223 void Preconditioner_LocalAMG_setClassicProlongation_Block(SparseMatrix_ptr P_p, SparseMatrix_ptr A_p, const index_t* offset_S, const dim_t* degree_S, const index_t* S, const index_t *counter_C);
00224 index_t Preconditioner_LocalAMG_getMaxLevel(const Preconditioner_LocalAMG * in);
00225 double Preconditioner_LocalAMG_getCoarseLevelSparsity(const Preconditioner_LocalAMG * in);
00226 dim_t Preconditioner_LocalAMG_getNumCoarseUnknowns(const Preconditioner_LocalAMG * in);
00227 void Preconditioner_LocalAMG_enforceFFConnectivity(const dim_t n, const index_t* offset_S, const dim_t* degree_S, const index_t* S, AMGBlockSelect*split_marker);
00228 
00229 
00230 struct Preconditioner_AMG_Root
00231 {
00232   bool is_local;
00233   Preconditioner_AMG* amg;
00234   Preconditioner_LocalAMG* localamg;
00235   Preconditioner_BoomerAMG* boomeramg;
00236   dim_t sweeps;
00237   Preconditioner_Smoother* amgsubstitute;
00238 };
00239 
00240 Preconditioner_AMG_Root* Preconditioner_AMG_Root_alloc(SystemMatrix_ptr A, Options* options);
00241 void Preconditioner_AMG_Root_free(Preconditioner_AMG_Root * in);
00242 void Preconditioner_AMG_Root_solve(SystemMatrix_ptr A, Preconditioner_AMG_Root * amg, double * x, double * b);
00243 
00245 struct Solver_ILU
00246 {
00247     double* factors;
00248 };
00249 
00251 struct Solver_RILU
00252 {
00253     dim_t n;
00254     dim_t n_block;
00255     dim_t n_F;
00256     dim_t n_C;
00257     double* inv_A_FF;
00258     index_t* A_FF_pivot;
00259     SparseMatrix_ptr A_FC;
00260     SparseMatrix_ptr A_CF;
00261     index_t* rows_in_F;
00262     index_t* rows_in_C;
00263     index_t* mask_F;
00264     index_t* mask_C;
00265     double* x_F;
00266     double* b_F;
00267     double* x_C;
00268     double* b_C;
00269     Solver_RILU* RILU_of_Schur;
00270 };
00271 
00272 void Solver_ILU_free(Solver_ILU * in);
00273 Solver_ILU* Solver_getILU(SparseMatrix_ptr A, bool verbose);
00274 void Solver_solveILU(SparseMatrix_ptr A, Solver_ILU* ilu, double* x, const double* b);
00275 
00276 void Solver_RILU_free(Solver_RILU* in);
00277 Solver_RILU* Solver_getRILU(SparseMatrix_ptr A, bool verbose);
00278 void Solver_solveRILU(Solver_RILU* rilu, double* x, double* b);
00279 
00280 void Solver_updateIncompleteSchurComplement(SparseMatrix_ptr A_CC,
00281         SparseMatrix_ptr A_CF, double* invA_FF, index_t* A_FF_pivot,
00282         SparseMatrix_ptr A_FC);
00283 
00284 } // namespace paso
00285 
00286 #endif // __PASO_PRECONDITIONER_H__
00287