System Preprocessors
|
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 PetscBool 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,"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; PetscBool 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; PetscBool 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,PETSC_OWN_POINTER,&new_part); CHKERRQ(ierr); 00155 } else { 00156 SETERRQ(PETSC_COMM_SELF,1,"Hm....."); 00157 } 00158 goto scatter; 00159 } 00160 00161 SETERRQ1(PETSC_COMM_SELF,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,PETSC_OWN_POINTER,&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,PetscBool overwrite, 00187 NumericalProblem inproblem,NumericalProblem *outproblem, 00188 void *gctx,void **ctx, PetscBool *success) 00189 { 00190 LinearSystem problem = (LinearSystem)inproblem,newproblem; 00191 PetscBool 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(PETSC_COMM_SELF,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,PetscBool 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 }