System Preprocessors
approximating.c
Go to the documentation of this file.
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; PetscTruth 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; PetscTruth 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,PetscTruth overwrite,
00221  NumericalProblem inproblem,NumericalProblem *outproblem,
00222  void *gctx,void **ctx, PetscTruth *success)
00223 {
00224   LinearSystem problem = (LinearSystem)inproblem, newproblem;
00225   PetscTruth 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(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(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,PetscTruth 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 }