System Preprocessors
|
00001 /*! \file approximating.c \ingroup linear */ 00002 00003 /*! \page approximation Approximate the coefficient matrix 00004 00005 A preconditioner need not be derived from the coefficient matrix. 00006 For instance, in the case of a higher order finite element matrix, 00007 incomplete factorization preconditioner are better derived from a 00008 linear element discretization of the same problem, since this matrix 00009 will be an M-matrix. 00010 00011 This preprocessor can perform the following approximations: 00012 - symmetric: take the symmetric part of the coefficient matrix 00013 - gustafsson: apply the Gustafsson modified element matrix transformation 00014 (see reference [GUS] below). 00015 00016 \verbatim 00017 [GUS] 00018 @article{Gu:modified_element, 00019 author = {Ivar Gustafsson}, 00020 title = {An Incomplete Factorization Preconditioning Method 00021 based on Modification of Element Matrices}, 00022 journal = {BIT}, 00023 year = {1996}, 00024 volume = {36}, 00025 pages = {86--100} 00026 } 00027 \endverbatim 00028 */ 00029 #include <stdlib.h> 00030 #include <stdio.h> 00031 #include "petsc.h" 00032 #include "petscis.h" 00033 #include "petsc.h" 00034 #include "syspro.h" 00035 #include "sysprotransform.h" 00036 #include "sysprolinear.h" 00037 #include "linear_impl.h" 00038 00039 #define PREPROCESSOR "approximation" 00040 00041 #undef __FUNCT__ 00042 #define __FUNCT__ "MatSymmetricPart" 00043 static PetscErrorCode MatSymmetricPart 00044 (NumericalProblem inproblem,NumericalProblem outproblem) 00045 { 00046 LinearSystem problem = (LinearSystem)inproblem, 00047 newproblem = (LinearSystem)outproblem; 00048 Mat A,At,B; int first,last,nun,i; PetscScalar half=0.5,zero=0.,one=1.; 00049 PetscReal avg; PetscBool flg; PetscErrorCode ierr; 00050 00051 PetscFunctionBegin; 00052 ierr = LinearSystemGetParts 00053 (problem,&A,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 00054 ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr); 00055 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&At); CHKERRQ(ierr); 00056 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B); CHKERRQ(ierr); 00057 ierr = SysProComputeQuantity 00058 (inproblem,"structure","n-struct-unsymm",(void*)&nun,PETSC_NULL, 00059 &flg); CHKERRQ(ierr); 00060 if (flg && nun==0) { 00061 ierr = MatAXPY(B,one,At,SAME_NONZERO_PATTERN); CHKERRQ(ierr); 00062 } else { 00063 ierr = MatAXPY(B,one,At,DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr); 00064 } 00065 ierr = MatScale(B,half); CHKERRQ(ierr); 00066 for (i=first; i<last; i++) { 00067 ierr = MatSetValues(B,1,&i,1,&i,&zero,ADD_VALUES); CHKERRQ(ierr); 00068 } 00069 ierr = MatDestroy(&At); CHKERRQ(ierr); 00070 ierr = SysProRetrieveQuantity 00071 (inproblem,"variance","diagonal-average",(void*)&avg,PETSC_NULL, 00072 &flg); CHKERRQ(ierr); 00073 if (flg) { 00074 int first,last,row; 00075 avg *= 1.e-12; 00076 ierr = MatGetOwnershipRange(B,&first,&last); CHKERRQ(ierr); 00077 for (row=first; row<last; row++) { 00078 ierr = MatSetValues 00079 (B,1,&row,1,&row,&avg,ADD_VALUES); CHKERRQ(ierr); 00080 } 00081 } 00082 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 00083 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 00084 ierr = LinearSystemSetParts 00085 (newproblem,PETSC_NULL,B,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 00086 PetscFunctionReturn(0); 00087 } 00088 00089 #undef __FUNCT__ 00090 #define __FUNCT__ "MatGustafssonMod" 00091 static PetscErrorCode MatGustafssonMod 00092 (NumericalProblem inproblem,NumericalProblem outproblem) 00093 { 00094 LinearSystem problem = (LinearSystem)inproblem, 00095 newproblem = (LinearSystem)outproblem; 00096 Mat A,B; PetscScalar *a; int first,last,row,ierr; 00097 PetscFunctionBegin; 00098 00099 ierr = LinearSystemGetParts 00100 (problem,&A,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 00101 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B); CHKERRQ(ierr); 00102 ierr = MatGetOwnershipRange(B,&first,&last); CHKERRQ(ierr); 00103 for (row=first; row<last; row++) { 00104 int icol,ncols; const int *cols; const PetscScalar *vals; 00105 ierr = MatGetRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr); 00106 for (icol=0; icol<ncols; icol++) { 00107 if (vals[icol]<0) { 00108 PetscScalar v = -vals[icol]; 00109 ierr = MatSetValues 00110 (B,1,&row,1,cols+icol,&v,ADD_VALUES); CHKERRQ(ierr); 00111 ierr = MatSetValues 00112 (B,1,&row,1,&row,&v,ADD_VALUES); CHKERRQ(ierr); 00113 } 00114 } 00115 ierr = MatRestoreRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr); 00116 } 00117 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 00118 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 00119 00120 { 00121 /* 00122 * Enforce non-negative diagonal 00123 */ 00124 Vec V; 00125 ierr = VecCreateSeq(MPI_COMM_SELF,last-first,&V); CHKERRQ(ierr); 00126 ierr = MatGetDiagonal(B,V); CHKERRQ(ierr); 00127 ierr = VecGetArray(V,&a); CHKERRQ(ierr); 00128 for (row=first; row<last; row++) { 00129 if (a[row-first]<0) { 00130 PetscScalar v=0; 00131 ierr = MatSetValues(B,1,&row,1,&row,&v,INSERT_VALUES); CHKERRQ(ierr); 00132 } 00133 } 00134 ierr = VecRestoreArray(V,&a); CHKERRQ(ierr); 00135 ierr = VecDestroy(&V); CHKERRQ(ierr); 00136 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 00137 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 00138 } 00139 ierr = LinearSystemSetParts 00140 (newproblem,PETSC_NULL,B,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 00141 PetscFunctionReturn(0); 00142 } 00143 00144 #undef __FUNCT__ 00145 #define __FUNCT__ "setup_approximation_choices" 00146 static PetscErrorCode setup_approximation_choices() 00147 { 00148 SalsaTransform tf; SalsaTransformObject cur; PetscErrorCode ierr; 00149 00150 PetscFunctionBegin; 00151 ierr = TransformGetByName(PREPROCESSOR,&tf); CHKERRQ(ierr); 00152 00153 ierr = NewTransformObject(PREPROCESSOR,"none",&cur); CHKERRQ(ierr); 00154 ierr = TransformObjectSetExplanation 00155 (cur,"use coefficient matrix"); CHKERRQ(ierr); 00156 00157 ierr = NewTransformObject(PREPROCESSOR,"symmetric",&cur); CHKERRQ(ierr); 00158 ierr = TransformObjectSetExplanation 00159 (cur,"symmetric part"); CHKERRQ(ierr); 00160 00161 ierr = NewTransformObject(PREPROCESSOR,"gustafsson",&cur); CHKERRQ(ierr); 00162 ierr = TransformObjectSetExplanation 00163 (cur,"gustafsson M-matrix"); CHKERRQ(ierr); 00164 00165 PetscFunctionReturn(0); 00166 } 00167 00168 #undef __FUNCT__ 00169 #define __FUNCT__ "specific_approximation_choices" 00170 /*! 00171 This is the 'specific setup' phase of the approximation preprocessor. 00172 See \ref modes for details. 00173 00174 This routine eliminates the Gustafsson approximation for diagonally dominant 00175 systems, and the symmetric for symmetric systems. 00176 */ 00177 static PetscErrorCode specific_approximation_choices 00178 (NumericalProblem inproblem,SalsaTransform transform) 00179 { 00180 LinearSystem problem = (LinearSystem)inproblem; Mat A; 00181 PetscReal dd,anorm; PetscBool flg; PetscErrorCode ierr; 00182 PetscFunctionBegin; 00183 ierr = LinearSystemGetParts 00184 (problem,&A,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 00185 00186 /* 00187 * No Gustafsson approximation for diagonally dominant matrices 00188 */ 00189 ierr = SysProRetrieveQuantity 00190 (inproblem,"simple","diagonal-dominance",(void*)&dd,PETSC_NULL, 00191 &flg); CHKERRQ(ierr); 00192 if (flg && dd>=0.) { 00193 SalsaTransformObject tf; 00194 ierr = TransformObjectGetByName 00195 (PREPROCESSOR,"gustafsson",&tf); CHKERRQ(ierr); 00196 ierr = TransformObjectMark(tf); CHKERRQ(ierr); 00197 } 00198 00199 /* 00200 * symmetric approximations of a symmetric matrices 00201 * are a silly idea 00202 */ 00203 ierr = SysProRetrieveQuantity 00204 (inproblem,"simple","symmetry-anorm",(void*)&anorm,PETSC_NULL, 00205 &flg); CHKERRQ(ierr); 00206 if (flg && anorm==0.) { 00207 SalsaTransformObject tf; 00208 ierr = TransformObjectGetByName 00209 (PREPROCESSOR,"symmetric",&tf); CHKERRQ(ierr); 00210 ierr = TransformObjectMark(tf); CHKERRQ(ierr); 00211 } 00212 00213 00214 PetscFunctionReturn(0); 00215 } 00216 00217 #undef __FUNCT__ 00218 #define __FUNCT__ "approximate_system" 00219 static PetscErrorCode approximate_system 00220 (const char *type,int nopt,PetscBool overwrite, 00221 NumericalProblem inproblem,NumericalProblem *outproblem, 00222 void *gctx,void **ctx, PetscBool *success) 00223 { 00224 LinearSystem problem = (LinearSystem)inproblem, newproblem; 00225 PetscBool flg; PetscErrorCode ierr; 00226 PetscFunctionBegin; 00227 00228 { 00229 Mat A,B; 00230 ierr = LinearSystemGetParts 00231 (problem,&A,&B,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 00232 if (A!=B) SETERRQ(PETSC_COMM_SELF,1,"A and B should be the same before approximation"); 00233 } 00234 ierr = LinearSystemDuplicatePointers(problem,&newproblem); CHKERRQ(ierr); 00235 00236 *success = PETSC_TRUE; 00237 00238 ierr = PetscStrcmp("none",type,&flg); CHKERRQ(ierr); 00239 if (flg) { 00240 goto done; 00241 } 00242 00243 ierr = PetscStrcmp("symmetric",type,&flg); CHKERRQ(ierr); 00244 if (flg) { 00245 ierr = MatSymmetricPart(inproblem,(NumericalProblem)newproblem); CHKERRQ(ierr); 00246 goto done; 00247 } 00248 00249 ierr = PetscStrcmp("gustafsson",type,&flg); CHKERRQ(ierr); 00250 if (flg) { 00251 ierr = MatGustafssonMod(inproblem,(NumericalProblem)newproblem); CHKERRQ(ierr); 00252 goto done; 00253 } 00254 00255 SETERRQ1(PETSC_COMM_SELF,1,"Unknown matrix approximation <%s>",type); 00256 00257 done: 00258 *outproblem = (NumericalProblem)newproblem; 00259 PetscFunctionReturn(0); 00260 } 00261 00262 #undef __FUNCT__ 00263 #define __FUNCT__ "unapproximate_system" 00264 static PetscErrorCode unapproximate_system 00265 (const char *type,PetscBool overwrite,void *gctx,void *ctx, 00266 NumericalProblem problem,NumericalProblem nextproblem, 00267 NumericalSolution before,NumericalSolution after) 00268 { 00269 PetscErrorCode ierr; 00270 PetscFunctionBegin; 00271 ierr = LinearSolutionCopy 00272 ((LinearSolution)before,(LinearSolution)after); CHKERRQ(ierr); 00273 /*ierr = LinearSolutionDelete((LinearSolution)before); CHKERRQ(ierr);*/ 00274 //ierr = DeleteLinearSystem((LinearSystem)problem); CHKERRQ(ierr); 00275 PetscFunctionReturn(0); 00276 } 00277 00278 #undef __FUNCT__ 00279 #define __FUNCT__ "DeclareApproximationPreprocessor" 00280 PetscErrorCode DeclareApproximationPreprocessor(void) 00281 { 00282 PetscErrorCode ierr; 00283 PetscFunctionBegin; 00284 ierr = DeclarePreprocessor 00285 ("approximation", 00286 &setup_approximation_choices,&specific_approximation_choices,0,0, 00287 0,0, 00288 &approximate_system,&unapproximate_system); CHKERRQ(ierr); 00289 ierr = PreprocessorSetPreservedCategories 00290 (PREPROCESSOR,"laplace"); CHKERRQ(ierr); 00291 PetscFunctionReturn(0); 00292 }