System Preprocessors
u12.c
Go to the documentation of this file.
00001 #include <stdlib.h>
00002 #include "syspro.h"
00003 #include "sysprolinear.h"
00004 
00005 /* test for the syspro linear package
00006    test for problem solving without preprocessors:
00007    this uses PETSc's KSP for actually solving a linear system
00008 */
00009 
00010 #undef __FUNCT__
00011 #define __FUNCT__ "solvelinear"
00012 static PetscErrorCode solvelinear
00013 (NumericalProblem problem,void *dum,NumericalSolution *rsol)
00014 {
00015   LinearSystem linear = (LinearSystem)problem;
00016   KSP solver; Mat A; Vec rhs,sol; PetscErrorCode ierr;
00017   PetscFunctionBegin;
00018   ierr = LinearSystemGetParts(linear,&A,NULL,&rhs,NULL,NULL); CHKERRQ(ierr);
00019   ierr = VecDuplicate(rhs,&sol); CHKERRQ(ierr);
00020   ierr = KSPCreate(MPI_COMM_SELF,&solver); CHKERRQ(ierr);
00021   ierr = KSPSetOperators(solver,A,A,(MatStructure)0); CHKERRQ(ierr);
00022   ierr = KSPSolve(solver,rhs,sol); CHKERRQ(ierr);
00023   ierr = KSPDestroy(solver); CHKERRQ(ierr);
00024   *(Vec*)rsol = sol;
00025   PetscFunctionReturn(0);
00026 }
00027 
00028 int main(int argc,char **argv) {
00029   NumericalProblem problem; NumericalSolution solution;
00030   Mat A; Vec rhs,sol; PetscErrorCode ierr;
00031   PetscInitialize(&argc,&argv,0,0);
00032   ierr = SysProInitialize(); CHKERRQ(ierr);
00033   ierr = SysProDeclareFunctions
00034     (NULL,NULL,NULL,solvelinear,NULL,NULL,NULL,NULL,NULL,NULL,NULL); CHKERRQ(ierr);
00035   
00036   /* set up a simple problem */
00037   {
00038     LinearSystem linearproblem;
00039     int n=10,i; PetscScalar two=2.,mone=-1.,one=1.;
00040 
00041 #include "testmat.c"
00042     /*
00043      * now create the linear system to contain them
00044      */
00045     ierr = CreateLinearSystem(MPI_COMM_SELF,&linearproblem); CHKERRQ(ierr);
00046     ierr = LinearSystemSetParts
00047       (linearproblem,A,NULL,rhs,NULL,NULL); CHKERRQ(ierr);
00048     problem = (NumericalProblem)linearproblem;
00049   }
00050 
00051   /* solve this problem */
00052   ierr = PreprocessedProblemSolving(problem,&solution); CHKERRQ(ierr);
00053 
00054   /* check for correctness */
00055   {
00056     Vec sol = (Vec)solution,asol; int mone=-1; PetscReal nrm;
00057     ierr = VecDuplicate(sol,&asol); CHKERRQ(ierr);
00058     ierr = MatMult(A,sol,asol); CHKERRQ(ierr);
00059     ierr = VecAXPY(asol,mone,rhs); CHKERRQ(ierr);
00060     ierr = VecNorm(asol,NORM_2,&nrm); CHKERRQ(ierr);
00061     if (nrm>1.e-6) {
00062       SETERRQ1(1,"Computed wrong solution: norm residual is %e",nrm);
00063     } else {
00064       printf("Solution computed to accuracy\n");
00065     }
00066     ierr = VecDestroy(asol); CHKERRQ(ierr);
00067     ierr = VecDestroy(sol); CHKERRQ(ierr);
00068   }
00069 
00070   ierr = MatDestroy(A); CHKERRQ(ierr);
00071   ierr = VecDestroy(rhs); CHKERRQ(ierr);
00072   ierr = VecDestroy(sol); CHKERRQ(ierr);
00073   ierr = PetscFree(problem); CHKERRQ(ierr);
00074   ierr = SysProFinalize(); CHKERRQ(ierr);
00075   PetscFinalize();
00076   PetscFunctionReturn(0);
00077 }