System Preprocessors
u13.c
Go to the documentation of this file.
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 }