System Preprocessors
|
00001 #include <stdlib.h> 00002 #include "syspro.h" 00003 #include "sysprotransform.h" 00004 #include "sysprolinear.h" 00005 00006 /* test for the syspro linear package 00007 test for problem solving without preprocessors: 00008 this uses PETSc's KSP for actually solving a linear system 00009 */ 00010 00011 /* this example features a preprocessor that we can loop over: 00012 preconditioners 00013 */ 00014 #undef __FUNCT__ 00015 #define __FUNCT__ "create_solver" 00016 /*! Create a solver and 00017 install a monitor that dynamically increases the maximum 00018 number of iterations. 00019 */ 00020 static PetscErrorCode create_solver(NumericalProblem prob,void **ctx) 00021 { 00022 MPI_Comm comm; KSP solver; PetscErrorCode ierr; 00023 PetscFunctionBegin; 00024 ierr = NumericalProblemGetComm(prob,&comm); CHKERRQ(ierr); 00025 ierr = KSPCreate(comm,&solver); CHKERRQ(ierr); 00026 *(KSP*)ctx = solver; 00027 PetscFunctionReturn(0); 00028 } 00029 00030 #undef __FUNCT__ 00031 #define __FUNCT__ "destroy_solver" 00032 static PetscErrorCode destroy_solver(void *ctx) 00033 { 00034 KSP solver; PetscErrorCode ierr; 00035 PetscFunctionBegin; 00036 solver = (KSP)ctx; 00037 ierr = KSPDestroy(solver); CHKERRQ(ierr); 00038 PetscFunctionReturn(0); 00039 } 00040 00041 #undef __FUNCT__ 00042 #define __FUNCT__ "setup_pc_choices" 00043 static PetscErrorCode setup_pc_choices() 00044 { 00045 PetscErrorCode ierr; 00046 PetscFunctionBegin; 00047 ierr = NewTransformObject("pc",PCNONE,PETSC_NULL); CHKERRQ(ierr); 00048 ierr = NewTransformObject("pc",PCJACOBI,PETSC_NULL); CHKERRQ(ierr); 00049 ierr = NewTransformObject("pc",PCILU,PETSC_NULL); CHKERRQ(ierr); 00050 PetscFunctionReturn(0); 00051 } 00052 00053 #undef __FUNCT__ 00054 #define __FUNCT__ "setup_pc" 00055 static PetscErrorCode setup_pc 00056 (const char *type,int pcv,PetscTruth overwrite, 00057 NumericalProblem inproblem,NumericalProblem *outproblem, 00058 void *gctx,void **ctx,PetscTruth *success) 00059 { 00060 LinearSystem problem = (LinearSystem)inproblem; 00061 KSP solver; PC pc; Mat A; 00062 PetscErrorCode ierr; 00063 00064 PetscFunctionBegin; 00065 *outproblem = inproblem; 00066 00067 /* 00068 * the fun part is setting up the Petsc PC 00069 */ 00070 ierr = LinearSystemGetParts 00071 (problem,&A,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 00072 solver = (KSP)gctx; 00073 ierr = KSPSetOperators(solver,A,A,SAME_NONZERO_PATTERN); CHKERRQ(ierr); 00074 ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr); 00075 ierr = PCSetType(pc,type); CHKERRQ(ierr); 00076 *success = PETSC_TRUE; 00077 00078 PetscFunctionReturn(0); 00079 } 00080 00081 #if 0 00082 #undef __FUNCT__ 00083 #define __FUNCT__ "unset_pc" 00084 static PetscErrorCode unset_pc 00085 (const char *type,PetscTruth overwrite,void *ctx,void *gctx, 00086 NumericalProblem thisproblem,NumericalProblem upproblem, 00087 NumericalSolution old,NumericalSolution new) 00088 { 00089 PetscFunctionBegin; 00090 00091 ierr = LinearSolutionCopy 00092 ((LinearSolution)old,(LinearSolution)new); CHKERRQ(ierr); 00093 00094 PetscFunctionReturn(0); 00095 } 00096 #endif 00097 00098 00099 #undef __FUNCT__ 00100 #define __FUNCT__ "solvelinear" 00101 static PetscErrorCode solvelinear 00102 (NumericalProblem problem,void *dum,NumericalSolution *rsol) 00103 { 00104 LinearSystem linear = (LinearSystem)problem; 00105 KSP solver; Mat A; Vec rhs,sol; 00106 KSPConvergedReason reason; PetscInt it; PetscErrorCode ierr; 00107 PetscFunctionBegin; 00108 ierr = LinearSystemGetParts(linear,&A,NULL,&rhs,NULL,NULL); CHKERRQ(ierr); 00109 ierr = VecDuplicate(rhs,&sol); CHKERRQ(ierr); 00110 ierr = KSPCreate(MPI_COMM_SELF,&solver); CHKERRQ(ierr); 00111 ierr = KSPSetOperators(solver,A,A,DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr); 00112 ierr = KSPSolve(solver,rhs,sol); CHKERRQ(ierr); 00113 ierr = KSPGetConvergedReason(solver,&reason); CHKERRQ(ierr); 00114 if (reason<=0) SETERRQ(1,"This should have converged"); 00115 ierr = KSPGetIterationNumber(solver,&it); CHKERRQ(ierr); 00116 printf("Convergence in %d iterations\n",it); 00117 ierr = KSPDestroy(solver); CHKERRQ(ierr); 00118 *(Vec*)rsol = sol; 00119 PetscFunctionReturn(0); 00120 } 00121 00122 #undef __FUNCT__ 00123 #define __FUNCT__ "destroysolution" 00124 static PetscErrorCode destroysolution(NumericalSolution sol) 00125 { 00126 PetscErrorCode ierr; 00127 PetscFunctionBegin; 00128 ierr = VecDestroy((Vec)sol); CHKERRQ(ierr); 00129 PetscFunctionReturn(0); 00130 } 00131 00132 #undef __FUNCT__ 00133 #define __FUNCT__ "main" 00134 int main(int argc,char **argv) { 00135 NumericalProblem problem; NumericalSolution solution; 00136 Mat A; Vec rhs,sol; PetscErrorCode ierr; 00137 00138 PetscInitialize(&argc,&argv,0,0); 00139 00140 ierr = SysProInitialize(); CHKERRQ(ierr); 00141 00142 ierr = SysProDeclareFunctions 00143 (NULL,NULL,NULL, /* class functions */ 00144 solvelinear,NULL, /* problem solver & problem delete */ 00145 NULL,NULL,destroysolution, /* solution create, copy, delete */ 00146 NULL,NULL,NULL); CHKERRQ(ierr); 00147 00148 /* set up a simple problem */ 00149 { 00150 LinearSystem linearproblem; 00151 int n=100,i; PetscScalar two=2.,mone=-1.,one=1.; 00152 00153 #include "testmat.c" 00154 /* 00155 * now create the linear system to contain them 00156 */ 00157 ierr = CreateLinearSystem(MPI_COMM_SELF,&linearproblem); CHKERRQ(ierr); 00158 ierr = LinearSystemSetParts 00159 (linearproblem,A,NULL,rhs,NULL,NULL); CHKERRQ(ierr); 00160 problem = (NumericalProblem)linearproblem; 00161 } 00162 00163 /* declare a preconditioner loop */ 00164 ierr = DeclarePreprocessor 00165 (/* name */ "pc", 00166 /* one-time creation */ &setup_pc_choices, 00167 /* creation and deletion for specific choices */ NULL,NULL,NULL, 00168 /* create/delete context that will live for the duration of this 00169 preprocessor */ 00170 &create_solver,&destroy_solver, 00171 /* set and unset general/specific ? choice */ 00172 &setup_pc,NULL /*&unset_pc*/ ); CHKERRQ(ierr); 00173 /* and make sure you loop over it */ 00174 ierr = PetscOptionsSetValue("-syspro_pc","exhaustive"); CHKERRQ(ierr); 00175 ierr = PreprocessorsOptionsHandling(); CHKERRQ(ierr); 00176 00177 /* give some feedback on what's happening */ 00178 ierr = SysProDeclareTraceFunction(SysProDefaultTrace); CHKERRQ(ierr); 00179 CHKMEMQ; 00180 00181 /* solve this problem */ 00182 ierr = PreprocessedProblemSolving(problem,&solution); CHKERRQ(ierr); 00183 CHKMEMQ; 00184 00185 /* check for correctness */ 00186 { 00187 Vec sol = (Vec)solution,asol; int mone=-1; PetscReal nrm; 00188 ierr = VecDuplicate(sol,&asol); CHKERRQ(ierr); 00189 ierr = MatMult(A,sol,asol); CHKERRQ(ierr); 00190 ierr = VecAXPY(asol,mone,rhs); CHKERRQ(ierr); 00191 ierr = VecNorm(asol,NORM_2,&nrm); CHKERRQ(ierr); 00192 if (nrm>1.e-6) { 00193 SETERRQ1(1,"Computed wrong solution: norm residual is %e",nrm); 00194 } else printf("Solved within 10e-6 precision\n"); 00195 ierr = VecDestroy(asol); CHKERRQ(ierr); 00196 } 00197 CHKMEMQ; 00198 00199 ierr = MatDestroy(A); CHKERRQ(ierr); 00200 ierr = VecDestroy(rhs); CHKERRQ(ierr); 00201 ierr = VecDestroy(sol); CHKERRQ(ierr); 00202 ierr = VecDestroy((Vec)solution); CHKERRQ(ierr); 00203 ierr = PetscFree(problem); CHKERRQ(ierr); 00204 ierr = SysProFinalize(); CHKERRQ(ierr); 00205 CHKMEMQ; 00206 PetscFinalize(); 00207 PetscFunctionReturn(0); 00208 }