System Preprocessors
|
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(PETSC_COMM_SELF,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 }