System Preprocessors
scaling.c
Go to the documentation of this file.
00001 /*! \file scaling.c \ingroup linear */
00002 
00003 /*! \page scaling Scale a linear system
00004 
00005   This preprocessor can perform pointwise left, right, and symmetric scalings
00006   of a linear system by the diagonal of its coefficient matrix.
00007 */
00008 #include <stdlib.h>
00009 #include <stdio.h>
00010 #include "petsc.h"
00011 #include "syspro.h"
00012 #include "sysprotransform.h"
00013 #include "sysprolinear.h"
00014 #include "petscmat.h"
00015 
00016 #define PREPROCESSOR "scaling"
00017 
00018 #undef __FUNCT__
00019 #define __FUNCT__ "set_intelligent_scaling"
00020 static PetscErrorCode set_intelligent_scaling
00021 (NumericalProblem theproblem,const char **type,const char **reason)
00022 {
00023   PetscBool f1,f2,flg; PetscReal r1,r2; PetscErrorCode ierr;
00024 
00025   PetscFunctionBegin;
00026 
00027   ierr = SysProRetrieveQuantity
00028     (theproblem,"variance","row-variability",(void*)&r1,PETSC_NULL,
00029      &f1); CHKERRQ(ierr);
00030   ierr = SysProRetrieveQuantity
00031     (theproblem,"variance","col-variability",(void*)&r2,PETSC_NULL,
00032      &f2); CHKERRQ(ierr);
00033   flg = TRUTH(f1 && f2);
00034   if (flg) {
00035     if (r1>r2+10) {
00036       *type = "right";
00037       *reason = "trying to tame extreme variability within rows";
00038     } else if (r2>r1+10) {
00039       *type = "left";
00040       *reason = "trying to tame extreme variability within columns";
00041     } else {
00042       *type = "symmetric";
00043       *reason = "rows and columns are already pretty equilibrated";
00044     }
00045   } else {
00046     ierr = SysProRetrieveQuantity
00047       (theproblem,"variance","diagonal-average",(void*)&r1,PETSC_NULL,
00048        &f1); CHKERRQ(ierr);
00049     ierr = SysProRetrieveQuantity
00050       (theproblem,"variance","diagonal-variance",(void*)&r2,PETSC_NULL,
00051        &f2); CHKERRQ(ierr);
00052     if (f1 && f2 && r2<r1/20) {
00053       *type = "none";
00054       *reason = "not enough variability to justify scaling";
00055     } else {
00056       *type = "symmetric";
00057       *reason = "symmetric scaling is always better than nothing";
00058     }
00059   }
00060 
00061   PetscFunctionReturn(0);
00062 }
00063 
00064 #undef __FUNCT__
00065 #define __FUNCT__ "scale_system"
00066 static PetscErrorCode scale_system
00067 (const char *type,int nopt,PetscBool overwrite,
00068  NumericalProblem inproblem,NumericalProblem *outproblem,
00069  void *gctx,void **ctx, PetscBool *success)
00070 {
00071   LinearSystem newproblem,problem = (LinearSystem)inproblem;
00072   Vec D; PetscBool flg; PetscErrorCode ierr; 
00073   PetscFunctionBegin;
00074 
00075   ierr = LinearSystemDuplicate(problem,&newproblem); CHKERRQ(ierr);
00076   ierr = LinearSystemCopy(problem,newproblem); CHKERRQ(ierr);
00077 
00078   *success = PETSC_TRUE;
00079   
00080 #if 0
00081   /*
00082    * Scaling leaves the structure intact
00083    */
00084   /*
00085   {
00086     NMD_metadata nmd_in,nmd_out;
00087 
00088     ierr = LinearSystemGetMetadata(problem,&nmd_in); CHKERRQ(ierr);
00089     ierr = LinearSystemGetMetadata(newproblem,&nmd_out); CHKERRQ(ierr);
00090     ierr = SysProTraceMessage("Copying structure categories\n");
00091     ierr = NMDHasCategory(nmd_in,"structure",(int*)&flg); CHKERRQ(ierr);
00092     if (flg) {
00093       ierr = NMDCopyCategory(nmd_in,nmd_out,"structure"); CHKERRQ(ierr);}
00094     ierr = NMDHasCategory(nmd_in,"iprs",(int*)&flg); CHKERRQ(ierr);
00095     if (flg) {
00096       ierr = NMDCopyCategory(nmd_in,nmd_out,"iprs"); CHKERRQ(ierr);}
00097     ierr = NMDHasCategory(nmd_in,"jpl",(int*)&flg); CHKERRQ(ierr);
00098     if (flg) {
00099       ierr = NMDCopyCategory(nmd_in,nmd_out,"jpl"); CHKERRQ(ierr);}
00100     ierr = NMDHasCategory(nmd_in,"icmk",(int*)&flg); CHKERRQ(ierr);
00101     if (flg) {
00102       ierr = NMDCopyCategory(nmd_in,nmd_out,"icmk"); CHKERRQ(ierr);}
00103   }
00104   */
00105 #endif
00106 
00107   ierr = PetscStrcmp(type,"none",&flg); CHKERRQ(ierr);
00108   if (flg)
00109     goto done;
00110 
00111   /*
00112    * construct the scaling
00113    */
00114   {
00115     Mat A; int local_size; MPI_Comm comm;
00116 
00117     ierr = LinearSystemGetParts
00118       (problem,&A,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00119     ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00120     ierr = MatGetLocalSize(A,&local_size,&ierr); CHKERRQ(ierr);
00121     ierr = VecCreateMPI(comm,local_size,PETSC_DECIDE,&D); CHKERRQ(ierr);
00122     ierr = MatGetDiagonal(A,D); CHKERRQ(ierr);
00123     { /* prevent problems with zero diagonal elements */
00124       PetscScalar *a; int local_size,i;
00125       ierr = VecGetArray(D,&a); CHKERRQ(ierr);
00126       ierr = MatGetLocalSize(A,&local_size,&i); CHKERRQ(ierr);
00127       for (i=0; i<local_size; i++)
00128         if (a[i]==0.) a[i]=1.;
00129       ierr = VecRestoreArray(D,&a); CHKERRQ(ierr);
00130     }
00131     ierr = VecReciprocal(D); CHKERRQ(ierr);
00132     *(Vec*)ctx = D;
00133   }
00134 
00135   /*
00136    * do the actual scaling
00137    */
00138   {
00139     Mat A,B; Vec rhs,sol,init;
00140     ierr = LinearSystemGetParts
00141       (newproblem,&A,&B,&rhs,&sol,&init); CHKERRQ(ierr);
00142 
00143     ierr = PetscStrcmp(type,"symmetric",&flg); CHKERRQ(ierr);
00144     if (flg) {
00145       ierr = VecSqrtAbs(D); CHKERRQ(ierr);
00146       ierr = MatDiagonalScale(A,D,D); CHKERRQ(ierr);
00147       if (B && B!=A) {
00148         ierr = MatDiagonalScale(B,D,D); CHKERRQ(ierr);
00149       }
00150       ierr = VecPointwiseMult(rhs,D,rhs); CHKERRQ(ierr);
00151       ierr = VecPointwiseDivide(sol,sol,D); CHKERRQ(ierr);
00152       if (init) {
00153         ierr = VecPointwiseDivide(init,init,D); CHKERRQ(ierr);}
00154       goto done;
00155     }
00156     
00157     ierr = PetscStrcmp(type,"left",&flg); CHKERRQ(ierr);
00158     if (flg) {
00159       ierr = MatDiagonalScale(A,D,PETSC_NULL); CHKERRQ(ierr);
00160       if (B && B!=A) {
00161         ierr = MatDiagonalScale(B,D,PETSC_NULL); CHKERRQ(ierr);
00162       }
00163       ierr = VecPointwiseMult(rhs,D,rhs); CHKERRQ(ierr);
00164       goto done;
00165     }
00166     
00167     ierr = PetscStrcmp(type,"right",&flg); CHKERRQ(ierr);
00168     if (flg) {
00169       ierr = MatDiagonalScale(A,PETSC_NULL,D); CHKERRQ(ierr);
00170       if (B && B!=A) {
00171         ierr = MatDiagonalScale(B,PETSC_NULL,D); CHKERRQ(ierr);
00172       }
00173       ierr = VecPointwiseDivide(sol,sol,D); CHKERRQ(ierr);
00174       if (init) {
00175         ierr = VecPointwiseDivide(init,init,D); CHKERRQ(ierr);}
00176       goto done;
00177     }
00178     
00179     SETERRQ1(PETSC_COMM_SELF,1,"scaling type not implemented: %s",type);
00180   }
00181 
00182  done:
00183   *outproblem = (NumericalProblem)newproblem;
00184   PetscFunctionReturn(0);
00185 }
00186 
00187 #undef __FUNCT__
00188 #define __FUNCT__ "unscale_system"
00189 static PetscErrorCode unscale_system
00190 (const char *scaling_type,PetscBool overwrite,void *gctx,void *ctx,
00191  NumericalProblem problem,NumericalProblem nextproblem,
00192  NumericalSolution scaled,NumericalSolution unscaled)
00193 {
00194   LinearSolution old = (LinearSolution)scaled,lnew = (LinearSolution)unscaled;
00195   Vec D = (Vec)ctx;
00196   Vec X,Xscaled; PetscBool fnone,f1,f2; PetscErrorCode ierr;
00197 
00198   PetscFunctionBegin;
00199   ierr = LinearSolutionCopyStats(old,lnew); CHKERRQ(ierr);
00200   ierr = LinearSolutionGetVector(old,&Xscaled); CHKERRQ(ierr);
00201   ierr = LinearSolutionGetVector(lnew,&X); CHKERRQ(ierr);
00202   ierr = VecCopy(Xscaled,X); CHKERRQ(ierr);
00203 
00204   ierr = PetscStrcmp(scaling_type,"none",&f1); CHKERRQ(ierr);
00205   fnone = f1;
00206   ierr = PetscStrcmp(scaling_type,"left",&f2); CHKERRQ(ierr);
00207   if (f1 || f2) {
00208     ierr = VecCopy(Xscaled,X); CHKERRQ(ierr);
00209   } else {
00210     ierr = PetscStrcmp(scaling_type,"right",&f1); CHKERRQ(ierr);
00211     ierr = PetscStrcmp(scaling_type,"symmetric",&f2); CHKERRQ(ierr);
00212     if (f1 || f2) {
00213       ierr = VecPointwiseMult(X,Xscaled,D); CHKERRQ(ierr);
00214     } else {
00215       SETERRQ1(PETSC_COMM_SELF,1,"Unknown scaling <%s>",scaling_type); 
00216     }
00217   }
00218   if (!fnone) {
00219     ierr = VecDestroy(&D); CHKERRQ(ierr);}
00220 
00221   /*ierr = LinearSolutionDelete((LinearSolution)scaled); CHKERRQ(ierr);*/
00222   //ierr = DeleteLinearSystem((LinearSystem)problem); CHKERRQ(ierr);
00223 
00224   PetscFunctionReturn(0);
00225 }
00226 
00227 #undef __FUNCT__
00228 #define __FUNCT__ "setup_scaling_choices"
00229 /*! This routine is called by DeclarePreprocessor() */
00230 static PetscErrorCode setup_scaling_choices()
00231 {
00232   SalsaTransform tf; SalsaTransformObject cur; PetscErrorCode ierr;
00233 
00234   PetscFunctionBegin;
00235   ierr = TransformGetByName(PREPROCESSOR,&tf); CHKERRQ(ierr);
00236   ierr = SysProDefineIntAnnotation(PREPROCESSOR,"symmetric"); CHKERRQ(ierr);
00237 
00238   ierr = NewTransformObject(PREPROCESSOR,"none",&cur); CHKERRQ(ierr);
00239   ierr = TransformObjectSetExplanation(cur,"no scaling"); CHKERRQ(ierr);
00240   ierr = TransformObjectIntAnnotate(cur,"symmetric",1); CHKERRQ(ierr);
00241 
00242   ierr = NewTransformObject(PREPROCESSOR,"symmetric",&cur); CHKERRQ(ierr);
00243   ierr = TransformObjectSetExplanation(cur,"symmetric scaling"); CHKERRQ(ierr);
00244   ierr = TransformObjectIntAnnotate(cur,"symmetric",1); CHKERRQ(ierr);
00245 
00246   ierr = NewTransformObject(PREPROCESSOR,"left",&cur); CHKERRQ(ierr);
00247   ierr = TransformObjectSetExplanation(cur,"left scaling"); CHKERRQ(ierr);
00248 
00249   ierr = NewTransformObject(PREPROCESSOR,"right",&cur); CHKERRQ(ierr);
00250   ierr = TransformObjectSetExplanation(cur,"right scaling"); CHKERRQ(ierr);
00251 
00252   PetscFunctionReturn(0);
00253 }
00254 
00255 #undef __FUNCT__
00256 #define __FUNCT__ "specific_scaling_choices"
00257 /*!
00258   This is the 'specific setup' phase of the scaling preprocessor.
00259   See \ref modes for details.
00260 
00261   This routine eliminates unsymmetric scalings if we are dealing with a 
00262   symmetric system.
00263 */
00264 static PetscErrorCode specific_scaling_choices
00265 (NumericalProblem theproblem,SalsaTransform scaling)
00266 {
00267   int isymm; PetscBool flg; PetscErrorCode ierr;
00268   PetscFunctionBegin;
00269 
00270   ierr = SysProRetrieveQuantity
00271     (theproblem,"structure","symmetry",(void*)&isymm,PETSC_NULL,
00272      &flg); CHKERRQ(ierr);
00273   if (flg && isymm>0) {
00274     int i,n; SalsaTransformObject *objs;
00275     ierr = TransformGetObjects(scaling,&n,&objs); CHKERRQ(ierr);
00276     for (i=0; i<n; i++) {
00277       SalsaTransformObject cur = objs[i]; 
00278       int v; PetscBool f;
00279       ierr = TransformObjectGetIntAnnotation
00280         (cur,"symmetric",&v,&f); CHKERRQ(ierr);
00281       if (f && !v) {
00282         ierr = TransformObjectMark(cur); CHKERRQ(ierr);
00283       }
00284     }
00285   }
00286 
00287   PetscFunctionReturn(0);
00288 }
00289 
00290 #undef __FUNCT__
00291 #define __FUNCT__ "DeclareScalingPreprocessor"
00292 PetscErrorCode DeclareScalingPreprocessor(void)
00293 {
00294   PetscErrorCode ierr;
00295   PetscFunctionBegin;
00296   ierr = DeclarePreprocessor
00297     (PREPROCESSOR,
00298      &setup_scaling_choices, /* declare all scalings */
00299      &specific_scaling_choices, /* metadata inspection to set intell choice */
00300      0,0, /* local and global unset */
00301      0,0, /* context create and destroy */
00302      &scale_system, /* start function */
00303      &unscale_system /* end function */); CHKERRQ(ierr);
00304   ierr = DeclarePreprocessorIntelligentChoice
00305     ("scaling",&set_intelligent_scaling); CHKERRQ(ierr);
00306   ierr = PreprocessorSetPreservedCategories
00307     (PREPROCESSOR,"laplace,structure,iprs,icmk"); CHKERRQ(ierr);
00308   PetscFunctionReturn(0);
00309 }