00001
00002
00003
00004
00005
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
00038
00039
00040
00041 ierr = SysProDefineIntAnnotation(PREPROCESSOR,"notrans"); CHKERRQ(ierr);
00042
00043 ierr = SysProDefineIntAnnotation(PREPROCESSOR,"noparallel"); CHKERRQ(ierr);
00044
00045 ierr = SysProDefineIntAnnotation(PREPROCESSOR,"directsolve"); CHKERRQ(ierr);
00046
00047 ierr = SysProDefineIntAnnotation(PREPROCESSOR,"marked"); CHKERRQ(ierr);
00048
00049 #if defined(PETSC_HAVE_BLOCKSOLVE)
00050
00051
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
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
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
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
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
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
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
00140
00141
00142
00143
00144
00145
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
00160
00161
00162
00163
00164
00165
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
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
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
00236
00237
00238
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,&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
00271
00272
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
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 (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
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(type,B,&Buse); CHKERRQ(ierr);
00331 ierr = LinearSystemSetParts
00332 (newproblem,PETSC_NULL,Buse,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00333
00334 ierr = PetscObjectStateIncrease((PetscObject)A); CHKERRQ(ierr);
00335
00336 *outproblem = (NumericalProblem)newproblem;
00337
00338
00339
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,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 (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
00381 ierr = PetscObjectStateDecrease((PetscObject)A); CHKERRQ(ierr);
00382
00383 ierr = LinearSolutionCopy
00384 ((LinearSolution)old,(LinearSolution)nnew); CHKERRQ(ierr);
00385
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 }