System Preprocessors
|
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 }