System Preprocessors
|
00001 /*! \file singleton.c \ingroup linear */ 00002 00003 /*! \page singleton Eliminate fully determined (singleton) rows from a matrix 00004 00005 In a linear system, any variable whose matrix row has only a single 00006 element is independent of the rest of the system, in the sense that the 00007 other variables can be solved independently of it. The singleton 00008 preprocessor eliminates such rows, and performs the backsubstitution. 00009 */ 00010 #include <stdlib.h> 00011 #include <stdio.h> 00012 #include "petsc.h" 00013 #include "syspro.h" 00014 #include "sysprotransform.h" 00015 #include "sysprolinear.h" 00016 #include "linear_impl.h" 00017 00018 typedef struct {int n,t; VecScatter extractor;} singleton_struct; 00019 #define PREPROCESSOR "singleton" 00020 00021 #undef __FUNCT__ 00022 #define __FUNCT__ "eliminate_singletons" 00023 static PetscErrorCode eliminate_singletons 00024 (const char *type,int nopt,PetscTruth overwrite, 00025 NumericalProblem inproblem,NumericalProblem *outproblem, 00026 void *gctx,void **ctx,PetscTruth *success) 00027 { 00028 LinearSystem problem = (LinearSystem)inproblem, newproblem; 00029 PetscTruth flg; PetscErrorCode ierr; 00030 00031 PetscFunctionBegin; 00032 00033 *success = PETSC_TRUE; *(Vec*)ctx = NULL; 00034 ierr = LinearSystemDuplicatePointers(problem,&newproblem); CHKERRQ(ierr); 00035 00036 ierr = PetscStrcmp(type,"identity",&flg); CHKERRQ(ierr); 00037 if (flg) { 00038 goto done; 00039 } 00040 00041 00042 ierr = PetscStrcmp(type,"eliminate",&flg); CHKERRQ(ierr); 00043 if (flg) { 00044 Mat A,Ause, B,Buse; Vec X,Xuse, V,Vuse, I,Iuse; 00045 MPI_Comm comm; 00046 IS dummy_rows,all_rows,true_rows,gtrue_rows,compact_rows; 00047 VecScatter compact_extractor; 00048 int ndummy,dummy_kind,first,last,local_size,compact_size,nd,*dd; 00049 singleton_struct *savestruct; 00050 00051 ierr = LinearSystemGetParts(problem,&A,&B,&V,&X,&I); CHKERRQ(ierr); 00052 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00053 ierr = PetscMalloc(sizeof(singleton_struct),&savestruct); CHKERRQ(ierr); 00054 savestruct->n = 1; 00055 00056 /* 00057 * Get dummy row information 00058 */ 00059 ierr = SysProRetrieveQuantity 00060 (inproblem,"structure","dummy-rows",(void*)&dd,PETSC_NULL, 00061 &flg); CHKERRQ(ierr); 00062 if (!flg) goto no_compute; 00063 ierr = SysProRetrieveQuantity 00064 (inproblem,"structure","n-dummy-rows",(void*)&ndummy,PETSC_NULL, 00065 &flg); CHKERRQ(ierr); 00066 if (!flg) goto no_compute; 00067 ierr = SysProRetrieveQuantity 00068 (inproblem,"structure","dummy-rows-kind",(void*)&dummy_kind,PETSC_NULL, 00069 &flg); CHKERRQ(ierr); 00070 if (!flg) goto no_compute; 00071 savestruct->t = dummy_kind; 00072 goto eliminate; 00073 00074 no_compute: 00075 /* 00076 * No dummy row information; take default action, then exit 00077 * WRONG this will crash in the back substitution 00078 */ 00079 goto done; 00080 eliminate: 00081 /* 00082 * Eliminate dummy rows 00083 */ 00084 #if 0 00085 if (ndummy!=dd[0]) 00086 SETERRQ2(1,"Dummy rows size mismatch",ndummy,dd[0]); 00087 dd++; /* skip past zero element which is array size */ 00088 #endif 00089 ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr); 00090 local_size = last-first; 00091 ierr = ISCreateStride 00092 (comm,local_size,first,1,&all_rows); CHKERRQ(ierr); 00093 for (nd=0; nd<ndummy; nd++) if (dd[nd]>=first) break; 00094 dd = dd+nd; ndummy -= nd; 00095 for (nd=ndummy-1; nd>=0; nd--) if (dd[nd]>=last) ndummy--; else break; 00096 compact_size = local_size-ndummy; 00097 ierr = ISCreateGeneral 00098 (comm,ndummy,dd,&dummy_rows); CHKERRQ(ierr); 00099 printf("dummy rows\n"); ISView(dummy_rows,0); 00100 ierr = ISDifference(all_rows,dummy_rows,&true_rows); CHKERRQ(ierr); 00101 ierr = ISAllGather(true_rows,>rue_rows); CHKERRQ(ierr); 00102 00103 /* 00104 * Extract the true submatrix 00105 */ 00106 ierr = MatGetSubMatrix 00107 (A,true_rows,gtrue_rows,MAT_INITIAL_MATRIX, 00108 &Ause); CHKERRQ(ierr); 00109 if (B && B!=A) { 00110 ierr = MatGetSubMatrix 00111 (B,true_rows,gtrue_rows,MAT_INITIAL_MATRIX, 00112 &Buse); CHKERRQ(ierr); 00113 } else Buse = Ause; 00114 ierr = ISDestroy(gtrue_rows); CHKERRQ(ierr); 00115 /* and find out its distribution */ 00116 ierr = MatGetOwnershipRange(Ause,&first,&last); CHKERRQ(ierr); 00117 00118 /* 00119 * Extract the true vector 00120 */ 00121 /* create compact vector */ 00122 ierr = VecCreateMPI 00123 (comm,compact_size,PETSC_DECIDE,&Xuse); CHKERRQ(ierr); 00124 ierr = VecDuplicate(Xuse,&Vuse); CHKERRQ(ierr); 00125 ierr = VecDuplicate(Xuse,&Iuse); CHKERRQ(ierr); 00126 00127 /* create scatter */ 00128 ierr = ISCreateStride 00129 (MPI_COMM_SELF,compact_size,first,1,&compact_rows); CHKERRQ(ierr); 00130 ierr = VecScatterCreate 00131 (X,true_rows,Xuse,compact_rows,&compact_extractor); CHKERRQ(ierr); 00132 savestruct->extractor = compact_extractor; 00133 ierr = ISDestroy(compact_rows); CHKERRQ(ierr); 00134 ierr = ISDestroy(dummy_rows); CHKERRQ(ierr); 00135 ierr = ISDestroy(all_rows); CHKERRQ(ierr); 00136 ierr = ISDestroy(true_rows); CHKERRQ(ierr); 00137 00138 /* compact rhs */ 00139 ierr = VecScatterBegin 00140 (compact_extractor,V,Vuse,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr); 00141 ierr = VecScatterEnd 00142 (compact_extractor,V,Vuse,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr); 00143 /* compact true solution */ 00144 ierr = VecScatterBegin 00145 (compact_extractor,X,Xuse,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr); 00146 ierr = VecScatterEnd 00147 (compact_extractor,X,Xuse,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr); 00148 /* compact starting guess */ 00149 if (I) { 00150 ierr = VecScatterBegin 00151 (compact_extractor,I,Iuse,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 00152 ierr = VecScatterEnd 00153 (compact_extractor,I,Iuse,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 00154 } 00155 ierr = LinearSystemSetParts 00156 (newproblem,Ause,Buse,Vuse,Xuse,Iuse); CHKERRQ(ierr); 00157 { 00158 void *info; 00159 ierr = LinearSystemGetContext(problem,&info); CHKERRQ(ierr); 00160 ierr = LinearSystemSetContext(newproblem,info); CHKERRQ(ierr); 00161 } 00162 *(singleton_struct**)ctx = savestruct; 00163 goto done; 00164 } 00165 00166 SETERRQ1(1,"Unknown singleton preprocessor <%s>",type); 00167 done: 00168 *outproblem = (NumericalProblem)newproblem; 00169 PetscFunctionReturn(0); 00170 } 00171 00172 #undef __FUNCT__ 00173 #define __FUNCT__ "back_singleton" 00174 static PetscErrorCode back_singleton 00175 (const char *singleton_type,PetscTruth overwrite,void *gctx,void *ctx, 00176 NumericalProblem compactproblem,NumericalProblem fullproblem, 00177 NumericalSolution compactvector,NumericalSolution fullvector) 00178 { 00179 LinearSolution compactsolution = (LinearSolution)compactvector, 00180 fullsolution = (LinearSolution)fullvector; 00181 Vec Xcompact,Xfull,Rhs; 00182 singleton_struct* savestruct; int dummy_kind; 00183 VecScatter compact_extractor; 00184 PetscTruth flg; PetscErrorCode ierr; 00185 00186 PetscFunctionBegin; 00187 00188 ierr = LinearSolutionCopyStats(compactsolution,fullsolution); CHKERRQ(ierr); 00189 ierr = LinearSolutionGetVector(compactsolution,&Xcompact); CHKERRQ(ierr); 00190 ierr = LinearSolutionGetVector(fullsolution,&Xfull); CHKERRQ(ierr); 00191 00192 ierr = PetscStrcmp(singleton_type,"identity",&flg); CHKERRQ(ierr); 00193 if (flg) { 00194 ierr = VecCopy(Xcompact,Xfull); CHKERRQ(ierr); 00195 } else { 00196 ierr = LinearSystemGetParts 00197 ((LinearSystem)fullproblem,PETSC_NULL,PETSC_NULL,&Rhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 00198 00199 savestruct = (singleton_struct*) ctx; 00200 dummy_kind = savestruct->t; 00201 if (dummy_kind>0) 00202 printf("\n\n>>>> WARNING Cannot properly handle dummy rows of kind %d <<<<\n\n",dummy_kind); 00203 compact_extractor = savestruct->extractor; 00204 00205 /* 00206 * first copy dummy rows from rhs to solution 00207 */ 00208 ierr = VecDuplicate(Rhs,&Xfull); CHKERRQ(ierr); 00209 ierr = VecCopy(Rhs,Xfull); CHKERRQ(ierr); 00210 /* then insert computed solution on non-dummy variables */ 00211 ierr = VecScatterBegin(compact_extractor,Xcompact,Xfull,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr); 00212 ierr = VecScatterEnd (compact_extractor,Xcompact,Xfull,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr); 00213 ierr = LinearSolutionSetVector(fullsolution,Xfull); CHKERRQ(ierr); 00214 00215 ierr = VecScatterDestroy(compact_extractor); CHKERRQ(ierr); 00216 ierr = PetscFree(savestruct); CHKERRQ(ierr); 00217 } 00218 00219 PetscFunctionReturn(0); 00220 } 00221 00222 #undef __FUNCT__ 00223 #define __FUNCT__ "setup_singleton_choices" 00224 /*! This routine is only called when the singleton preprocessor 00225 is created by DeclarePreprocessor() inside DeclareSingletonPreprocessor() 00226 */ 00227 static PetscErrorCode setup_singleton_choices() 00228 { 00229 SalsaTransform tf; SalsaTransformObject cur; PetscErrorCode ierr; 00230 00231 PetscFunctionBegin; 00232 ierr = TransformGetByName(PREPROCESSOR,&tf); CHKERRQ(ierr); 00233 00234 ierr = NewTransformObject(PREPROCESSOR,"identity",&cur); CHKERRQ(ierr); 00235 ierr = TransformObjectSetExplanation 00236 (cur,"no singletons"); CHKERRQ(ierr); 00237 00238 ierr = NewTransformObject(PREPROCESSOR,"eliminate",&cur); CHKERRQ(ierr); 00239 ierr = TransformObjectSetExplanation 00240 (cur,"eliminate singletons"); CHKERRQ(ierr); 00241 00242 PetscFunctionReturn(0); 00243 } 00244 00245 #undef __FUNCT__ 00246 #define __FUNCT__ "specific_singleton_choices" 00247 /*! 00248 This is the 'specific setup' phase of the singleton preprocessor. 00249 See \ref modes for details. 00250 00251 It disables either the identity or the elimination routine, to 00252 leave only the one applicable to this particular system. 00253 00254 Maybe if we'd ever want to prove how effective singleton elimination 00255 is, we could leave identity in place for systems with singletons. 00256 */ 00257 static PetscErrorCode specific_singleton_choices 00258 (NumericalProblem theproblem,SalsaTransform singleton) 00259 { 00260 SalsaTransformObject tf; PetscTruth flg; int ndummy,nnz; 00261 PetscErrorCode ierr; 00262 PetscFunctionBegin; 00263 00264 /* 00265 * Mark (as unsuitable) the identity if there are dummy rows, 00266 * and the matrix is not diagonal. 00267 */ 00268 ierr = SysProComputeQuantity 00269 (theproblem,"structure","n-dummy-rows",(void*)&ndummy,PETSC_NULL, 00270 &flg); CHKERRQ(ierr); 00271 if (flg && ndummy>0) { 00272 ierr = SysProComputeQuantity 00273 (theproblem,"structure","nnzeros",(void*)&nnz,PETSC_NULL, 00274 &flg); CHKERRQ(ierr); 00275 if (flg && nnz>ndummy) { 00276 ierr = TransformObjectGetByName 00277 (PREPROCESSOR,"identity",&tf); CHKERRQ(ierr); 00278 ierr = TransformObjectMark(tf); CHKERRQ(ierr); 00279 goto done; 00280 } 00281 } 00282 00283 /* Otherwise, mark elimination for exclusion */ 00284 ierr = TransformObjectGetByName(PREPROCESSOR,"eliminate",&tf); CHKERRQ(ierr); 00285 ierr = TransformObjectMark(tf); CHKERRQ(ierr); 00286 done: 00287 00288 PetscFunctionReturn(0); 00289 } 00290 00291 static PetscErrorCode singleton_specific_unset(NumericalProblem theproblem) 00292 { 00293 PetscErrorCode ierr; PetscTruth flg; 00294 PetscFunctionBegin; 00295 ierr = SysProRemoveQuantity 00296 (theproblem,"structure","n-dummy-rows",&flg); CHKERRQ(ierr); 00297 ierr = SysProRemoveQuantity 00298 (theproblem,"structure","nnzeros",&flg); CHKERRQ(ierr); 00299 PetscFunctionReturn(0); 00300 } 00301 00302 #undef __FUNCT__ 00303 #define __FUNCT__ "DeclareSingletonPreprocessor" 00304 PetscErrorCode DeclareSingletonPreprocessor(void) 00305 { 00306 PetscErrorCode ierr; 00307 PetscFunctionBegin; 00308 ierr = DeclarePreprocessor 00309 (PREPROCESSOR, 00310 &setup_singleton_choices,&specific_singleton_choices, 00311 &singleton_specific_unset,0, 00312 0,0, 00313 &eliminate_singletons,&back_singleton); CHKERRQ(ierr); 00314 ierr = PreprocessorSetPreservedCategories 00315 (PREPROCESSOR,"laplace"); CHKERRQ(ierr); 00316 PetscFunctionReturn(0); 00317 }