System Preprocessors
ksp.c
Go to the documentation of this file.
00001 /*! \page ksp The iterative method
00002 
00003   The iterative method is not really a transformation, but it is
00004   the last choice made in a preprocessor loop before the final
00005   solver is called. This means that no new matrix analysis is performed 
00006   after applying this transformation.
00007  */
00008 
00009 #include <stdlib.h>
00010 #include <stdio.h>
00011 #include "string.h"
00012 #include "syspro.h"
00013 #include "sysprotransform.h"
00014 #include "sysprolinear.h"
00015 #include "sysprosuit.h"
00016 #include "anamod.h"
00017 #include "linksp.h"
00018 #include "petscmat.h"
00019 #include "petscpc.h"
00020 #include "petscksp.h"
00021 
00022 #define PREPROCESSOR "ksp"
00023 int gmrescycleid;
00024 
00025 #undef __FUNCT__
00026 #define __FUNCT__ "is_gmres_method"
00027 static PetscErrorCode is_gmres_method(KSPType kspt,PetscTruth *f)
00028 {
00029   SalsaTransformObject ksp; int is; PetscErrorCode ierr;
00030 
00031   PetscFunctionBegin;
00032   ierr = TransformObjectGetByName("ksp",(char*)kspt,&ksp); CHKERRQ(ierr);
00033   ierr = TransformObjectGetIntAnnotation(ksp,"is-gmres",&is,f); CHKERRQ(ierr);
00034   *f = TRUTH(*f && is==1);
00035   PetscFunctionReturn(0);
00036 }
00037 
00038 #undef __FUNCT__
00039 #define __FUNCT__ "setup_ksp_choices"
00040 static PetscErrorCode setup_ksp_choices()
00041 {
00042   SalsaTransform ksp; SalsaTransformObject cur;
00043   PetscErrorCode ierr;
00044 
00045   PetscFunctionBegin;
00046 
00047   ierr = PetscObjectComposedDataRegister(&gmrescycleid); CHKERRQ(ierr);
00048   ierr = TransformGetByName(PREPROCESSOR,&ksp); CHKERRQ(ierr);
00049   ierr = SetSuitabilityContextTester(ksp,&AnaModCheckValidFeatureSet); CHKERRQ(ierr);
00050 
00051   /* symmetric substitutions */
00052   ierr = SysProDefineCharAnnotation(PREPROCESSOR,"symm"); CHKERRQ(ierr);
00053   /* requires transpose */
00054   ierr = SysProDefineIntAnnotation(PREPROCESSOR,"trans"); CHKERRQ(ierr);
00055   /* has restart stuff */
00056   ierr = SysProDefineIntAnnotation(PREPROCESSOR,"is-gmres"); CHKERRQ(ierr);
00057 
00058   /*
00059    * Define our repertoire of KSPs
00060    */
00061 
00062   /*
00063    * GMRES
00064    */
00065   ierr = NewTransformObject(PREPROCESSOR,KSPGMRES,&cur); CHKERRQ(ierr);
00066   ierr = TransformObjectSetExplanation
00067     (cur,"Generalized Minimum Residual Method"); CHKERRQ(ierr);
00068   ierr = TransformObjectIntAnnotate(cur,"is-gmres",1); CHKERRQ(ierr);
00069   ierr = TransformObjectDefineOption(cur,"restart length"); CHKERRQ(ierr);
00070   ierr = TransformObjectAddOption(cur,5); CHKERRQ(ierr);
00071   ierr = TransformObjectAddOption(cur,20); CHKERRQ(ierr);
00072   ierr = TransformObjectAddOption(cur,50); CHKERRQ(ierr);
00073 
00074   /*
00075    * fGMRES
00076    */
00077   ierr = NewTransformObject(PREPROCESSOR,KSPFGMRES,&cur); CHKERRQ(ierr);
00078   ierr = TransformObjectSetExplanation(cur,"flexible GMRES"); CHKERRQ(ierr);
00079   ierr = TransformObjectIntAnnotate(cur,"is-gmres",1); CHKERRQ(ierr);
00080   ierr = TransformObjectDefineOption(cur,"restart length"); CHKERRQ(ierr);
00081   ierr = TransformObjectAddOption(cur,5); CHKERRQ(ierr);
00082   ierr = TransformObjectAddOption(cur,20); CHKERRQ(ierr);
00083   ierr = TransformObjectAddOption(cur,50); CHKERRQ(ierr);
00084   
00085   /*
00086    * lGMRES
00087    */
00088   ierr = NewTransformObject(PREPROCESSOR,KSPLGMRES,&cur); CHKERRQ(ierr);
00089   ierr = TransformObjectSetExplanation(cur,"loose GMRES"); CHKERRQ(ierr);
00090   ierr = TransformObjectIntAnnotate(cur,"is-gmres",1); CHKERRQ(ierr);
00091   ierr = TransformObjectDefineOption(cur,"restart length"); CHKERRQ(ierr);
00092   ierr = TransformObjectAddOption(cur,5); CHKERRQ(ierr);
00093   ierr = TransformObjectAddOption(cur,20); CHKERRQ(ierr);
00094   ierr = TransformObjectAddOption(cur,50); CHKERRQ(ierr);
00095   
00096   /* other methods in no particular order */
00097   ierr = NewTransformObject(PREPROCESSOR,KSPTFQMR,&cur); CHKERRQ(ierr);
00098   ierr = TransformObjectSetExplanation(cur,"Transpose-free QMR"); CHKERRQ(ierr);
00099 
00100   ierr = NewTransformObject(PREPROCESSOR,KSPCGS,&cur); CHKERRQ(ierr);
00101   ierr = TransformObjectSetExplanation
00102     (cur,"Conjugate Gradients Squared"); CHKERRQ(ierr);
00103 
00104   ierr = NewTransformObject(PREPROCESSOR,KSPBICG,&cur); CHKERRQ(ierr);
00105   ierr = TransformObjectSetExplanation
00106     (cur,"Bi-Conjugate Gradients"); CHKERRQ(ierr);
00107   ierr = TransformObjectIntAnnotate(cur,"trans",1); CHKERRQ(ierr);
00108 
00109   ierr = NewTransformObject(PREPROCESSOR,KSPBCGS,&cur); CHKERRQ(ierr);
00110   ierr = TransformObjectSetExplanation
00111     (cur,"Bi-conjugate Gradients Stabilized"); CHKERRQ(ierr);
00112 
00113   ierr = NewTransformObject(PREPROCESSOR,KSPBCGSL,&cur); CHKERRQ(ierr);
00114   ierr = TransformObjectSetExplanation
00115     (cur,"Bi-conjugate Gradients Stabilized(l)"); CHKERRQ(ierr);
00116   ierr = TransformObjectDefineOption(cur,"gmres length"); CHKERRQ(ierr);
00117   ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
00118   ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
00119   ierr = TransformObjectAddOption(cur,4); CHKERRQ(ierr);
00120   ierr = TransformObjectAddOption(cur,5); CHKERRQ(ierr);
00121   ierr = TransformObjectAddOption(cur,6); CHKERRQ(ierr);
00122 
00123   ierr = NewTransformObject(PREPROCESSOR,KSPCGNE,&cur); CHKERRQ(ierr);
00124   ierr = TransformObjectSetExplanation
00125     (cur,"CG on the normal equations"); CHKERRQ(ierr);
00126   ierr = TransformObjectIntAnnotate(cur,"trans",1); CHKERRQ(ierr);
00127   
00128   ierr = NewTransformObject(PREPROCESSOR,KSPCG,&cur); CHKERRQ(ierr);
00129   ierr = TransformObjectSetExplanation
00130     (cur,"plain conjugate gradients"); CHKERRQ(ierr);
00131   {
00132     FeatureSet symmetry;
00133     ierr = NewFeatureSet(&symmetry); CHKERRQ(ierr);
00134     ierr = AddToFeatureSet
00135       (symmetry,"simple","symmetry-snorm",PETSC_NULL); CHKERRQ(ierr);
00136     ierr = AddToFeatureSet
00137       (symmetry,"simple","symmetry-anorm",PETSC_NULL); CHKERRQ(ierr);
00138     ierr = TransformObjectSetSuitabilityFunction
00139       (cur,(void*)symmetry,&onlyforsymmetricproblem); CHKERRQ(ierr);
00140   }
00141 
00142   ierr = NewTransformObject(PREPROCESSOR,KSPPREONLY,&cur); CHKERRQ(ierr);
00143   ierr = TransformObjectSetExplanation
00144     (cur,"Preconditioner application"); CHKERRQ(ierr);
00145 
00146   PetscFunctionReturn(0);
00147 }
00148 
00149 #undef __FUNCT__
00150 #define __FUNCT__ "unset_ksps"
00151 static PetscErrorCode unset_ksps()
00152 {
00153   SalsaTransformObject to; void *ctx;
00154   PetscErrorCode ierr;
00155   PetscFunctionBegin;
00156   ierr = TransformObjectGetByName(PREPROCESSOR,KSPCG,&to); CHKERRQ(ierr);
00157   ierr = TransformObjectGetSuitabilityFunction
00158     (to,&ctx,PETSC_NULL); CHKERRQ(ierr);
00159   ierr = DeleteFeatureSet((FeatureSet)ctx); CHKERRQ(ierr);
00160   PetscFunctionReturn(0);
00161 }
00162 
00163 #undef __FUNCT__
00164 #define __FUNCT__ "disable_ksps"
00165 static PetscErrorCode disable_ksps
00166 (NumericalProblem theproblem,SalsaTransform ksp)
00167 {
00168   SalsaTransformObject curpc;
00169   PetscTruth has_annot,flg; PetscErrorCode ierr;
00170   PetscFunctionBegin;
00171 
00172   {
00173     const char *pct; int pcv;
00174     ierr = PreprocessorGetSetting("pc",&pct,&pcv); CHKERRQ(ierr);
00175     ierr = TransformObjectGetByName("pc",pct,&curpc); CHKERRQ(ierr);
00176   }
00177 
00178   /* we use all KSPs .... */
00179   ierr = TransformObjectsUnmarkAll(ksp); CHKERRQ(ierr);
00180 
00181   /* ... except if the current PC has no transpose ... */
00182   ierr = TransformObjectGetIntAnnotation
00183     (curpc,"notrans",(int*)&flg,&has_annot); CHKERRQ(ierr);
00184   if (flg && has_annot) {
00185     int iksp,n; SalsaTransformObject *objs;
00186     ierr = TransformGetObjects(ksp,&n,&objs); CHKERRQ(ierr);
00187     for (iksp=0; iksp<n; iksp++) {
00188       /* ... then we mark the KSPs that need a transpose */
00189       SalsaTransformObject to = objs[iksp]; int trans;
00190       ierr = TransformObjectGetIntAnnotation
00191         (to,"trans",&trans,&flg); CHKERRQ(ierr);
00192       if (flg && trans) {
00193         ierr = TransformObjectMark(to); CHKERRQ(ierr);
00194       }
00195     }
00196   }
00197 
00198   /* ... and if the pc is LU, then we use only preonly and vv */
00199   ierr = TransformObjectGetIntAnnotation
00200     (curpc,"directsolve",(int*)&flg,&has_annot); CHKERRQ(ierr);
00201   if (has_annot && flg) {
00202     int iksp,n; SalsaTransformObject *objs;
00203     ierr = TransformGetObjects(ksp,&n,&objs); CHKERRQ(ierr);
00204     for (iksp=0; iksp<n; iksp++) {
00205       /* ... and we mark everything but preonly  */
00206       SalsaTransformObject to = objs[iksp]; const char *name;
00207       ierr = TransformObjectGetName(to,&name); CHKERRQ(ierr);
00208       if (strcmp(name,KSPPREONLY)) {
00209         ierr = TransformObjectMark(to); CHKERRQ(ierr);
00210       }
00211     }
00212   } else {
00213     /* if it's not a direct solve */
00214     int iksp,n; SalsaTransformObject *objs;
00215     ierr = TransformGetObjects(ksp,&n,&objs); CHKERRQ(ierr);
00216     for (iksp=0; iksp<n; iksp++) {
00217       /* we disable preonly */
00218       SalsaTransformObject to = objs[iksp]; const char *name;
00219       ierr = TransformObjectGetName(to,&name); CHKERRQ(ierr);
00220       if (!strcmp(name,KSPPREONLY)) {
00221         ierr = TransformObjectMark(to); CHKERRQ(ierr);
00222       }
00223     }
00224   }
00225 
00226   PetscFunctionReturn(0);
00227 }
00228 
00229 #undef __FUNCT__
00230 #define __FUNCT__ "set_ksp_options"
00231 static PetscErrorCode set_ksp_options(SalsaTransformObject tf, int kspv)
00232 {
00233   KSPType kspt;;
00234   int g; PetscTruth f; PetscErrorCode ierr; const char *opt;
00235 
00236   PetscFunctionBegin;
00237 
00238   ierr = TransformObjectGetName(tf,(const char**)&kspt); CHKERRQ(ierr);
00239 
00240   opt = "-ksp_type";
00241   ierr = PetscOptionsHasName(PETSC_NULL,opt,&f); CHKERRQ(ierr);
00242   if (f) {ierr = PetscOptionsClearValue(opt); CHKERRQ(ierr);}
00243   ierr = PetscOptionsSetValue(opt,kspt); CHKERRQ(ierr);
00244 
00245   ierr = TransformObjectGetIntAnnotation(tf,"is-gmres",&g,&f); CHKERRQ(ierr);
00246   if (f && g) {
00247     char val[5];
00248     sprintf(val,"%d",kspv);
00249     opt = "-ksp_gmres_restart";
00250     ierr = PetscOptionsHasName(PETSC_NULL,opt,&f); CHKERRQ(ierr);
00251     if (f) {ierr = PetscOptionsClearValue(opt); CHKERRQ(ierr);}
00252     ierr = PetscOptionsSetValue(opt,val); CHKERRQ(ierr);
00253     goto method_set;
00254   }
00255 
00256   ierr = PetscStrcmp(kspt,KSPBCGSL,&f); CHKERRQ(ierr);
00257   if (f) {
00258     char val[5];
00259     sprintf(val,"%d",kspv);
00260     opt = "-ksp_bcgsl_ell";
00261     ierr = PetscOptionsHasName(PETSC_NULL,opt,&f); CHKERRQ(ierr);
00262     if (f) {ierr = PetscOptionsClearValue(opt); CHKERRQ(ierr);}
00263     ierr = PetscOptionsSetValue(opt,val); CHKERRQ(ierr);
00264     ierr = PetscOptionsSetValue("-ksp_bcgsl_xres","0"); CHKERRQ(ierr);
00265     goto method_set;
00266   }
00267 
00268  method_set:
00269   /*
00270    * GMRES only works with preconditioned norm
00271    */
00272   ierr = PetscStrcmp(kspt,KSPGMRES,&f); CHKERRQ(ierr);
00273   if (f) {
00274     ierr = PetscOptionsSetValue
00275       ("-ksp_norm_type","preconditioned"); CHKERRQ(ierr);
00276   } else {
00277     ierr = PetscOptionsSetValue
00278       ("-ksp_norm_type","unpreconditioned"); CHKERRQ(ierr);
00279   }
00280   /*
00281    * FGMRES only works with right pc
00282    */
00283   ierr = PetscOptionsClearValue("-ksp_left_pc"); CHKERRQ(ierr);
00284   ierr = PetscOptionsClearValue("-ksp_right_pc"); CHKERRQ(ierr);
00285   ierr = PetscOptionsClearValue("-ksp_symmetric_pc"); CHKERRQ(ierr);
00286   ierr = PetscStrcmp(kspt,KSPFGMRES,&f); CHKERRQ(ierr);
00287   if (f) {
00288     ierr = PetscOptionsSetValue("-ksp_right_pc",PETSC_NULL); CHKERRQ(ierr);
00289   } else {
00290     ierr = PetscOptionsSetValue("-ksp_left_pc",PETSC_NULL); CHKERRQ(ierr);
00291   }
00292 
00293   PetscFunctionReturn(0);
00294 }
00295 
00296 #undef __FUNCT__
00297 #define __FUNCT__ "setup_ksp"
00298 static PetscErrorCode setup_ksp
00299 (const char* kspt,int kspv,PetscTruth overwrite,
00300  NumericalProblem inproblem,NumericalProblem *outproblem,
00301  void *gctx,void **ctx,PetscTruth *success)
00302 {
00303   LinearSystem newproblem,problem = (LinearSystem)inproblem;
00304   SalsaTransformObject ksp; KSP method; PetscTruth isgmres; PetscErrorCode ierr;
00305   PetscFunctionBegin;
00306   ierr = TransformObjectGetByName("ksp",(KSPType)kspt,&ksp); CHKERRQ(ierr);
00307   ierr = set_ksp_options(ksp,kspv); CHKERRQ(ierr);
00308 
00309   /* get the ksp context which was created by the pc preprocessor */
00310   ierr = PreprocessorGetContext("pc",(void**)&method); CHKERRQ(ierr);
00311 #if 0
00312   {
00313     PC pc; MPI_Comm comm;
00314     ierr = PreprocessorGetContext("pc",(void**)&pc); CHKERRQ(ierr);
00315     ierr = PetscObjectGetComm((PetscObject)pc,&comm); CHKERRQ(ierr);
00316     ierr = KSPCreate(comm,&method); CHKERRQ(ierr);
00317     ierr = SysProLinearInstallCustomKSPMonitor(method); CHKERRQ(ierr);
00318     /*ierr = KSPSetPC(method,pc); CHKERRQ(ierr);*/
00319     *(KSP*)ctx = method;
00320   }
00321 #endif
00322   ierr = KSPSetFromOptions(method); CHKERRQ(ierr);
00323   ierr = is_gmres_method((KSPType)kspt,&isgmres); CHKERRQ(ierr);
00324   if (isgmres) {
00325     ierr = PetscObjectComposedDataSetInt
00326       ((PetscObject)method,gmrescycleid,kspv); CHKERRQ(ierr);
00327   } else {
00328     ierr = PetscObjectComposedDataSetInt
00329       ((PetscObject)method,gmrescycleid,-1); CHKERRQ(ierr);
00330   }
00331   ierr = PetscPushErrorHandler(&PetscReturnErrorHandler,NULL); CHKERRQ(ierr);
00332   ierr = KSPSetUp(method);
00333   if (ierr) {
00334     *success = PETSC_FALSE ; 
00335   } else {
00336     *success = PETSC_TRUE;
00337     ierr = LinearSystemDuplicatePointers(problem,&newproblem); CHKERRQ(ierr);
00338     *outproblem = (NumericalProblem)newproblem;
00339   }
00340   ierr = PetscPopErrorHandler(); CHKERRQ(ierr)
00341   PetscFunctionReturn(0);
00342 }
00343 
00344 #undef __FUNCT__
00345 #define __FUNCT__ "unset_ksp"
00346 static PetscErrorCode unset_ksp
00347 (const char *kspt,PetscTruth overwrite,void *gctx,void *ctx,
00348  NumericalProblem inproblem,NumericalProblem nextproblem,
00349  NumericalSolution old,NumericalSolution nnew)
00350 {
00351   PetscErrorCode ierr;
00352   PetscFunctionBegin;
00353 #if 0
00354   {
00355     KSP solver = (KSP)ctx;
00356     ierr = KSPDestroy(solver); CHKERRQ(ierr);
00357   }
00358 #endif
00359   ierr = LinearSolutionCopy
00360     ((LinearSolution)old,(LinearSolution)nnew); CHKERRQ(ierr);
00361   /*ierr = LinearSolutionDelete((LinearSolution)old); CHKERRQ(ierr);*/
00362   //ierr = DeleteLinearSystem((LinearSystem)inproblem); CHKERRQ(ierr);
00363   PetscFunctionReturn(0);
00364 }
00365 
00366 #undef __FUNCT__
00367 #define __FUNCT__ "DeclareKSPPreprocessor"
00368 PetscErrorCode DeclareKSPPreprocessor(void)
00369 {
00370   PetscErrorCode ierr;
00371   PetscFunctionBegin;
00372   ierr = DeclarePreprocessor
00373     ("ksp",
00374      &setup_ksp_choices,/* this_preprocessor_setup */
00375      &disable_ksps, /* specific_setup */
00376      0, /* specific unset */
00377      &unset_ksps, /* global_unset */
00378      0,0, /* context create/delete */
00379      &setup_ksp,&unset_ksp /* start/end function */
00380      ); CHKERRQ(ierr);
00381   ierr = PreprocessorSetPreservedCategories
00382     (PREPROCESSOR,"laplace"); CHKERRQ(ierr);
00383   /*
00384   ierr = DeclarePreprocessorIntelligentChoice
00385     ("scaling",&set_intelligent_scaling); CHKERRQ(ierr);
00386   */
00387   PetscFunctionReturn(0);
00388 }
00389 
00390 PetscErrorCode (*custommonitor)(KSP,int,PetscReal,void*) = NULL;
00391 void *monitordata = NULL;
00392 
00393 #undef __FUNCT__
00394 #define __FUNCT__ "SysProLinearInstallCustomKSPMonitor"
00395 PetscErrorCode SysProLinearInstallCustomKSPMonitor(KSP solver)
00396 {
00397   PetscErrorCode ierr;
00398   PetscFunctionBegin;
00399   if (custommonitor) {
00400     ierr = KSPMonitorSet
00401       (solver,custommonitor,monitordata,PETSC_NULL); CHKERRQ(ierr);
00402   }
00403   PetscFunctionReturn(0);
00404 }
00405 
00406 #undef __FUNCT__
00407 #define __FUNCT__ "SysProLinearDeclareCustomKSPMonitor"
00408 PetscErrorCode SysProLinearDeclareCustomKSPMonitor
00409 (PetscErrorCode(*monitor)(KSP,int,PetscReal,void*),void *data)
00410 {
00411   PetscFunctionBegin;
00412   custommonitor = monitor;
00413   if (data) monitordata = data;
00414   PetscFunctionReturn(0);
00415 }
00416