00001
00002
00003
00004
00005
00006
00007
00008
00009
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 (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
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); 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
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
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 (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,global_to_global;
00212 ierr = sans_partition
00213 (type,inproblem,6,&local_to_global,&perm); CHKERRQ(ierr);
00214 ierr = ISAllGather(local_to_global,&global_to_global); CHKERRQ(ierr);
00215 ierr = ISSetPermutation(local_to_global); CHKERRQ(ierr);
00216 ierr = ISSetPermutation(global_to_global); CHKERRQ(ierr);
00217 ierr = MatPermute
00218 (A,local_to_global,global_to_global,&Ause); CHKERRQ(ierr);
00219 if (B && B!=A) {
00220 ierr = MatPermute
00221 (B,local_to_global,local_to_global,&Buse); CHKERRQ(ierr);
00222 } else Buse = Ause;
00223 ierr = ISDestroy(local_to_global); CHKERRQ(ierr);
00224 ierr = ISDestroy(global_to_global); CHKERRQ(ierr);
00225 }
00226 *(VecScatter*)ctx = perm;
00227
00228 {
00229 int local_size; MPI_Comm comm;
00230 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00231 ierr = MatGetLocalSize(Ause,&local_size,PETSC_NULL); CHKERRQ(ierr);
00232 ierr = VecCreateMPI(comm,local_size,PETSC_DECIDE,&Xuse); CHKERRQ(ierr);
00233 ierr = VecDuplicate(Xuse,&Iuse); CHKERRQ(ierr);
00234 ierr = VecDuplicate(Xuse,&Vuse); CHKERRQ(ierr);
00235 }
00236
00237
00238 ierr = VecScatterBegin
00239 (perm,V,Vuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00240 ierr = VecScatterEnd
00241 (perm,V,Vuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00242
00243 ierr = VecScatterBegin
00244 (perm,X,Xuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00245 ierr = VecScatterEnd
00246 (perm,X,Xuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00247
00248 ierr = VecScatterBegin
00249 (perm,I,Iuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00250 ierr = VecScatterEnd
00251 (perm,I,Iuse,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
00252
00253 ierr = LinearSystemSetParts
00254 (newproblem,Ause,Buse,Vuse,Xuse,Iuse); CHKERRQ(ierr);
00255
00256 goto done;
00257 }
00258
00259 SETERRQ1(1,"Undefined distribution <%s>",type);
00260
00261 done:
00262 *outproblem = (NumericalProblem)newproblem;
00263 PetscFunctionReturn(0);
00264 }
00265
00266 #undef __FUNCT__
00267 #define __FUNCT__ "undistribute_system"
00268 static PetscErrorCode undistribute_system
00269 (char *scaling_type,PetscTruth overwrite,void *gctx,void *ctx,
00270 NumericalProblem problem,NumericalProblem nextproblem,
00271 NumericalSolution before,NumericalSolution after)
00272 {
00273 LinearSolution shuffled = (LinearSolution)before,
00274 rectified = (LinearSolution)after;
00275 VecScatter perm = (VecScatter)ctx; Vec Xfore,Xback;
00276 PetscErrorCode ierr;
00277 PetscFunctionBegin;
00278
00279 ierr = LinearSolutionCopyStats(shuffled,rectified); CHKERRQ(ierr);
00280 ierr = LinearSolutionGetVector(shuffled,&Xfore); CHKERRQ(ierr);
00281 ierr = LinearSolutionGetVector(rectified,&Xback); CHKERRQ(ierr);
00282 if (perm) {
00283 ierr = VecScatterBegin
00284 (perm,Xfore,Xback,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00285 ierr = VecScatterEnd
00286 (perm,Xfore,Xback,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00287 ierr = VecScatterDestroy(perm); CHKERRQ(ierr);
00288 } else {
00289 ierr = VecCopy(Xfore,Xback); CHKERRQ(ierr);
00290 }
00291
00292
00293 ierr = DeleteLinearSystem((LinearSystem)problem); CHKERRQ(ierr);
00294 PetscFunctionReturn(0);
00295 }
00296
00297 #undef __FUNCT__
00298 #define __FUNCT__ "DeclareDistributionPreprocessor"
00299 PetscErrorCode DeclareDistributionPreprocessor(void)
00300 {
00301 PetscErrorCode ierr;
00302 PetscFunctionBegin;
00303 ierr = DeclarePreprocessor
00304 ("distribution",
00305 &setup_distribution_choices,&specific_distribution_choices,0,0,
00306 0,0,
00307 &distribute_system,&undistribute_system); CHKERRQ(ierr);
00308 ierr = PreprocessorSetPreservedCategories
00309 (PREPROCESSOR,"laplace"); CHKERRQ(ierr);
00310 PetscFunctionReturn(0);
00311 }