System Preprocessors
singleton.c
Go to the documentation of this file.
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,&gtrue_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 }