System Preprocessors
|
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,PetscBool *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 PetscBool 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; PetscBool 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,PetscBool overwrite, 00300 NumericalProblem inproblem,NumericalProblem *outproblem, 00301 void *gctx,void **ctx,PetscBool *success) 00302 { 00303 LinearSystem newproblem,problem = (LinearSystem)inproblem; 00304 SalsaTransformObject ksp; KSP method; PetscBool 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,PetscBool 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