00001 #include <stdlib.h>
00002 #include "syspro.h"
00003 #include "sysprolinear.h"
00004
00005
00006
00007
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
00037 {
00038 LinearSystem linearproblem;
00039 int n=10,i; PetscScalar two=2.,mone=-1.,one=1.;
00040
00041 #include "testmat.c"
00042
00043
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
00052 ierr = PreprocessedProblemSolving(problem,&solution); CHKERRQ(ierr);
00053
00054
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 }