00001
00002
00003
00004
00005
00006
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,char **type,char **reason)
00022 {
00023 PetscTruth 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 (char *type,int nopt,PetscTruth overwrite,
00068 NumericalProblem inproblem,NumericalProblem *outproblem,
00069 void *gctx,void **ctx, PetscTruth *success)
00070 {
00071 LinearSystem newproblem,problem = (LinearSystem)inproblem;
00072 Vec D; PetscTruth 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
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 #endif
00106
00107 ierr = PetscStrcmp(type,"none",&flg); CHKERRQ(ierr);
00108 if (flg)
00109 goto done;
00110
00111
00112
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 {
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
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 = VecSqrt(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(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 (char *scaling_type,PetscTruth 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; PetscTruth 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(1,"Unknown scaling <%s>",scaling_type);
00216 }
00217 }
00218 if (!fnone) {
00219 ierr = VecDestroy(D); CHKERRQ(ierr);}
00220
00221
00222 ierr = DeleteLinearSystem((LinearSystem)problem); CHKERRQ(ierr);
00223
00224 PetscFunctionReturn(0);
00225 }
00226
00227 #undef __FUNCT__
00228 #define __FUNCT__ "setup_scaling_choices"
00229
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
00259
00260
00261
00262
00263
00264 static PetscErrorCode specific_scaling_choices
00265 (NumericalProblem theproblem,SalsaTransform scaling)
00266 {
00267 int isymm; PetscTruth 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; PetscTruth 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,
00299 &specific_scaling_choices,
00300 0,0,
00301 0,0,
00302 &scale_system,
00303 &unscale_system ); 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 }