00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058 #include <stdlib.h>
00059 #include "anamod.h"
00060 #include "anamodsalsamodules.h"
00061 #include "anamatrix.h"
00062 #include "petscerror.h"
00063 #include "petscmat.h"
00064
00065 #include "petscis.h"
00066 #include "petscmat.h"
00067 #include "src/mat/impls/aij/mpi/mpiaij.h"
00068 #include "src/mat/impls/aij/seq/aij.h"
00069
00070 #undef __FUNCT__
00071 #define __FUNCT__ "LookAtUncolouredVar"
00072
00073
00074
00075
00076 static PetscErrorCode LookAtUncolouredVar
00077 (int compare,int this_var,Mat A,
00078 PetscScalar *clr_array,PetscScalar *ran_array,PetscReal this_rand,
00079 int *neighb,int *nneigh, PetscTruth *colour_now)
00080 {
00081
00082 int nnz,j; const int *idx; PetscErrorCode ierr;
00083
00084 PetscFunctionBegin;
00085
00086 *colour_now = PETSC_TRUE;
00087
00088 ierr = MatGetRow(A,this_var,&nnz,&idx,PETSC_NULL); CHKERRQ(ierr);
00089 for (j=0; j<nnz; j++) {
00090
00091 int other,clr;
00092 other=idx[j]; clr = (int)(clr_array[other]);
00093
00094
00095 if (other==compare) continue;
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 if ( (clr==0 && ran_array[other]>this_rand) ||
00106 (ran_array[other]==this_rand && other>this_var) ) {
00107 *colour_now = PETSC_FALSE;
00108 goto exit;
00109 }
00110
00111
00112
00113
00114 neighb[(*nneigh)++] = clr;
00115 }
00116 exit:
00117 ierr = MatRestoreRow(A,this_var,&nnz,&idx,PETSC_NULL); CHKERRQ(ierr);
00118
00119 PetscFunctionReturn(0);
00120 }
00121
00122 #undef __FUNCT__
00123 #define __FUNCT__ "JonesPlassmannColouring"
00124
00125
00126
00127
00128
00129 static PetscErrorCode JonesPlassmannColouring(Mat A,PetscTruth *flg)
00130 {
00131 MPI_Comm comm; const MatType type;
00132 Mat dia,off; Vec rand,ran_bord, clrs,clrs_bkp,clr_bord; VecScatter sctr;
00133 PetscScalar *ran_array,*bran_array,*clr_array,*clr_bkp_array,*bclr_array;
00134 int *neighb,ncoloured,total_coloured,pass,max_clr,
00135 global_maxclr,local_size,total_size, id;
00136 PetscTruth parallel,sequential; PetscErrorCode ierr;
00137
00138 PetscFunctionBegin;
00139
00140 *flg = PETSC_FALSE;
00141
00142
00143
00144
00145 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00146 ierr = MatGetType(A,&type); CHKERRQ(ierr);
00147 ierr = PetscStrcmp(type,MATSEQAIJ,&sequential); CHKERRQ(ierr);
00148 ierr = PetscStrcmp(type,MATMPIAIJ,¶llel); CHKERRQ(ierr);
00149 if (!(sequential || parallel)) {
00150 *flg = PETSC_FALSE;
00151 PetscFunctionReturn(0);
00152 }
00153
00154
00155
00156
00157 if (parallel) {
00158 Mat_MPIAIJ *aij = (Mat_MPIAIJ *) A->data;
00159 dia=aij->A; off=aij->B; clr_bord = aij->lvec; sctr = aij->Mvctx;
00160 ierr = VecDuplicate(aij->lvec,&ran_bord); CHKERRQ(ierr);
00161 } else {
00162 dia = A; off = NULL; clr_bord = NULL; sctr = NULL;
00163 }
00164 ierr = MatGetSize(A,&total_size,PETSC_NULL); CHKERRQ(ierr);
00165 ierr = MatGetLocalSize(A,&local_size,PETSC_NULL); CHKERRQ(ierr);
00166
00167
00168 {
00169 int l1,l2,len; PetscTruth f;
00170 ierr = MatrixComputeQuantity
00171 (dia,"structure","max-nnzeros-per-row",(AnalysisItem*)&l1,PETSC_NULL,
00172 &f); CHKERRQ(ierr);
00173 if (parallel) {
00174 ierr = MatrixComputeQuantity
00175 (off,"structure","max-nnzeros-per-row",(AnalysisItem*)&l2,PETSC_NULL,
00176 &f); CHKERRQ(ierr);
00177 if (f) len=l1+l2; else len = local_size;
00178 } else len = l1;
00179 ierr = PetscMalloc((len+1)*sizeof(int),&neighb); CHKERRQ(ierr);
00180 }
00181
00182
00183
00184
00185 {
00186 if (parallel) {
00187 ierr = VecCreateMPI(comm,local_size,PETSC_DECIDE,&rand); CHKERRQ(ierr);
00188 } else {
00189 ierr = VecCreateSeq(MPI_COMM_SELF,local_size,&rand); CHKERRQ(ierr);
00190 }
00191 {
00192 PetscRandom rctx;
00193 ierr = PetscRandomCreate(MPI_COMM_SELF,&rctx); CHKERRQ(ierr);
00194 ierr = PetscRandomSetType(rctx,PETSCRAND48); CHKERRQ(ierr);
00195 ierr = VecSetRandom(rand,rctx); CHKERRQ(ierr);
00196 ierr = PetscRandomDestroy(rctx); CHKERRQ(ierr);
00197 }
00198 ierr = VecDuplicate(rand,&clrs); CHKERRQ(ierr);
00199 ierr = VecDuplicate(rand,&clrs_bkp); CHKERRQ(ierr);
00200 ierr = VecSet(clrs,0.); CHKERRQ(ierr);
00201 if (parallel) {
00202 ierr = VecScatterBegin
00203 (sctr,rand,ran_bord,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00204 ierr = VecScatterEnd
00205 (sctr,rand,ran_bord,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00206 ierr = VecGetArray(ran_bord,&bran_array); CHKERRQ(ierr);
00207 ierr = VecGetArray(clr_bord,&bclr_array); CHKERRQ(ierr);
00208 }
00209 ierr = VecGetArray(rand,&ran_array); CHKERRQ(ierr);
00210 ierr = VecGetArray(clrs,&clr_array); CHKERRQ(ierr);
00211 ierr = VecGetArray(clrs_bkp,&clr_bkp_array); CHKERRQ(ierr);
00212 }
00213
00214
00215
00216
00217
00218 ncoloured = 0; max_clr = 0;
00219 for (pass=0; ; pass++) {
00220 PetscTruth busy = PETSC_FALSE,gbusy; int var;
00221
00222 if (max_clr>pass)
00223 SETERRQ2(1,"This can not happen: maxclr=%d in pass %d",max_clr,pass);
00224 ierr = PetscMemcpy
00225 (clr_bkp_array,clr_array,local_size*sizeof(PetscScalar)); CHKERRQ(ierr);
00226
00227
00228
00229
00230 for (var=0; var<local_size; var++) {
00231 int nn = 0; PetscTruth now;
00232
00233
00234 if (clr_array[var]>0) continue;
00235
00236
00237 ierr = LookAtUncolouredVar
00238 (var,var,dia,clr_bkp_array,ran_array,ran_array[var],
00239 neighb,&nn,&now); CHKERRQ(ierr);
00240 if (now && parallel) {
00241 if (parallel) {
00242 ierr = LookAtUncolouredVar
00243 (-1,var,off,bclr_array,bran_array,ran_array[var],
00244 neighb,&nn,&now); CHKERRQ(ierr);
00245 }
00246 }
00247
00248
00249
00250 if (now) {
00251 int clr=1,i;
00252 attempt:
00253 for (i=0; i<nn; i++)
00254 if (neighb[i]==clr) {clr++; goto attempt;}
00255 clr_array[var] = (PetscScalar) clr;
00256 if (clr>max_clr) max_clr = clr;
00257 ncoloured++; busy = PETSC_TRUE;
00258
00259 } else {
00260
00261 }
00262 }
00263
00264
00265 MPI_Allreduce
00266 ((void*)&ncoloured,(void*)&total_coloured,1,MPI_INT,MPI_SUM,comm);
00267 if (total_coloured==total_size) break;
00268
00269
00270 MPI_Allreduce
00271 ((void*)&busy,(void*)&gbusy,1,MPI_INT,MPI_MAX,comm);
00272 if (!gbusy) SETERRQ(1,"We are stuck");
00273
00274 if (parallel) {
00275 ierr = VecScatterBegin
00276 (sctr,clrs,clr_bord,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00277 ierr = VecScatterEnd
00278 (sctr,clrs,clr_bord,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
00279 }
00280 }
00281
00282
00283
00284
00285 MPI_Allreduce
00286 ((void*)&max_clr,(void*)&global_maxclr,1,MPI_INT,MPI_MAX,comm);
00287 ierr = GetDataID("jpl","n-colours",&id,PETSC_NULL); CHKERRQ(ierr);
00288 ierr = PetscObjectComposedDataSetInt
00289 ((PetscObject)A,id,global_maxclr); CHKERRQ(ierr);
00290
00291
00292
00293
00294 {
00295 int *colour_count,*colour_offsets,*ptr,*colours,clr,i;
00296 ierr = PetscMalloc
00297 ((global_maxclr+1)*sizeof(int),&colour_count); CHKERRQ(ierr);
00298 colour_count[0] = global_maxclr;
00299 for (clr=0; clr<global_maxclr; clr++)
00300 colour_count[clr+1] = 0;
00301 for (i=0; i<local_size; i++)
00302 colour_count[(int)(clr_array[i])]++;
00303
00304
00305
00306
00307
00308
00309
00310
00311 ierr = GetDataID("jpl","colour-set-sizes",&id,PETSC_NULL); CHKERRQ(ierr);
00312 ierr = PetscObjectComposedDataSetIntstar
00313 ((PetscObject)A,id,colour_count); CHKERRQ(ierr);
00314
00315 ierr = PetscMalloc
00316 ((global_maxclr+1)*sizeof(int),&colour_offsets); CHKERRQ(ierr);
00317 colour_offsets[0] = global_maxclr;
00318 colour_offsets[1] = 0;
00319 for (clr=1; clr<global_maxclr; clr++)
00320 colour_offsets[1+clr] = colour_offsets[clr]+colour_count[clr];
00321
00322
00323
00324
00325
00326 ierr = GetDataID("jpl","colour-offsets",&id,PETSC_NULL); CHKERRQ(ierr);
00327 ierr = PetscObjectComposedDataSetIntstar
00328 ((PetscObject)A,id,colour_offsets); CHKERRQ(ierr);
00329
00330 ierr = PetscMalloc((local_size+1)*sizeof(int),&colours); CHKERRQ(ierr);
00331 colours[0] = local_size;
00332 ierr = PetscMalloc(global_maxclr*sizeof(int),&ptr); CHKERRQ(ierr);
00333 PetscMemcpy(ptr,colour_offsets+1,global_maxclr*sizeof(int));
00334 for (i=0; i<local_size; i++) {
00335 int c=(int)(clr_array[i]);
00336 colours[1+ptr[c-1]++] = i;
00337 }
00338 ierr = PetscFree(ptr); CHKERRQ(ierr);
00339
00340
00341
00342
00343
00344 ierr = GetDataID("jpl","colours",&id,PETSC_NULL); CHKERRQ(ierr);
00345 ierr = PetscObjectComposedDataSetIntstar
00346 ((PetscObject)A,id,colours); CHKERRQ(ierr);
00347 }
00348
00349
00350 ierr = PetscFree(neighb); CHKERRQ(ierr);
00351 ierr = VecRestoreArray(rand,&ran_array); CHKERRQ(ierr);
00352 ierr = VecDestroy(rand); CHKERRQ(ierr);
00353 ierr = VecRestoreArray(clrs,&clr_array); CHKERRQ(ierr);
00354 ierr = VecDestroy(clrs); CHKERRQ(ierr);
00355 ierr = VecRestoreArray(clrs_bkp,&clr_bkp_array); CHKERRQ(ierr);
00356 ierr = VecDestroy(clrs_bkp); CHKERRQ(ierr);
00357 if (parallel) {
00358 ierr = VecRestoreArray(ran_bord,&bran_array); CHKERRQ(ierr);
00359 ierr = VecDestroy(ran_bord); CHKERRQ(ierr);
00360 ierr = VecRestoreArray(clr_bord,&bclr_array); CHKERRQ(ierr);
00361
00362 }
00363
00364 *flg = PETSC_TRUE;
00365
00366
00367 PetscFunctionReturn(0);
00368 }
00369
00370 #undef __FUNCT__
00371 #define __FUNCT__ "NColours"
00372
00373
00374
00375
00376 static PetscErrorCode NColours
00377 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00378 {
00379 Mat A = (Mat)prob;
00380 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00381 PetscFunctionBegin;
00382 ierr = GetDataID("jpl","n-colours",&id,&has); CHKERRQ(ierr);
00383 HASTOEXIST(has);
00384 ierr = PetscObjectComposedDataGetInt
00385 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00386 if (*flg) rv->i = v;
00387 else {
00388 ierr = JonesPlassmannColouring(A,flg); CHKERRQ(ierr);
00389 if (*flg) {
00390 ierr = PetscObjectComposedDataGetInt
00391 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00392 if (*flg) rv->i = v;
00393 }
00394 }
00395
00396 PetscFunctionReturn(0);
00397 }
00398
00399 #undef __FUNCT__
00400 #define __FUNCT__ "ColourSizes"
00401
00402
00403
00404
00405
00406
00407 static PetscErrorCode ColourSizes
00408 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00409 {
00410 Mat A = (Mat)prob;
00411 int llv=0,*v = NULL; PetscTruth has; int id,id2; PetscErrorCode ierr;
00412 PetscFunctionBegin;
00413
00414 ierr = GetDataID("jpl","n-colours",&id2,&has); CHKERRQ(ierr);
00415 HASTOEXIST(has);
00416 ierr = PetscObjectComposedDataGetInt
00417 ((PetscObject)A,id2,llv,*flg); CHKERRQ(ierr);
00418 if (*flg) {
00419 if (lv) *lv = llv;
00420 ierr = GetDataID("jpl","colour-set-sizes",&id,&has); CHKERRQ(ierr);
00421 HASTOEXIST(has);
00422 ierr = PetscObjectComposedDataGetIntstar
00423 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00424 }
00425 if (*flg) rv->ii = v;
00426 else {
00427 ierr = JonesPlassmannColouring(A,flg); CHKERRQ(ierr);
00428 if (*flg) {
00429 ierr = PetscObjectComposedDataGetIntstar
00430 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00431 if (*flg) rv->ii = v;
00432 ierr = PetscObjectComposedDataGetInt
00433 ((PetscObject)A,id2,llv,*flg); CHKERRQ(ierr);
00434 if (*flg)
00435 if (lv) *lv = llv;
00436 }
00437 }
00438
00439 PetscFunctionReturn(0);
00440 }
00441
00442 #undef __FUNCT__
00443 #define __FUNCT__ "ColourOffsets"
00444
00445
00446
00447
00448
00449
00450 static PetscErrorCode ColourOffsets
00451 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00452 {
00453 Mat A = (Mat)prob;
00454 int llv=0,*v = NULL; PetscTruth has; int id,id2; PetscErrorCode ierr;
00455 PetscFunctionBegin;
00456
00457 ierr = GetDataID("jpl","n-colours",&id2,&has); CHKERRQ(ierr);
00458 HASTOEXIST(has);
00459 ierr = PetscObjectComposedDataGetInt
00460 ((PetscObject)A,id2,llv,*flg); CHKERRQ(ierr);
00461 if (*flg) {
00462 if (lv) *lv = llv;
00463 ierr = GetDataID("jpl","colour-offsets",&id,&has); CHKERRQ(ierr);
00464 HASTOEXIST(has);
00465 ierr = PetscObjectComposedDataGetIntstar
00466 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00467 }
00468 if (*flg) rv->ii = v;
00469 else {
00470 ierr = JonesPlassmannColouring(A,flg); CHKERRQ(ierr);
00471 if (*flg) {
00472 ierr = PetscObjectComposedDataGetIntstar
00473 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00474 if (*flg) rv->ii = v;
00475 ierr = PetscObjectComposedDataGetInt
00476 ((PetscObject)A,id2,llv,*flg); CHKERRQ(ierr);
00477 if (*flg)
00478 if (lv) *lv = llv;
00479 }
00480 }
00481
00482 PetscFunctionReturn(0);
00483 }
00484
00485 #undef __FUNCT__
00486 #define __FUNCT__ "Colours"
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496 static PetscErrorCode Colours
00497 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00498 {
00499 Mat A = (Mat)prob; int N;
00500 int *v = NULL; PetscTruth has; int id; PetscErrorCode ierr;
00501 PetscFunctionBegin;
00502 ierr = MatGetLocalSize(A,&N,PETSC_NULL); CHKERRQ(ierr);
00503 if (lv) *lv = N;
00504 ierr = GetDataID("jpl","colours",&id,&has); CHKERRQ(ierr);
00505 HASTOEXIST(has);
00506 ierr = PetscObjectComposedDataGetIntstar
00507 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00508 if (*flg) rv->ii = v;
00509 else {
00510 ierr = JonesPlassmannColouring(A,flg); CHKERRQ(ierr);
00511 if (*flg) {
00512 ierr = PetscObjectComposedDataGetIntstar
00513 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00514 if (*flg) rv->ii = v;
00515 }
00516 }
00517
00518 PetscFunctionReturn(0);
00519 }
00520
00521 #undef __FUNCT__
00522 #define __FUNCT__ "RegisterJPLModules"
00523 PetscErrorCode RegisterJPLModules(void)
00524 {
00525 PetscErrorCode ierr;
00526 PetscFunctionBegin;
00527
00528 ierr = RegisterModule
00529 ("jpl","n-colours",ANALYSISINTEGER,&NColours); CHKERRQ(ierr);
00530 ierr = RegisterModule
00531 ("jpl","colour-set-sizes",ANALYSISINTARRAY,&ColourSizes); CHKERRQ(ierr);
00532 ierr = RegisterModule
00533 ("jpl","colour-offsets",ANALYSISINTARRAY,&ColourOffsets); CHKERRQ(ierr);
00534 ierr = RegisterModule
00535 ("jpl","colours",ANALYSISINTARRAY,&Colours); CHKERRQ(ierr);
00536
00537 PetscFunctionReturn(0);
00538 }