System Preprocessors
distribution.c
Go to the documentation of this file.
00001 /*! \file distribution.c \ingroup linear */
00002 
00003 /*! \page distribution Permutation / load balancing
00004 
00005   Linear system solving is sensitive in several ways to permutations
00006   and load balancing applied to the system. This dependency typically comes
00007   through the preconditioner: incomplete factorizations are sensitive
00008   to permutations, and block Jacobi and Schwarz preconditioners are sensitive 
00009   to load distributions, even without any permutation applied.
00010  */
00011 #include <stdlib.h>
00012 #include "petscmat.h"
00013 #include "petscconf.h"
00014 #include "syspro.h"
00015 #include "sysprotransform.h"
00016 #include "sysprolinear.h"
00017 #include "linear_impl.h"
00018 #include "anamod.h"
00019 
00020 #define PREPROCESSOR "distribution"
00021 extern int SpectrumComputeUnpreconditionedSpectrum();
00022 
00023 #undef __FUNCT__
00024 #define __FUNCT__ "setup_distribution_choices"
00025 static PetscErrorCode setup_distribution_choices()
00026 {
00027   SalsaTransform tf; SalsaTransformObject cur;
00028   PetscTruth flg; PetscErrorCode ierr;
00029 
00030   PetscFunctionBegin;
00031   ierr = TransformGetByName(PREPROCESSOR,&tf); CHKERRQ(ierr);
00032 
00033   ierr = SysProDefineIntAnnotation
00034     (PREPROCESSOR,"preserve-simple"); CHKERRQ(ierr);
00035   ierr = SysProDefineIntAnnotation
00036     (PREPROCESSOR,"strictly-parallel"); CHKERRQ(ierr);
00037 
00038   ierr = NewTransformObject(PREPROCESSOR,"induced",&cur); CHKERRQ(ierr);
00039   ierr = TransformObjectSetExplanation
00040     (cur,"default distribution"); CHKERRQ(ierr);
00041   ierr = TransformObjectIntAnnotate(cur,"preserve-simple",1); CHKERRQ(ierr);
00042 
00043   ierr = HasComputeModule("icmk","splits",&flg); CHKERRQ(ierr);
00044   if (flg) {
00045     ierr = NewTransformObject(PREPROCESSOR,"icmk",&cur); CHKERRQ(ierr);
00046     ierr = TransformObjectSetExplanation
00047       (cur,"inverse cuthill-mckee"); CHKERRQ(ierr);
00048     ierr = TransformObjectIntAnnotate(cur,"strictly-parallel",1); CHKERRQ(ierr);
00049   } else printf("no icmk available\n");
00050 
00051 #if defined(PETSC_HAVE_PARMETIS)
00052   ierr = NewTransformObject
00053     (PREPROCESSOR,MAT_PARTITIONING_PARMETIS,&cur); CHKERRQ(ierr);
00054   ierr = TransformObjectSetExplanation
00055     (cur,"parmetis package"); CHKERRQ(ierr);
00056 #endif
00057 
00058 #if defined(PETSC_HAVE_CHACO)
00059   ierr = NewTransformObject(PREPROCESSOR,"chaco",&cur); CHKERRQ(ierr);
00060   ierr = TransformObjectSetExplanation
00061     (cur,"chaco package"); CHKERRQ(ierr);
00062   ierr = TransformObjectIntAnnotate(cur,"strictly-parallel",1); CHKERRQ(ierr);
00063 #endif
00064 
00065   PetscFunctionReturn(0);
00066 }
00067 
00068 static PetscErrorCode specific_distribution_choices
00069 (NumericalProblem problem,SalsaTransform tf)
00070 {
00071   MPI_Comm comm = problem->comm; int ntids;
00072   int i,n; SalsaTransformObject *objs; PetscErrorCode ierr;
00073   PetscFunctionBegin;
00074   MPI_Comm_size(comm,&ntids);
00075   ierr = TransformGetObjects(tf,&n,&objs); CHKERRQ(ierr);
00076   for (i=0; i<n; i++) {
00077     SalsaTransformObject cur = objs[i]; 
00078     int v; PetscTruth f;
00079     ierr = TransformObjectGetIntAnnotation
00080       (cur,"strictly-parallel",&v,&f); CHKERRQ(ierr);
00081     if (f && v) {
00082       ierr = TransformObjectMark(cur); CHKERRQ(ierr);
00083     }
00084   }
00085   PetscFunctionReturn(0);
00086 }
00087 
00088 #undef __FUNCT__
00089 #define __FUNCT__ "sans_partition"
00090 static PetscErrorCode sans_partition
00091 (const char *type,NumericalProblem inproblem,int nparts,
00092  IS *local_to_global,VecScatter *perm)
00093 {
00094   LinearSystem problem = (LinearSystem)inproblem;
00095   Mat A; Vec X; MPI_Comm comm; PetscTruth flg;
00096   MatPartitioning Part;
00097   IS local,local_to_proc,new_part;
00098   int *part_sizes,local_size,first,last,ntids,mytid, ierr;
00099   
00100   PetscFunctionBegin;
00101   ierr = LinearSystemGetParts
00102     (problem,&A,PETSC_NULL,PETSC_NULL,&X,PETSC_NULL); CHKERRQ(ierr);
00103   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00104   ierr = MPI_Comm_size(comm,&ntids); CHKERRQ(ierr);
00105   ierr = MPI_Comm_rank(comm,&mytid); CHKERRQ(ierr);
00106   ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00107   local_size = last-first;
00108   ierr = VecCreateMPI(comm,local_size,PETSC_DECIDE,&X); CHKERRQ(ierr);
00109   ierr = PetscMalloc((nparts+1)*sizeof(int),&part_sizes); CHKERRQ(ierr);
00110   ierr = ISCreateStride(comm,local_size,first,1,&local); CHKERRQ(ierr);
00111 
00112   flg = PETSC_FALSE;
00113 #if defined(PETSC_HAVE_PARMETIS)
00114   if (!flg) {
00115     ierr = PetscStrcmp(type,"parmetis",&flg); CHKERRQ(ierr);}
00116 #endif
00117 #if defined(PETSC_HAVE_CHACO)
00118   if (!flg) {
00119     ierr = PetscStrcmp(type,"chaco",&flg); CHKERRQ(ierr);}
00120 #endif
00121   if (flg) {
00122     IS tis; PetscInt local_size;
00123 
00124     /*ierr = MatConvert(A,MATMPIADJ,MAT_INITIAL_MATRIX,&Adj); CHKERRQ(ierr);*/
00125     ierr = MatPartitioningCreate(comm,&Part); CHKERRQ(ierr);
00126     ierr = MatPartitioningRegisterAll(PETSC_NULL); CHKERRQ(ierr);
00127     ierr = MatPartitioningSetType(Part,type); CHKERRQ(ierr);
00128     ierr = MatPartitioningSetAdjacency(Part,A/* Adj? */); CHKERRQ(ierr);
00129 
00130     ierr = MatPartitioningSetNParts(Part,nparts); CHKERRQ(ierr);
00131     ierr = MatPartitioningApply(Part,&local_to_proc); CHKERRQ(ierr);
00132     ierr = MatPartitioningDestroy(Part); CHKERRQ(ierr);
00133     /*ierr = MatDestroy(Adj); CHKERRQ(ierr);*/
00134 
00135     ierr = ISPartitioningToNumbering(local_to_proc,&tis); CHKERRQ(ierr);
00136     ierr = ISPartitioningCount(local_to_proc,nparts,part_sizes); CHKERRQ(ierr);
00137     ierr = ISDestroy(local_to_proc); CHKERRQ(ierr);
00138     ierr = MatGetLocalSize(A,&local_size,PETSC_NULL); CHKERRQ(ierr);
00139     ierr = ISInvertPermutation(tis,local_size,local_to_global); CHKERRQ(ierr);
00140     /*printf("mapping\n"); ISView(local_to_global,0);*/
00141     
00142     goto size_to_part;
00143   }
00144 
00145   ierr = PetscStrcmp(type,"icmk",&flg); CHKERRQ(ierr);
00146   if (flg) {
00147     int *ss;
00148     *local_to_global = local;
00149     ierr = SysProComputeQuantity
00150       (inproblem,"icmk","splits",(void*)&ss,PETSC_NULL,
00151        &flg); CHKERRQ(ierr);
00152     if (flg) {
00153       ierr = ISCreateGeneral
00154         (MPI_COMM_SELF,ss[0],ss+1,&new_part); CHKERRQ(ierr);
00155     } else {
00156       SETERRQ(1,"Hm.....");
00157     }
00158     goto scatter;
00159   }    
00160 
00161   SETERRQ1(1,"Unknown distribution <%s>",type);
00162 
00163  size_to_part:
00164   {
00165     int i,t,s=0;
00166     for (i=0; i<nparts; i++) {
00167       t = part_sizes[i]; part_sizes[i] = s; s+=t;}
00168     part_sizes[nparts] = s;
00169   }
00170 
00171   ierr = ISCreateGeneral
00172     (MPI_COMM_SELF,nparts+1,part_sizes,&new_part); CHKERRQ(ierr);
00173   ierr = PetscFree(part_sizes); CHKERRQ(ierr);
00174  scatter:
00175   ierr = VecScatterCreate(X,local,X,*local_to_global,perm); CHKERRQ(ierr);
00176   if (local!=*local_to_global) {
00177     ierr = ISDestroy(local); CHKERRQ(ierr);
00178   }
00179 
00180   PetscFunctionReturn(0);
00181 }
00182 
00183 #undef __FUNCT__
00184 #define __FUNCT__ "distribute_system"
00185 static PetscErrorCode distribute_system
00186 (const char *type,int nopt,PetscTruth overwrite,
00187  NumericalProblem inproblem,NumericalProblem *outproblem,
00188  void *gctx,void **ctx, PetscTruth *success)
00189 {
00190   LinearSystem problem = (LinearSystem)inproblem,newproblem;
00191   PetscTruth flg; PetscErrorCode ierr; 
00192 
00193   PetscFunctionBegin;
00194   *success = PETSC_TRUE;
00195 
00196   ierr = PetscStrcmp(type,"induced",&flg); CHKERRQ(ierr);
00197   if (flg) {
00198     ierr = LinearSystemDuplicatePointers(problem,&newproblem); CHKERRQ(ierr);
00199     *(VecScatter*)ctx = NULL;
00200     goto done;
00201   }
00202 
00203   {
00204     Mat A,Ause,B,Buse; Vec X,Xuse,I,Iuse,V,Vuse;
00205     VecScatter perm;
00206 
00207     ierr = LinearSystemDuplicatePointers(problem,&newproblem); CHKERRQ(ierr);
00208     ierr = LinearSystemGetParts(problem,&A,&B,&V,&X,&I); CHKERRQ(ierr);
00209 
00210     {
00211       IS local_to_global=0,global_to_global=0;
00212       ierr = sans_partition(type,inproblem,6,&local_to_global,&perm); CHKERRQ(ierr);
00213       ierr = ISAllGather(local_to_global,&global_to_global); CHKERRQ(ierr);
00214       ierr = ISSetPermutation(local_to_global); CHKERRQ(ierr);
00215       ierr = ISSetPermutation(global_to_global); CHKERRQ(ierr);
00216       ierr = MatPermute
00217         (A,local_to_global,global_to_global,&Ause); CHKERRQ(ierr);
00218       if (B && B!=A) {
00219         ierr = MatPermute
00220           (B,local_to_global,local_to_global,&Buse); CHKERRQ(ierr);
00221       } else Buse = Ause;
00222       ierr = ISDestroy(local_to_global); CHKERRQ(ierr);
00223       ierr = ISDestroy(global_to_global); CHKERRQ(ierr);
00224     }
00225     *(VecScatter*)ctx = perm;
00226 
00227     {
00228       int local_size; MPI_Comm comm;
00229       ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00230       ierr = MatGetLocalSize(Ause,&local_size,PETSC_NULL); CHKERRQ(ierr);
00231       ierr = VecCreateMPI(comm,local_size,PETSC_DECIDE,&Xuse); CHKERRQ(ierr);
00232       ierr = VecDuplicate(Xuse,&Iuse); CHKERRQ(ierr);
00233       ierr = VecDuplicate(Xuse,&Vuse); CHKERRQ(ierr);
00234     }
00235 
00236     /* permute the rhs */
00237     ierr = VecScatterBegin
00238       (perm,V,Vuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00239     ierr = VecScatterEnd
00240       (perm,V,Vuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00241     /* permute known solution */
00242       ierr = VecScatterBegin
00243         (perm,X,Xuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00244     ierr = VecScatterEnd
00245       (perm,X,Xuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00246     /* permute initial guess */
00247     ierr = VecScatterBegin
00248       (perm,I,Iuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00249     ierr = VecScatterEnd
00250       (perm,I,Iuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00251   
00252     ierr = LinearSystemSetParts
00253       (newproblem,Ause,Buse,Vuse,Xuse,Iuse); CHKERRQ(ierr);
00254 
00255     goto done;
00256   }
00257 
00258   SETERRQ1(1,"Undefined distribution <%s>",type);
00259 
00260  done:
00261   *outproblem = (NumericalProblem)newproblem;
00262   PetscFunctionReturn(0);
00263 }
00264 
00265 #undef __FUNCT__
00266 #define __FUNCT__ "undistribute_system"
00267 static PetscErrorCode undistribute_system
00268 (const char *scaling_type,PetscTruth overwrite,void *gctx,void *ctx,
00269  NumericalProblem problem,NumericalProblem nextproblem,
00270  NumericalSolution before,NumericalSolution after)
00271 {
00272   LinearSolution shuffled = (LinearSolution)before, 
00273     rectified = (LinearSolution)after;
00274   VecScatter perm = (VecScatter)ctx; Vec Xfore,Xback; 
00275   PetscErrorCode ierr;
00276   PetscFunctionBegin;
00277 
00278   ierr = LinearSolutionCopyStats(shuffled,rectified); CHKERRQ(ierr);
00279   ierr = LinearSolutionGetVector(shuffled,&Xfore); CHKERRQ(ierr);
00280   ierr = LinearSolutionGetVector(rectified,&Xback); CHKERRQ(ierr);
00281   if (perm) {
00282     ierr = VecScatterBegin
00283       (perm,Xfore,Xback,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00284     ierr = VecScatterEnd
00285       (perm,Xfore,Xback,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00286     ierr = VecScatterDestroy(perm); CHKERRQ(ierr);
00287   } else {
00288     ierr = VecCopy(Xfore,Xback); CHKERRQ(ierr);
00289   }
00290 
00291   /*ierr = LinearSolutionDelete((LinearSolution)before); CHKERRQ(ierr);*/
00292   //ierr = DeleteLinearSystem((LinearSystem)problem); CHKERRQ(ierr);
00293   PetscFunctionReturn(0);
00294 }
00295 
00296 #undef __FUNCT__
00297 #define __FUNCT__ "DeclareDistributionPreprocessor"
00298 PetscErrorCode DeclareDistributionPreprocessor(void)
00299 {
00300   PetscErrorCode ierr;
00301   PetscFunctionBegin;
00302   ierr = DeclarePreprocessor
00303     ("distribution",
00304      &setup_distribution_choices,&specific_distribution_choices,0,0,
00305      0,0,
00306      &distribute_system,&undistribute_system); CHKERRQ(ierr);
00307   ierr = PreprocessorSetPreservedCategories
00308     (PREPROCESSOR,"laplace"); CHKERRQ(ierr);
00309   PetscFunctionReturn(0);
00310 }