00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
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
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
00172
00173
00174
00175
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
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
00201
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 (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 (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
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 }