00001
00002
00003
00004
00005
00006
00007
00008
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 (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
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
00077
00078
00079 goto done;
00080 eliminate:
00081
00082
00083
00084 #if 0
00085 if (ndummy!=dd[0])
00086 SETERRQ2(1,"Dummy rows size mismatch",ndummy,dd[0]);
00087 dd++;
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
00105
00106 ierr = MatGetSubMatrix
00107 (A,true_rows,gtrue_rows,compact_size,MAT_INITIAL_MATRIX,
00108 &Ause); CHKERRQ(ierr);
00109 if (B && B!=A) {
00110 ierr = MatGetSubMatrix
00111 (B,true_rows,gtrue_rows,compact_size,MAT_INITIAL_MATRIX,
00112 &Buse); CHKERRQ(ierr);
00113 } else Buse = Ause;
00114
00115 ierr = MatGetOwnershipRange(Ause,&first,&last); CHKERRQ(ierr);
00116
00117
00118
00119
00120
00121 ierr = VecCreateMPI
00122 (comm,compact_size,PETSC_DECIDE,&Xuse); CHKERRQ(ierr);
00123 ierr = VecDuplicate(Xuse,&Vuse); CHKERRQ(ierr);
00124 ierr = VecDuplicate(Xuse,&Iuse); CHKERRQ(ierr);
00125
00126
00127 ierr = ISCreateStride
00128 (MPI_COMM_SELF,compact_size,first,1,&compact_rows); CHKERRQ(ierr);
00129 ierr = VecScatterCreate
00130 (X,true_rows,Xuse,compact_rows,&compact_extractor); CHKERRQ(ierr);
00131 savestruct->extractor = compact_extractor;
00132 ierr = ISDestroy(compact_rows); CHKERRQ(ierr);
00133 ierr = ISDestroy(dummy_rows); CHKERRQ(ierr);
00134 ierr = ISDestroy(all_rows); CHKERRQ(ierr);
00135 ierr = ISDestroy(true_rows); CHKERRQ(ierr);
00136
00137
00138 ierr = VecScatterBegin
00139 (compact_extractor,V,Vuse,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00140 ierr = VecScatterEnd
00141 (compact_extractor,V,Vuse,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00142
00143 ierr = VecScatterBegin
00144 (compact_extractor,X,Xuse,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00145 ierr = VecScatterEnd
00146 (compact_extractor,X,Xuse,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00147
00148 if (I) {
00149 ierr = VecScatterBegin
00150 (compact_extractor,I,Iuse,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
00151 ierr = VecScatterEnd
00152 (compact_extractor,I,Iuse,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
00153 }
00154 ierr = LinearSystemSetParts
00155 (newproblem,Ause,Buse,Vuse,Xuse,Iuse); CHKERRQ(ierr);
00156 {
00157 void *info;
00158 ierr = LinearSystemGetContext(problem,&info); CHKERRQ(ierr);
00159 ierr = LinearSystemSetContext(newproblem,info); CHKERRQ(ierr);
00160 }
00161 *(singleton_struct**)ctx = savestruct;
00162 goto done;
00163 }
00164
00165 SETERRQ1(1,"Unknown singleton preprocessor <%s>",type);
00166 done:
00167 *outproblem = (NumericalProblem)newproblem;
00168 PetscFunctionReturn(0);
00169 }
00170
00171 #undef __FUNCT__
00172 #define __FUNCT__ "back_singleton"
00173 static PetscErrorCode back_singleton
00174 (char *singleton_type,PetscTruth overwrite,void *gctx,void *ctx,
00175 NumericalProblem compactproblem,NumericalProblem fullproblem,
00176 NumericalSolution compactvector,NumericalSolution fullvector)
00177 {
00178 LinearSolution compactsolution = (LinearSolution)compactvector,
00179 fullsolution = (LinearSolution)fullvector;
00180 Vec Xcompact,Xfull,Rhs;
00181 singleton_struct* savestruct; int dummy_kind;
00182 VecScatter compact_extractor;
00183 PetscTruth flg; PetscErrorCode ierr;
00184
00185 PetscFunctionBegin;
00186
00187 ierr = LinearSolutionCopyStats(compactsolution,fullsolution); CHKERRQ(ierr);
00188 ierr = LinearSolutionGetVector(compactsolution,&Xcompact); CHKERRQ(ierr);
00189 ierr = LinearSolutionGetVector(fullsolution,&Xfull); CHKERRQ(ierr);
00190
00191 ierr = PetscStrcmp(singleton_type,"identity",&flg); CHKERRQ(ierr);
00192 if (flg) {
00193 ierr = VecCopy(Xcompact,Xfull); CHKERRQ(ierr);
00194 } else {
00195
00196 ierr = LinearSystemGetParts
00197 ((LinearSystem)fullproblem,
00198 PETSC_NULL,PETSC_NULL,&Rhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00199
00200 savestruct = (singleton_struct*) ctx;
00201 dummy_kind = savestruct->t;
00202 if (dummy_kind>0)
00203 printf("\n\n>>>> WARNING Cannot properly handle dummy rows of kind %d <<<<\n\n",
00204 dummy_kind);
00205 compact_extractor = savestruct->extractor;
00206
00207
00208
00209
00210 ierr = VecDuplicate(Rhs,&Xfull); CHKERRQ(ierr);
00211 ierr = VecCopy(Rhs,Xfull); CHKERRQ(ierr);
00212
00213 ierr = VecScatterBegin
00214 (compact_extractor,Xcompact,Xfull,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00215 ierr = VecScatterEnd
00216 (compact_extractor,Xcompact,Xfull,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00217 ierr = LinearSolutionSetVector(fullsolution,Xfull); CHKERRQ(ierr);
00218
00219 ierr = VecScatterDestroy(compact_extractor); CHKERRQ(ierr);
00220 ierr = PetscFree(savestruct); CHKERRQ(ierr);
00221 }
00222
00223
00224 ierr = DeleteLinearSystem((LinearSystem)compactproblem); CHKERRQ(ierr);
00225
00226 PetscFunctionReturn(0);
00227 }
00228
00229 #undef __FUNCT__
00230 #define __FUNCT__ "setup_singleton_choices"
00231
00232
00233
00234 static PetscErrorCode setup_singleton_choices()
00235 {
00236 SalsaTransform tf; SalsaTransformObject cur; PetscErrorCode ierr;
00237
00238 PetscFunctionBegin;
00239 ierr = TransformGetByName(PREPROCESSOR,&tf); CHKERRQ(ierr);
00240
00241 ierr = NewTransformObject(PREPROCESSOR,"identity",&cur); CHKERRQ(ierr);
00242 ierr = TransformObjectSetExplanation
00243 (cur,"no singletons"); CHKERRQ(ierr);
00244
00245 ierr = NewTransformObject(PREPROCESSOR,"eliminate",&cur); CHKERRQ(ierr);
00246 ierr = TransformObjectSetExplanation
00247 (cur,"eliminate singletons"); CHKERRQ(ierr);
00248
00249 PetscFunctionReturn(0);
00250 }
00251
00252 #undef __FUNCT__
00253 #define __FUNCT__ "specific_singleton_choices"
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264 static PetscErrorCode specific_singleton_choices
00265 (NumericalProblem theproblem,SalsaTransform singleton)
00266 {
00267 SalsaTransformObject tf; PetscTruth flg; int ndummy,nnz;
00268 PetscErrorCode ierr;
00269 PetscFunctionBegin;
00270
00271
00272
00273
00274
00275 ierr = SysProComputeQuantity
00276 (theproblem,"structure","n-dummy-rows",(void*)&ndummy,PETSC_NULL,
00277 &flg); CHKERRQ(ierr);
00278 if (flg && ndummy>0) {
00279 ierr = SysProComputeQuantity
00280 (theproblem,"structure","nnzeros",(void*)&nnz,PETSC_NULL,
00281 &flg); CHKERRQ(ierr);
00282 if (flg && nnz>ndummy) {
00283 ierr = TransformObjectGetByName
00284 (PREPROCESSOR,"identity",&tf); CHKERRQ(ierr);
00285 ierr = TransformObjectMark(tf); CHKERRQ(ierr);
00286 goto done;
00287 }
00288 }
00289
00290
00291 ierr = TransformObjectGetByName(PREPROCESSOR,"eliminate",&tf); CHKERRQ(ierr);
00292 ierr = TransformObjectMark(tf); CHKERRQ(ierr);
00293 done:
00294
00295 PetscFunctionReturn(0);
00296 }
00297
00298 static PetscErrorCode singleton_specific_unset(NumericalProblem theproblem)
00299 {
00300 PetscErrorCode ierr; PetscTruth flg;
00301 PetscFunctionBegin;
00302 ierr = SysProRemoveQuantity
00303 (theproblem,"structure","n-dummy-rows",&flg); CHKERRQ(ierr);
00304 ierr = SysProRemoveQuantity
00305 (theproblem,"structure","nnzeros",&flg); CHKERRQ(ierr);
00306 PetscFunctionReturn(0);
00307 }
00308
00309 #undef __FUNCT__
00310 #define __FUNCT__ "DeclareSingletonPreprocessor"
00311 PetscErrorCode DeclareSingletonPreprocessor(void)
00312 {
00313 PetscErrorCode ierr;
00314 PetscFunctionBegin;
00315 ierr = DeclarePreprocessor
00316 (PREPROCESSOR,
00317 &setup_singleton_choices,&specific_singleton_choices,
00318 &singleton_specific_unset,0,
00319 0,0,
00320 &eliminate_singletons,&back_singleton); CHKERRQ(ierr);
00321 ierr = PreprocessorSetPreservedCategories
00322 (PREPROCESSOR,"laplace"); CHKERRQ(ierr);
00323 PetscFunctionReturn(0);
00324 }