00001 #include <stdlib.h>
00002 #include "syspro.h"
00003 #include "sysprotransform.h"
00004 #include "sysprolinear.h"
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #undef __FUNCT__
00015 #define __FUNCT__ "create_solver"
00016
00017
00018
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 (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
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 (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,
00144 solvelinear,NULL,
00145 NULL,NULL,destroysolution,
00146 NULL,NULL,NULL); CHKERRQ(ierr);
00147
00148
00149 {
00150 LinearSystem linearproblem;
00151 int n=100,i; PetscScalar two=2.,mone=-1.,one=1.;
00152
00153 #include "testmat.c"
00154
00155
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
00164 ierr = DeclarePreprocessor
00165 ( "pc",
00166 &setup_pc_choices,
00167 NULL,NULL,NULL,
00168
00169
00170 &create_solver,&destroy_solver,
00171
00172 &setup_pc,NULL ); CHKERRQ(ierr);
00173
00174 ierr = PetscOptionsSetValue("-syspro_pc","exhaustive"); CHKERRQ(ierr);
00175 ierr = PreprocessorsOptionsHandling(); CHKERRQ(ierr);
00176
00177
00178 ierr = SysProDeclareTraceFunction(SysProDefaultTrace); CHKERRQ(ierr);
00179 CHKMEMQ;
00180
00181
00182 ierr = PreprocessedProblemSolving(problem,&solution); CHKERRQ(ierr);
00183 CHKMEMQ;
00184
00185
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 }