System Preprocessors
pc.c
Go to the documentation of this file.
00001 /*! \page pc The preconditioner
00002 
00003   Choosing a preconditioner changes a linear system into a preconditioned
00004   system. However, this is not a transformation of any coefficient matrix,
00005   so this preprocessor is handled a bit differently from the previous ones.
00006  */
00007 
00008 #include <stdlib.h>
00009 #include <stdio.h>
00010 #include "syspro.h"
00011 #include "syspro_impl.h"
00012 #include "sysprolinear.h"
00013 #include "sysprotransform.h"
00014 #include "linear_impl.h"
00015 #include "linpc.h"
00016 #include "linksp.h"
00017 #include "petsc.h"
00018 #include "petscmat.h"
00019 #include "petscpc.h"
00020 #include "petscksp.h"
00021 #include "anamod.h"
00022 
00023 #define PREPROCESSOR "pc"
00024 
00025 #undef __FUNCT__
00026 #define __FUNCT__ "setup_pc_choices"
00027 static PetscErrorCode setup_pc_choices()
00028 {
00029   SalsaTransform pc; SalsaTransformObject cur;
00030   PetscErrorCode ierr;
00031 
00032   PetscFunctionBegin;
00033 
00034   ierr = TransformGetByName(PREPROCESSOR,&pc); CHKERRQ(ierr);
00035 
00036   /*
00037    * Set up various properties of PCs; these can be used
00038    * later to disable KSP choices
00039    */
00040   /* symmetric substitutions */
00041   ierr = SysProDefineIntAnnotation(PREPROCESSOR,"notrans"); CHKERRQ(ierr);
00042   /* requires transpose */
00043   ierr = SysProDefineIntAnnotation(PREPROCESSOR,"noparallel"); CHKERRQ(ierr);
00044   /* is direct solve */
00045   ierr = SysProDefineIntAnnotation(PREPROCESSOR,"directsolve"); CHKERRQ(ierr);
00046   /* marked as breaking down */
00047   ierr = SysProDefineIntAnnotation(PREPROCESSOR,"marked"); CHKERRQ(ierr);
00048   
00049 #if defined(PETSC_HAVE_BLOCKSOLVE)
00050   /*
00051    * BlockSolve 95
00052    */
00053   ierr = NewTransformObject(PREPROCESSOR,PCBS95,&cur); CHKERRQ(ierr);
00054   ierr = TransformObjectSetExplanation(cur,"BlockSolve 95"); CHKERRQ(ierr);
00055   ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
00056 #endif
00057   
00058 #if defined(PETSC_HAVE_HYPRE)
00059   /*
00060    * Hypre methods: boomeramg and pilut
00061    */
00062   ierr = NewTransformObject(PREPROCESSOR,PCBOOMERAMG,&cur); CHKERRQ(ierr);
00063   ierr = TransformObjectSetExplanation(cur,"Hypre BoomerAMG"); CHKERRQ(ierr);
00064   ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
00065 
00066   ierr = NewTransformObject(PREPROCESSOR,PCPILUT,&cur); CHKERRQ(ierr);
00067   ierr = TransformObjectSetExplanation(cur,"Hypre pILUt"); CHKERRQ(ierr);
00068   ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
00069 
00070   ierr = NewTransformObject(PREPROCESSOR,PCEUCLID,&cur); CHKERRQ(ierr);
00071   ierr = TransformObjectSetExplanation(cur,"Hypre Euclid"); CHKERRQ(ierr);
00072   ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
00073   ierr = TransformObjectDefineOption(cur,"fill"); CHKERRQ(ierr);
00074   ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
00075 
00076   ierr = NewTransformObject(PREPROCESSOR,PCPARASAILS,&cur); CHKERRQ(ierr);
00077   ierr = TransformObjectSetExplanation(cur,"Hypre ParaSails"); CHKERRQ(ierr);
00078   ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
00079   ierr = TransformObjectDefineOption(cur,"fill"); CHKERRQ(ierr);
00080   ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
00081   ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
00082   ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
00083 #endif
00084 
00085   /*
00086    * ILU sequential
00087    */
00088   ierr = NewTransformObject(PREPROCESSOR,PCILU,&cur); CHKERRQ(ierr);
00089   ierr = TransformObjectSetExplanation(cur,"Incomplete LU"); CHKERRQ(ierr);
00090   ierr = TransformObjectIntAnnotate(cur,"noparallel",1); CHKERRQ(ierr);
00091   ierr = TransformObjectDefineOption(cur,"fill"); CHKERRQ(ierr);
00092   ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
00093   ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
00094   ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
00095   /* shifted version */
00096   ierr = NewTransformObject(PREPROCESSOR,PCSILU,&cur); CHKERRQ(ierr);
00097   ierr = TransformObjectSetExplanation(cur,"shifted Incomplete LU"); CHKERRQ(ierr);
00098   ierr = TransformObjectIntAnnotate(cur,"noparallel",1); CHKERRQ(ierr);
00099   ierr = TransformObjectDefineOption(cur,"fill"); CHKERRQ(ierr);
00100   ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
00101   ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
00102   ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
00103   
00104 #if 0
00105 #if defined(PETSC_HAVE_SPAI)
00106   /*
00107    * SPAI
00108    */
00109   ierr = NewTransformObject(PREPROCESSOR,PCSPAI,&cur); CHKERRQ(ierr);
00110   ierr = TransformObjectSetExplanation(cur,"Sparse Approximate Inverse"); CHKERRQ(ierr);
00111   ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
00112   ierr = TransformObjectIntAnnotate(cur,"noparallel",1); CHKERRQ(ierr);
00113   ierr = TransformObjectDefineOption(cur,"blocksize"); CHKERRQ(ierr);
00114   ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
00115   ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
00116   ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
00117 #endif
00118 #endif
00119 
00120   /*
00121    * Additive Schwarz
00122    */
00123   ierr = NewTransformObject
00124     (PREPROCESSOR,PCASM,&cur); CHKERRQ(ierr);
00125   ierr = TransformObjectSetExplanation(cur,"Additive Schwarz method"); CHKERRQ(ierr);
00126   ierr = TransformObjectDefineOption(cur,"local-method"); CHKERRQ(ierr);
00127   /*
00128   ierr = TransformCurrentItemDefineOption(cur,"local method: -1: LU; $\\geq1$ fill levels$+1$; $\\leq-1$: negative of local CG iterations$+1$","local method"); CHKERRQ(ierr);
00129   */
00130   ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
00131   ierr = TransformObjectAddOptionExplanation(cur,1,"ILU(0)"); CHKERRQ(ierr); 
00132   ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
00133   ierr = TransformObjectAddOptionExplanation(cur,2,"ILU(1)"); CHKERRQ(ierr); 
00134   ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
00135   ierr = TransformObjectAddOptionExplanation(cur,3,"ILU(2)"); CHKERRQ(ierr); 
00136   ierr = TransformObjectAddOption(cur,-1); CHKERRQ(ierr);
00137   ierr = TransformObjectAddOptionExplanation(cur,-1,"LU"); CHKERRQ(ierr); 
00138   /*
00139   ierr = TransformObjectAddOption(cur,-7); CHKERRQ(ierr);
00140   ierr = TransformObjectAddOptionExplanation
00141     (cur,"6 iterations GMRES"); CHKERRQ(ierr); 
00142   */
00143 
00144   /*
00145    * Restricted Additive Schwarz
00146    */
00147   ierr = NewTransformObject(PREPROCESSOR,PCRASM,&cur); CHKERRQ(ierr);
00148   ierr = TransformObjectSetExplanation(cur,"Restrictied Additive Schwarz method"); CHKERRQ(ierr);
00149   ierr = TransformObjectDefineOption(cur,"local-method"); CHKERRQ(ierr);
00150   ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
00151   ierr = TransformObjectAddOptionExplanation(cur,1,"ILU(0)"); CHKERRQ(ierr); 
00152   ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
00153   ierr = TransformObjectAddOptionExplanation(cur,2,"ILU(1)"); CHKERRQ(ierr); 
00154   ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
00155   ierr = TransformObjectAddOptionExplanation(cur,3,"ILU(2)"); CHKERRQ(ierr); 
00156   ierr = TransformObjectAddOption(cur,-1); CHKERRQ(ierr);
00157   ierr = TransformObjectAddOptionExplanation(cur,-1,"LU"); CHKERRQ(ierr); 
00158   /*
00159   ierr = TransformObjectAddOption(cur,-7); CHKERRQ(ierr);
00160   ierr = TransformObjectAddOptionExplanation
00161     (cur,"6 iterations GMRES"); CHKERRQ(ierr); 
00162   */
00163 
00164   /*
00165    * Block Jacobi
00166    */
00167   ierr = NewTransformObject(PREPROCESSOR,PCBJACOBI,&cur); CHKERRQ(ierr);
00168   ierr = TransformObjectSetExplanation(cur,"Block Jacobi"); CHKERRQ(ierr);
00169   ierr = TransformObjectDefineOption(cur,"local-method"); CHKERRQ(ierr);
00170   ierr = TransformObjectIntAnnotate(cur,"notrans",1); CHKERRQ(ierr);
00171   ierr = TransformObjectAddOption(cur,1); CHKERRQ(ierr);
00172   ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
00173   ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
00174   ierr = TransformObjectAddOption(cur,-1); CHKERRQ(ierr);
00175   /*
00176   ierr = TransformObjectAddOption(cur,-7); CHKERRQ(ierr);
00177   */
00178 
00179   ierr = NewTransformObject(PREPROCESSOR,PCJACOBI,&cur); CHKERRQ(ierr);
00180   ierr = TransformObjectSetExplanation(cur,"Point Jacobi"); CHKERRQ(ierr);
00181   ierr = NewTransformObject(PREPROCESSOR,PCNONE,&cur); CHKERRQ(ierr);
00182   ierr = TransformObjectSetExplanation(cur,"None"); CHKERRQ(ierr);
00183   
00184   /*
00185    * Direct methods
00186    */
00187   ierr = NewTransformObject(PREPROCESSOR,PCLU,&cur); CHKERRQ(ierr);
00188   ierr = TransformObjectSetExplanation(cur,"Petsc LU"); CHKERRQ(ierr);
00189   ierr = TransformObjectIntAnnotate(cur,"directsolve",1); CHKERRQ(ierr);
00190 #if defined(PETSC_HAVE_MUMPS)
00191   ierr = NewTransformObject(PREPROCESSOR,PCMUMPS,&cur); CHKERRQ(ierr);
00192   ierr = TransformObjectSetExplanation(cur,"Mumps LU"); CHKERRQ(ierr);
00193   ierr = TransformObjectIntAnnotate(cur,"directsolve",1); CHKERRQ(ierr);
00194 #endif
00195 #if defined(PETSC_HAVE_SPOOLES)
00196   ierr = NewTransformObject(PREPROCESSOR,PCSPOOLES,&cur); CHKERRQ(ierr);
00197   ierr = TransformObjectSetExplanation(cur,"Spooles LU"); CHKERRQ(ierr);
00198   ierr = TransformObjectIntAnnotate(cur,"directsolve",1); CHKERRQ(ierr);
00199 #endif
00200 #if defined(PETSC_HAVE_SUPERLU_DIST)
00201   ierr = NewTransformObject(PREPROCESSOR,PCSUPERLU,&cur); CHKERRQ(ierr);
00202   ierr = TransformObjectSetExplanation(cur,"SuperLU LU"); CHKERRQ(ierr);
00203   ierr = TransformObjectIntAnnotate(cur,"directsolve",1); CHKERRQ(ierr);
00204 #endif
00205 #if defined(PETSC_HAVE_UMFPACK)
00206   ierr = NewTransformObject(PREPROCESSOR,PCUMFPACK,&cur); CHKERRQ(ierr);
00207   ierr = TransformObjectSetExplanation(cur,"Umfpack LU"); CHKERRQ(ierr);
00208   ierr = TransformObjectIntAnnotate(cur,"directsolve",1); CHKERRQ(ierr);
00209 #endif
00210 
00211   PetscFunctionReturn(0);
00212 }
00213 
00214 #undef __FUNCT__
00215 #define __FUNCT__ "disable_pcs"
00216 static PetscErrorCode disable_pcs
00217 (NumericalProblem theproblem,SalsaTransform pc)
00218 {
00219   SalsaTransformObject curpc; int ds;
00220   PetscTruth flg; PetscErrorCode ierr;
00221   PetscFunctionBegin;
00222 
00223   ierr = SysProRetrieveQuantity
00224     (theproblem,"variance","diagonal-sign",(void*)&ds,PETSC_NULL,
00225      &flg); CHKERRQ(ierr);
00226   ierr = TransformObjectGetByName("pc","euclid",&curpc); CHKERRQ(ierr);
00227   if (0) {
00228     ierr = TransformObjectMark(curpc);  CHKERRQ(ierr);}
00229 
00230   PetscFunctionReturn(0);
00231 }
00232 
00233 #undef __FUNCT__
00234 #define __FUNCT__ "pcoptionshandling"
00235 /*! Disable certain preconditioners based on commandline options.
00236   
00237   At the moment this is only disabling of direct solvers if the 
00238   user asks for iterative only.
00239 */
00240 static PetscErrorCode pcoptionshandling()
00241 {
00242   char **names; int nnames,iname; PetscTruth onlyiterative,onlydirect;
00243   PetscErrorCode ierr;
00244   PetscFunctionBegin;
00245   ierr = PetscOptionsHasName
00246     (PETSC_NULL,"-syspro_pc_iterative",&onlyiterative); CHKERRQ(ierr);
00247   ierr = PetscOptionsHasName
00248     (PETSC_NULL,"-syspro_pc_direct",&onlydirect); CHKERRQ(ierr);
00249   if (onlyiterative || onlydirect) { 
00250     ierr = RetrieveAllPreprocessorValues
00251       (PREPROCESSOR,(const char***)&names,&nnames); CHKERRQ(ierr);
00252     for (iname=0; iname<nnames; iname++) {
00253       SalsaTransformObject to; int isdirect; PetscTruth flg;
00254       ierr = TransformObjectGetByName
00255         (PREPROCESSOR,names[iname],&to); CHKERRQ(ierr);
00256       ierr = TransformObjectGetIntAnnotation
00257         (to,"directsolve",&isdirect,&flg); CHKERRQ(ierr);
00258       if ( (flg && isdirect && onlyiterative) ||
00259            ( onlydirect && ( !flg || !isdirect) ) 
00260            ){
00261         ierr = TransformObjectMark(to); CHKERRQ(ierr);
00262       }
00263     }
00264   }
00265   PetscFunctionReturn(0);
00266 }
00267 
00268 #undef __FUNCT__
00269 #define __FUNCT__ "create_solver"
00270 /*! Create a solver and 
00271   install a monitor that dynamically increases the maximum
00272   number of iterations.
00273 */
00274 static PetscErrorCode create_solver(NumericalProblem prob,void **ctx)
00275 {
00276   KSP solver; PetscErrorCode ierr;
00277   PetscFunctionBegin;
00278 
00279   ierr = KSPCreate(prob->comm,&solver); CHKERRQ(ierr);
00280   /*ierr = SysProLinearInstallCustomKSPMonitor(solver); CHKERRQ(ierr);*/
00281   *(KSP*)ctx = solver;
00282 #if 0
00283   { 
00284     PC pc;
00285     ierr = PCCreate(prob->comm,&pc); CHKERRQ(ierr);
00286     *(PC*)ctx = pc;
00287   }
00288 #endif
00289   PetscFunctionReturn(0);
00290 }
00291 
00292 #undef __FUNCT__
00293 #define __FUNCT__ "destroy_solver"
00294 static PetscErrorCode destroy_solver(void *ctx)
00295 {
00296   KSP solver; PetscErrorCode ierr;
00297   PetscFunctionBegin;
00298 
00299   solver = (KSP)ctx;
00300   ierr = KSPDestroy(solver); CHKERRQ(ierr);
00301 #if 0
00302   {
00303     PC pc;
00304     pc = (PC)ctx;
00305     ierr = PCDestroy(pc); CHKERRQ(ierr);
00306   }
00307 #endif
00308   PetscFunctionReturn(0);
00309 }
00310 
00311 #undef __FUNCT__
00312 #define __FUNCT__ "setup_pc"
00313 static PetscErrorCode setup_pc
00314 (const char *type,int pcv,PetscTruth overwrite,
00315  NumericalProblem inproblem,NumericalProblem *outproblem,
00316  void *gctx,void **ctx,PetscTruth *success)
00317 {
00318   LinearSystem newproblem,problem = (LinearSystem)inproblem;
00319   KSP solver; PC pc;  Mat A,B,Buse;
00320   PetscErrorCode ierr;
00321 
00322   PetscFunctionBegin;
00323   /*
00324    * handling the problem data is trivial
00325    */
00326   ierr = LinearSystemDuplicatePointers(problem,&newproblem); CHKERRQ(ierr);
00327 
00328   ierr = LinearSystemGetParts
00329     (problem,&A,&B,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00330   ierr = set_preconditioner_base_matrix((PCType)type,B,&Buse); CHKERRQ(ierr);
00331   ierr = LinearSystemSetParts
00332     (newproblem,PETSC_NULL,Buse,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00333   /*  ierr = PurgeAttachedArrays(A); CHKERRQ(ierr);*/
00334   ierr = PetscObjectStateIncrease((PetscObject)A); CHKERRQ(ierr);
00335 
00336   *outproblem = (NumericalProblem)newproblem;
00337   
00338   /*
00339    * the fun part is setting up the Petsc PC
00340    */
00341   solver = (KSP)gctx;
00342   ierr = KSPSetOperators(solver,A,Buse,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
00343   ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr);
00344 #if 0
00345   pc = (PC)gctx;
00346   ierr = PCSetOperators(pc,A,Buse,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
00347 #endif
00348   ierr = SetPetscOptionsForPC(pc,(PCType)type,pcv,6); CHKERRQ(ierr);
00349   {
00350     PetscLogDouble t1,t2,*tstore; PetscErrorCode ierr_save1,ierr_save2;
00351     ierr = PetscGetTime(&t1); CHKERRQ(ierr);
00352     ierr = PCSetFromOptions(pc); CHKERRQ(ierr);
00353     ierr = PetscPushErrorHandler(&PetscReturnErrorHandler,NULL); CHKERRQ(ierr);
00354     ierr = PCSetUp(pc); ierr_save1 = ierr;
00355     ierr = PCSetUpOnBlocks(pc); ierr_save2 = ierr;
00356     ierr = PetscGetTime(&t2); CHKERRQ(ierr);
00357     ierr = PetscPopErrorHandler(); CHKERRQ(ierr);
00358     *success = TRUTH(!ierr_save1 && !ierr_save2);
00359     ierr = PreprocessorGetContext("solution",(void**)&tstore); CHKERRQ(ierr);
00360     *tstore = t2-t1;
00361   }
00362   
00363   PetscFunctionReturn(0);
00364 }
00365 
00366 #undef __FUNCT__
00367 #define __FUNCT__ "unset_pc"
00368 static PetscErrorCode unset_pc
00369 (const char *type,PetscTruth overwrite,void *gctx,void *ctx,
00370  NumericalProblem thisproblem,NumericalProblem upproblem,
00371  NumericalSolution old,NumericalSolution nnew)
00372 {
00373   LinearSystem problem = (LinearSystem)thisproblem;
00374   Mat A;
00375   PetscErrorCode ierr;
00376 
00377   PetscFunctionBegin;
00378   ierr = LinearSystemGetParts
00379     (problem,&A,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00380   /*  ierr = PurgeAttachedArrays(A); CHKERRQ(ierr);*/
00381   ierr = PetscObjectStateDecrease((PetscObject)A); CHKERRQ(ierr);
00382 
00383   ierr = LinearSolutionCopy
00384     ((LinearSolution)old,(LinearSolution)nnew); CHKERRQ(ierr);
00385   /*ierr = LinearSolutionDelete((LinearSolution)old); CHKERRQ(ierr);*/
00386   //ierr = DeleteLinearSystem((LinearSystem)thisproblem); CHKERRQ(ierr);
00387   PetscFunctionReturn(0);
00388 }
00389 
00390 #undef __FUNCT__
00391 #define __FUNCT__ "DeclarePCPreprocessor"
00392 PetscErrorCode DeclarePCPreprocessor(void)
00393 {
00394   SystemPreprocessor pp; PetscErrorCode ierr;
00395   PetscFunctionBegin;
00396   ierr = DeclarePreprocessor
00397     (PREPROCESSOR,
00398      &setup_pc_choices,&disable_pcs,0,0,
00399      &create_solver,&destroy_solver,
00400      &setup_pc,&unset_pc); CHKERRQ(ierr);
00401   ierr = PreprocessorSetPreservedCategories
00402     (PREPROCESSOR,"laplace"); CHKERRQ(ierr);
00403   ierr = SystemPreprocessorGetByName(PREPROCESSOR,&pp); CHKERRQ(ierr);
00404   pp->optionshandling = &pcoptionshandling;
00405   PetscFunctionReturn(0);
00406 }