00001
00002
00003
00004
00005
00006
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
00050
00051 ierr = SysProDefineCharAnnotation(PREPROCESSOR,"symm"); CHKERRQ(ierr);
00052
00053 ierr = SysProDefineIntAnnotation(PREPROCESSOR,"trans"); CHKERRQ(ierr);
00054
00055 ierr = SysProDefineIntAnnotation(PREPROCESSOR,"is-gmres"); CHKERRQ(ierr);
00056
00057
00058
00059
00060
00061
00062
00063
00064 ierr = NewTransformObject(PREPROCESSOR,KSPGMRES,&cur); CHKERRQ(ierr);
00065 ierr = TransformObjectSetExplanation
00066 (cur,"Generalized Minimum Residual Method"); CHKERRQ(ierr);
00067 ierr = TransformObjectIntAnnotate(cur,"is-gmres",1); CHKERRQ(ierr);
00068 ierr = TransformObjectDefineOption(cur,"restart length"); CHKERRQ(ierr);
00069 ierr = TransformObjectAddOption(cur,5); CHKERRQ(ierr);
00070 ierr = TransformObjectAddOption(cur,20); CHKERRQ(ierr);
00071 ierr = TransformObjectAddOption(cur,50); CHKERRQ(ierr);
00072
00073
00074
00075
00076 ierr = NewTransformObject(PREPROCESSOR,KSPFGMRES,&cur); CHKERRQ(ierr);
00077 ierr = TransformObjectSetExplanation(cur,"flexible GMRES"); CHKERRQ(ierr);
00078 ierr = TransformObjectIntAnnotate(cur,"is-gmres",1); CHKERRQ(ierr);
00079 ierr = TransformObjectDefineOption(cur,"restart length"); CHKERRQ(ierr);
00080 ierr = TransformObjectAddOption(cur,5); CHKERRQ(ierr);
00081 ierr = TransformObjectAddOption(cur,20); CHKERRQ(ierr);
00082 ierr = TransformObjectAddOption(cur,50); CHKERRQ(ierr);
00083
00084
00085
00086
00087 ierr = NewTransformObject(PREPROCESSOR,KSPLGMRES,&cur); CHKERRQ(ierr);
00088 ierr = TransformObjectSetExplanation(cur,"loose GMRES"); CHKERRQ(ierr);
00089 ierr = TransformObjectIntAnnotate(cur,"is-gmres",1); CHKERRQ(ierr);
00090 ierr = TransformObjectDefineOption(cur,"restart length"); CHKERRQ(ierr);
00091 ierr = TransformObjectAddOption(cur,5); CHKERRQ(ierr);
00092 ierr = TransformObjectAddOption(cur,20); CHKERRQ(ierr);
00093 ierr = TransformObjectAddOption(cur,50); CHKERRQ(ierr);
00094
00095
00096 ierr = NewTransformObject(PREPROCESSOR,KSPTFQMR,&cur); CHKERRQ(ierr);
00097 ierr = TransformObjectSetExplanation(cur,"Transpose-free QMR"); CHKERRQ(ierr);
00098
00099 ierr = NewTransformObject(PREPROCESSOR,KSPCGS,&cur); CHKERRQ(ierr);
00100 ierr = TransformObjectSetExplanation
00101 (cur,"Conjugate Gradients Squared"); CHKERRQ(ierr);
00102
00103 ierr = NewTransformObject(PREPROCESSOR,KSPBICG,&cur); CHKERRQ(ierr);
00104 ierr = TransformObjectSetExplanation
00105 (cur,"Bi-Conjugate Gradients"); CHKERRQ(ierr);
00106 ierr = TransformObjectIntAnnotate(cur,"trans",1); CHKERRQ(ierr);
00107
00108 ierr = NewTransformObject(PREPROCESSOR,KSPBCGS,&cur); CHKERRQ(ierr);
00109 ierr = TransformObjectSetExplanation
00110 (cur,"Bi-conjugate Gradients Stabilized"); CHKERRQ(ierr);
00111
00112 ierr = NewTransformObject(PREPROCESSOR,KSPBCGSL,&cur); CHKERRQ(ierr);
00113 ierr = TransformObjectSetExplanation
00114 (cur,"Bi-conjugate Gradients Stabilized(l)"); CHKERRQ(ierr);
00115 ierr = TransformObjectDefineOption(cur,"gmres length"); CHKERRQ(ierr);
00116 ierr = TransformObjectAddOption(cur,2); CHKERRQ(ierr);
00117 ierr = TransformObjectAddOption(cur,3); CHKERRQ(ierr);
00118 ierr = TransformObjectAddOption(cur,4); CHKERRQ(ierr);
00119 ierr = TransformObjectAddOption(cur,5); CHKERRQ(ierr);
00120 ierr = TransformObjectAddOption(cur,6); CHKERRQ(ierr);
00121
00122 ierr = NewTransformObject(PREPROCESSOR,KSPCGNE,&cur); CHKERRQ(ierr);
00123 ierr = TransformObjectSetExplanation
00124 (cur,"CG on the normal equations"); CHKERRQ(ierr);
00125 ierr = TransformObjectIntAnnotate(cur,"trans",1); CHKERRQ(ierr);
00126
00127 ierr = NewTransformObject(PREPROCESSOR,KSPCG,&cur); CHKERRQ(ierr);
00128 ierr = TransformObjectSetExplanation
00129 (cur,"plain conjugate gradients"); CHKERRQ(ierr);
00130 {
00131 FeatureSet symmetry;
00132 ierr = NewFeatureSet(&symmetry); CHKERRQ(ierr);
00133 ierr = AddToFeatureSet
00134 (symmetry,"simple","symmetry-snorm",PETSC_NULL); CHKERRQ(ierr);
00135 ierr = AddToFeatureSet
00136 (symmetry,"simple","symmetry-anorm",PETSC_NULL); CHKERRQ(ierr);
00137 ierr = TransformObjectSetSuitabilityFunction
00138 (cur,(void*)symmetry,&onlyforsymmetricproblem); CHKERRQ(ierr);
00139 }
00140
00141 ierr = NewTransformObject(PREPROCESSOR,KSPPREONLY,&cur); CHKERRQ(ierr);
00142 ierr = TransformObjectSetExplanation
00143 (cur,"Preconditioner application"); CHKERRQ(ierr);
00144
00145 PetscFunctionReturn(0);
00146 }
00147
00148 #undef __FUNCT__
00149 #define __FUNCT__ "unset_ksps"
00150 static PetscErrorCode unset_ksps(NumericalProblem_ *dummy)
00151 {
00152 SalsaTransformObject to; void *ctx;
00153 PetscErrorCode ierr;
00154 PetscFunctionBegin;
00155 ierr = TransformObjectGetByName(PREPROCESSOR,KSPCG,&to); CHKERRQ(ierr);
00156 ierr = TransformObjectGetSuitabilityFunction
00157 (to,&ctx,PETSC_NULL); CHKERRQ(ierr);
00158 ierr = DeleteFeatureSet((FeatureSet)ctx); CHKERRQ(ierr);
00159 PetscFunctionReturn(0);
00160 }
00161
00162 #undef __FUNCT__
00163 #define __FUNCT__ "disable_ksps"
00164 static PetscErrorCode disable_ksps
00165 (NumericalProblem theproblem,SalsaTransform ksp)
00166 {
00167 SalsaTransformObject curpc;
00168 PetscTruth has_annot,flg; PetscErrorCode ierr;
00169 PetscFunctionBegin;
00170
00171 {
00172 char *pct; int pcv;
00173 ierr = PreprocessorGetSetting("pc",&pct,&pcv); CHKERRQ(ierr);
00174 ierr = TransformObjectGetByName("pc",pct,&curpc); CHKERRQ(ierr);
00175 }
00176
00177
00178 ierr = TransformObjectsUnmarkAll(ksp); CHKERRQ(ierr);
00179
00180
00181 ierr = TransformObjectGetIntAnnotation
00182 (curpc,"notrans",(int*)&flg,&has_annot); CHKERRQ(ierr);
00183 if (flg && has_annot) {
00184 int iksp,n; SalsaTransformObject *objs;
00185 ierr = TransformGetObjects(ksp,&n,&objs); CHKERRQ(ierr);
00186 for (iksp=0; iksp<n; iksp++) {
00187
00188 SalsaTransformObject to = objs[iksp]; int trans;
00189 ierr = TransformObjectGetIntAnnotation
00190 (to,"trans",&trans,&flg); CHKERRQ(ierr);
00191 if (flg && trans) {
00192 ierr = TransformObjectMark(to); CHKERRQ(ierr);
00193 }
00194 }
00195 }
00196
00197
00198 ierr = TransformObjectGetIntAnnotation
00199 (curpc,"directsolve",(int*)&flg,&has_annot); CHKERRQ(ierr);
00200 if (has_annot && flg) {
00201 int iksp,n; SalsaTransformObject *objs;
00202 ierr = TransformGetObjects(ksp,&n,&objs); CHKERRQ(ierr);
00203 for (iksp=0; iksp<n; iksp++) {
00204
00205 SalsaTransformObject to = objs[iksp]; char *name;
00206 ierr = TransformObjectGetName(to,&name); CHKERRQ(ierr);
00207 if (strcmp(name,KSPPREONLY)) {
00208 ierr = TransformObjectMark(to); CHKERRQ(ierr);
00209 }
00210 }
00211 } else {
00212
00213 int iksp,n; SalsaTransformObject *objs;
00214 ierr = TransformGetObjects(ksp,&n,&objs); CHKERRQ(ierr);
00215 for (iksp=0; iksp<n; iksp++) {
00216
00217 SalsaTransformObject to = objs[iksp]; char *name;
00218 ierr = TransformObjectGetName(to,&name); CHKERRQ(ierr);
00219 if (!strcmp(name,KSPPREONLY)) {
00220 ierr = TransformObjectMark(to); CHKERRQ(ierr);
00221 }
00222 }
00223 }
00224
00225 PetscFunctionReturn(0);
00226 }
00227
00228 #undef __FUNCT__
00229 #define __FUNCT__ "set_ksp_options"
00230 static PetscErrorCode set_ksp_options(SalsaTransformObject tf, int kspv)
00231 {
00232 KSPType kspt;;
00233 int g; PetscTruth f; PetscErrorCode ierr; char *opt;
00234
00235 PetscFunctionBegin;
00236
00237 ierr = TransformObjectGetName(tf,(char**)&kspt); CHKERRQ(ierr);
00238
00239 opt = "-ksp_type";
00240 ierr = PetscOptionsHasName(PETSC_NULL,opt,&f); CHKERRQ(ierr);
00241 if (f) {ierr = PetscOptionsClearValue(opt); CHKERRQ(ierr);}
00242 ierr = PetscOptionsSetValue(opt,kspt); CHKERRQ(ierr);
00243
00244 ierr = TransformObjectGetIntAnnotation(tf,"is-gmres",&g,&f); CHKERRQ(ierr);
00245 if (f && g) {
00246 char val[5];
00247 sprintf(val,"%d",kspv);
00248 opt = "-ksp_gmres_restart";
00249 ierr = PetscOptionsHasName(PETSC_NULL,opt,&f); CHKERRQ(ierr);
00250 if (f) {ierr = PetscOptionsClearValue(opt); CHKERRQ(ierr);}
00251 ierr = PetscOptionsSetValue(opt,val); CHKERRQ(ierr);
00252 goto method_set;
00253 }
00254
00255 ierr = PetscStrcmp(kspt,KSPBCGSL,&f); CHKERRQ(ierr);
00256 if (f) {
00257 char val[5];
00258 sprintf(val,"%d",kspv);
00259 opt = "-ksp_bcgsl_ell";
00260 ierr = PetscOptionsHasName(PETSC_NULL,opt,&f); CHKERRQ(ierr);
00261 if (f) {ierr = PetscOptionsClearValue(opt); CHKERRQ(ierr);}
00262 ierr = PetscOptionsSetValue(opt,val); CHKERRQ(ierr);
00263 ierr = PetscOptionsSetValue("-ksp_bcgsl_xres","0"); CHKERRQ(ierr);
00264 goto method_set;
00265 }
00266
00267 method_set:
00268
00269
00270
00271 ierr = PetscStrcmp(kspt,KSPGMRES,&f); CHKERRQ(ierr);
00272 if (f) {
00273 ierr = PetscOptionsSetValue
00274 ("-ksp_norm_type","preconditioned"); CHKERRQ(ierr);
00275 } else {
00276 ierr = PetscOptionsSetValue
00277 ("-ksp_norm_type","unpreconditioned"); CHKERRQ(ierr);
00278 }
00279
00280
00281
00282 ierr = PetscOptionsClearValue("-ksp_left_pc"); CHKERRQ(ierr);
00283 ierr = PetscOptionsClearValue("-ksp_right_pc"); CHKERRQ(ierr);
00284 ierr = PetscOptionsClearValue("-ksp_symmetric_pc"); CHKERRQ(ierr);
00285 ierr = PetscStrcmp(kspt,KSPFGMRES,&f); CHKERRQ(ierr);
00286 if (f) {
00287 ierr = PetscOptionsSetValue("-ksp_right_pc",PETSC_NULL); CHKERRQ(ierr);
00288 } else {
00289 ierr = PetscOptionsSetValue("-ksp_left_pc",PETSC_NULL); CHKERRQ(ierr);
00290 }
00291
00292 PetscFunctionReturn(0);
00293 }
00294
00295 #undef __FUNCT__
00296 #define __FUNCT__ "setup_ksp"
00297 static PetscErrorCode setup_ksp
00298 (char *kspt,int kspv,PetscTruth overwrite,
00299 NumericalProblem inproblem,NumericalProblem *outproblem,
00300 void *gctx,void **ctx,PetscTruth *success)
00301 {
00302 LinearSystem newproblem,problem = (LinearSystem)inproblem;
00303 SalsaTransformObject ksp; KSP method; PetscTruth isgmres; PetscErrorCode ierr;
00304 PetscFunctionBegin;
00305 ierr = TransformObjectGetByName("ksp",kspt,&ksp); CHKERRQ(ierr);
00306 ierr = set_ksp_options(ksp,kspv); CHKERRQ(ierr);
00307
00308
00309 ierr = PreprocessorGetContext("pc",(void**)&method); CHKERRQ(ierr);
00310 #if 0
00311 {
00312 PC pc; MPI_Comm comm;
00313 ierr = PreprocessorGetContext("pc",(void**)&pc); CHKERRQ(ierr);
00314 ierr = PetscObjectGetComm((PetscObject)pc,&comm); CHKERRQ(ierr);
00315 ierr = KSPCreate(comm,&method); CHKERRQ(ierr);
00316 ierr = SysProLinearInstallCustomKSPMonitor(method); CHKERRQ(ierr);
00317
00318 *(KSP*)ctx = method;
00319 }
00320 #endif
00321 ierr = KSPSetFromOptions(method); CHKERRQ(ierr);
00322 ierr = is_gmres_method(kspt,&isgmres); CHKERRQ(ierr);
00323 if (isgmres) {
00324 ierr = PetscObjectComposedDataSetInt
00325 ((PetscObject)method,gmrescycleid,kspv); CHKERRQ(ierr);
00326 } else {
00327 ierr = PetscObjectComposedDataSetInt
00328 ((PetscObject)method,gmrescycleid,-1); CHKERRQ(ierr);
00329 }
00330 ierr = PetscPushErrorHandler(&PetscReturnErrorHandler,NULL); CHKERRQ(ierr);
00331 ierr = KSPSetUp(method);
00332 if (ierr) {
00333 *success = PETSC_FALSE ;
00334 } else {
00335 *success = PETSC_TRUE;
00336 ierr = LinearSystemDuplicatePointers(problem,&newproblem); CHKERRQ(ierr);
00337 *outproblem = (NumericalProblem)newproblem;
00338 }
00339 ierr = PetscPopErrorHandler(); CHKERRQ(ierr)
00340 PetscFunctionReturn(0);
00341 }
00342
00343 #undef __FUNCT__
00344 #define __FUNCT__ "unset_ksp"
00345 static PetscErrorCode unset_ksp
00346 (char *kspt,PetscTruth overwrite,void *gctx,void *ctx,
00347 NumericalProblem inproblem,NumericalProblem nextproblem,
00348 NumericalSolution old,NumericalSolution nnew)
00349 {
00350 PetscErrorCode ierr;
00351 PetscFunctionBegin;
00352 #if 0
00353 {
00354 KSP solver = (KSP)ctx;
00355 ierr = KSPDestroy(solver); CHKERRQ(ierr);
00356 }
00357 #endif
00358 ierr = LinearSolutionCopy
00359 ((LinearSolution)old,(LinearSolution)nnew); CHKERRQ(ierr);
00360
00361 ierr = DeleteLinearSystem((LinearSystem)inproblem); CHKERRQ(ierr);
00362 PetscFunctionReturn(0);
00363 }
00364
00365 #undef __FUNCT__
00366 #define __FUNCT__ "DeclareKSPPreprocessor"
00367 PetscErrorCode DeclareKSPPreprocessor(void)
00368 {
00369 PetscErrorCode ierr;
00370 PetscFunctionBegin;
00371 ierr = DeclarePreprocessor
00372 ("ksp",
00373 &setup_ksp_choices,&disable_ksps,&unset_ksps,0,
00374 0,0,
00375 &setup_ksp,&unset_ksp); CHKERRQ(ierr);
00376 ierr = PreprocessorSetPreservedCategories
00377 (PREPROCESSOR,"laplace"); CHKERRQ(ierr);
00378
00379
00380
00381
00382 PetscFunctionReturn(0);
00383 }
00384
00385 PetscErrorCode (*custommonitor)(KSP,int,PetscReal,void*) = NULL;
00386 void *monitordata = NULL;
00387
00388 #undef __FUNCT__
00389 #define __FUNCT__ "SysProLinearInstallCustomKSPMonitor"
00390 PetscErrorCode SysProLinearInstallCustomKSPMonitor(KSP solver)
00391 {
00392 PetscErrorCode ierr;
00393 PetscFunctionBegin;
00394 if (custommonitor) {
00395 ierr = KSPMonitorSet
00396 (solver,custommonitor,monitordata,PETSC_NULL); CHKERRQ(ierr);
00397 }
00398 PetscFunctionReturn(0);
00399 }
00400
00401 #undef __FUNCT__
00402 #define __FUNCT__ "SysProLinearDeclareCustomKSPMonitor"
00403 PetscErrorCode SysProLinearDeclareCustomKSPMonitor
00404 (PetscErrorCode(*monitor)(KSP,int,PetscReal,void*),void *data)
00405 {
00406 PetscFunctionBegin;
00407 custommonitor = monitor;
00408 if (data) monitordata = data;
00409 PetscFunctionReturn(0);
00410 }
00411