System Preprocessors
|
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 PetscBool 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; PetscBool 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; PetscBool 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,PetscBool overwrite, 00315 NumericalProblem inproblem,NumericalProblem *outproblem, 00316 void *gctx,void **ctx,PetscBool *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,PetscBool 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 }