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 #undef __FUNCT__
00082 #define __FUNCT__ "unset_pc"
00083 static PetscErrorCode unset_pc
00084 (char *type,PetscTruth overwrite,void *ctx,void *gctx,
00085 NumericalProblem thisproblem,NumericalProblem upproblem,
00086 NumericalSolution old,NumericalSolution nnew)
00087 {
00088 #if 0
00089 LinearSystem problem = (LinearSystem)thisproblem;
00090 Mat A;
00091 PetscErrorCode ierr;
00092 #endif
00093 PetscFunctionBegin;
00094
00095 #if 0
00096 ierr = LinearSolutionCopy
00097 ((LinearSolution)old,(LinearSolution)nnew); CHKERRQ(ierr);
00098 #endif
00099
00100 PetscFunctionReturn(0);
00101 }
00102
00103
00104 #undef __FUNCT__
00105 #define __FUNCT__ "solvelinear"
00106 static PetscErrorCode solvelinear
00107 (NumericalProblem problem,void *dum,NumericalSolution *rsol)
00108 {
00109 LinearSystem linear = (LinearSystem)problem;
00110 KSP solver; Mat A; Vec rhs,sol;
00111 KSPConvergedReason reason; PetscInt it; PetscErrorCode ierr;
00112 PetscFunctionBegin;
00113 ierr = LinearSystemGetParts(linear,&A,NULL,&rhs,NULL,NULL); CHKERRQ(ierr);
00114 ierr = VecDuplicate(rhs,&sol); CHKERRQ(ierr);
00115 ierr = KSPCreate(MPI_COMM_SELF,&solver); CHKERRQ(ierr);
00116 ierr = KSPSetOperators(solver,A,A,DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
00117 ierr = KSPSolve(solver,rhs,sol); CHKERRQ(ierr);
00118 ierr = KSPGetConvergedReason(solver,&reason); CHKERRQ(ierr);
00119 if (reason<=0) SETERRQ(1,"This should have converged");
00120 ierr = KSPGetIterationNumber(solver,&it); CHKERRQ(ierr);
00121 printf("Convergence in %d iterations\n",it);
00122 ierr = KSPDestroy(solver); CHKERRQ(ierr);
00123 ierr = LinearCreateNumericalSolution(NULL,rsol); CHKERRQ(ierr);
00124 ierr = LinearSolutionSetVector((LinearSolution)*rsol,sol); CHKERRQ(ierr);
00125 PetscFunctionReturn(0);
00126 }
00127
00128 int main(int argc,char **argv) {
00129 NumericalProblem problem; NumericalSolution solution;
00130 Mat A; Vec rhs,sol; PetscErrorCode ierr;
00131 PetscInitialize(&argc,&argv,0,0);
00132 ierr = SysProInitialize(); CHKERRQ(ierr);
00133
00134
00135 {
00136 LinearSystem linearproblem;
00137 int n=100,i; PetscScalar two=2.,mone=-1.,one=1.;
00138
00139 #include "testmat.c"
00140
00141
00142
00143 ierr = CreateLinearSystem(MPI_COMM_SELF,&linearproblem); CHKERRQ(ierr);
00144 ierr = LinearSystemSetParts
00145 (linearproblem,A,NULL,rhs,sol,NULL); CHKERRQ(ierr);
00146 problem = (NumericalProblem)linearproblem;
00147 }
00148
00149
00150 ierr = DeclareScalingPreprocessor(); CHKERRQ(ierr);
00151 ierr = PetscOptionsSetValue("-syspro_scaling","exhaustive"); CHKERRQ(ierr);
00152
00153
00154 ierr = DeclarePreprocessor
00155 ( "pc",
00156 &setup_pc_choices,
00157 NULL,NULL,NULL,
00158
00159
00160 &create_solver,&destroy_solver,
00161
00162 &setup_pc,NULL ); CHKERRQ(ierr);
00163
00164 ierr = PetscOptionsSetValue("-syspro_pc","none"); CHKERRQ(ierr);
00165
00166
00167 ierr = SysProDeclareFunctions
00168 (NULL,NULL,NULL,
00169 solvelinear,NULL,
00170 &LinearCreateNumericalSolution,NULL,LinearDeleteNumericalSolution,
00171 NULL,NULL,NULL
00172 ); CHKERRQ(ierr);
00173
00174 ierr = PreprocessorsOptionsHandling(); CHKERRQ(ierr);
00175
00176
00177 ierr = SysProDeclareTraceFunction(SysProDefaultTrace); CHKERRQ(ierr);
00178
00179
00180 ierr = PreprocessedProblemSolving(problem,&solution); CHKERRQ(ierr);
00181
00182
00183 {
00184 Vec sol,asol; int mone=-1; PetscReal nrm;
00185 ierr = LinearSolutionGetVector
00186 ((LinearSolution)solution,&sol); CHKERRQ(ierr);
00187 ierr = VecDuplicate(sol,&asol); CHKERRQ(ierr);
00188 ierr = MatMult(A,sol,asol); CHKERRQ(ierr);
00189 ierr = VecAXPY(asol,mone,rhs); CHKERRQ(ierr);
00190 ierr = VecNorm(asol,NORM_2,&nrm); CHKERRQ(ierr);
00191 if (nrm>1.e-6)
00192 SETERRQ1(1,"Computed wrong solution: norm residual is %e",nrm);
00193 ierr = VecDestroy(asol); CHKERRQ(ierr);
00194 }
00195 ierr = LinearDeleteNumericalSolution(solution); CHKERRQ(ierr);
00196 ierr = MatDestroy(A); CHKERRQ(ierr);
00197 ierr = VecDestroy(rhs); CHKERRQ(ierr);
00198 ierr = VecDestroy(sol); CHKERRQ(ierr);
00199 ierr = PetscFree(problem); CHKERRQ(ierr);
00200 ierr = SysProFinalize(); CHKERRQ(ierr);
00201 PetscFinalize();
00202 PetscFunctionReturn(0);
00203 }