SALSA Analysis Modules
|
00001 /*! \file jpl.c \ingroup categories 00002 \brief Functions for computing the JPL multicolour structure of a matrix 00003 00004 \section jpl Jones-Plassmann multi-colouring 00005 00006 \section Usage 00007 00008 Activate this module with 00009 00010 PetscErrorCode RegisterJPLModules(); 00011 00012 Compute these elements with 00013 00014 ComputeQuantity("jpl","element",A,(void*)&res,&flg); 00015 00016 Available elements are: 00017 - "n-colours" : the number of colours computed 00018 - "colour-set-sizes" : an array containing in location i the number of 00019 points of colour i. (See \ref arraytypes for technicalities on arrays.) 00020 - "colour-offsets" : locations (zero-based) of the colour sets in the 00021 colours array 00022 - "colours" : points sorted first by colour, then increasing index 00023 00024 This routine works in parallel; the stored result is distributed 00025 over all participating processors. 00026 00027 \section References 00028 00029 \verbatim 00030 @article{JoPl:heuristic, 00031 author = {M.T. Jones and P.E. Plassmann}, 00032 title = {A parallel graph coloring heuristic}, 00033 journal = SJSSC, 00034 volume = {14}, 00035 year = {1993}, 00036 keywords = {incomplete factorization, multicolouring, matrix graph} 00037 } 00038 00039 @inproceedings{JoPl:parallelunstructured, 00040 author = {M.T. Jones and P.E. Plassmann}, 00041 title = {Parallel solution of unstructed, sparse systems of linear equations}, 00042 booktitle = {Proceedings of the Sixth SIAM conference on 00043 Parallel Processing for Scientific Computing}, 00044 editor = {R.F. Sincovec and D.E. Keyes and M.R. Leuze and L.R. Petzold 00045 and D.A. Reed}, 00046 publisher = {SIAM}, 00047 address = {Philadelphia}, 00048 pages = {471--475}, 00049 year = {1993} 00050 00051 \section Disclaimer 00052 00053 Unabashed use of British spelling of `colour' needs no disclaimer. 00054 } 00055 00056 \endverbatim 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 Investigate the local and remote connections of a variable 00074 to see whether it needs to be coloured, and if so, with what colour. 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, PetscBool *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 /* loop over all points connected to 'this_var' */ 00091 int other,clr; 00092 other=idx[j]; clr = (int)(clr_array[other]); 00093 00094 /* no need to look at ourselves */ 00095 if (other==compare) continue; 00096 /* this really shouldn't happen */ 00097 /* 00098 if (ran_array[other]==this_rand) 00099 SETERRQ3(MPI_COMM_WORLD,1,"Use a better random: %e @ %d,%d",this_rand,compare,other); 00100 */ 00101 /* 00102 * If the other variable is uncoloured and has a higher random value, 00103 * we can not colour the current variable. 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 * Otherwise, we take note of the colour of the other variable 00112 * for the computation of the colour of this variable later 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 Compute a multicolour ordering of a matrix, in parallel. 00126 00127 This is the main computational routine. 00128 */ 00129 static PetscErrorCode JonesPlassmannColouring(Mat A,PetscBool *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 PetscBool parallel,sequential; PetscErrorCode ierr; 00137 00138 PetscFunctionBegin; 00139 00140 *flg = PETSC_FALSE; 00141 00142 /* 00143 * What is the matrix type? We only handle AIJ (nonblocked) types 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 * Set up the matrices and vectors 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 /* A temporary array for storing neighbour information */ 00168 { 00169 int l1,l2,len; PetscBool 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 * Create the vectors for random and colours 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 * Here is the main loop. 00216 * Iterate until everything is coloured. 00217 */ 00218 ncoloured = 0; max_clr = 0; 00219 for (pass=0; ; pass++) { 00220 PetscBool busy = PETSC_FALSE,gbusy; int var; 00221 /*printf(".. pass %d :",pass);*/ 00222 if (max_clr>pass) 00223 SETERRQ2(MPI_COMM_WORLD,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 * Loop over all variables in every pass 00229 */ 00230 for (var=0; var<local_size; var++) { 00231 int nn = 0; PetscBool now; 00232 00233 /* Skip already coloured variables */ 00234 if (clr_array[var]>0) continue; 00235 00236 /* All others, look at all connections */ 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 /* If this variable has the highest value among uncoloured neighbours, 00249 * colour it now. */ 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 /*printf(" %d->%d",var,clr);*/ 00259 } else { 00260 /*printf(" %d.",var);*/ 00261 } 00262 } /*printf("\n");*/ 00263 00264 /* Exit if all points are coloured */ 00265 ierr = MPI_Allreduce 00266 ((void*)&ncoloured,(void*)&total_coloured,1,MPI_INT,MPI_SUM,comm); 00267 if (total_coloured==total_size) break; 00268 00269 /* If nothing was coloured in this pass, we're stuck */ 00270 ierr = MPI_Allreduce 00271 ((void*)&busy,(void*)&gbusy,1,MPI_INT,MPI_MAX,comm); 00272 if (!gbusy) SETERRQ(MPI_COMM_WORLD,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 * Store the number of colours 00284 */ 00285 ierr = 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 * Count how many points of each colour 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])]++; /* colours are from 1 */ 00303 /* 00304 printf("coloured:"); 00305 for (i=0; i<local_size; i++) printf(" %e",clr_array[i]); 00306 printf("\n"); 00307 printf("counts:"); 00308 for (i=0; i<colour_count[0]; i++) printf(" %d",colour_count[i+1]); 00309 printf("\n"); 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 printf("offsets:"); 00323 for (i=0; i<colour_offsets[0]; i++) printf(" %d",colour_offsets[i+1]); 00324 printf("\n"); 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 printf("colours:"); 00341 for (i=0; i<colours[0]; i++) printf(" %d",colours[i+1]); 00342 printf("\n"); 00343 */ 00344 ierr = GetDataID("jpl","colours",&id,PETSC_NULL); CHKERRQ(ierr); 00345 ierr = PetscObjectComposedDataSetIntstar 00346 ((PetscObject)A,id,colours); CHKERRQ(ierr); 00347 } 00348 00349 /* Clean up */ 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 /* clr_bord is the aij->lvec array; not to be destroyed */ 00362 } 00363 00364 *flg = PETSC_TRUE; 00365 00366 00367 PetscFunctionReturn(0); 00368 } 00369 00370 #undef __FUNCT__ 00371 #define __FUNCT__ "NColours" 00372 /*! 00373 Compute the global number of colours. Note that any given processor does 00374 not need to have all the colours. 00375 */ 00376 static PetscErrorCode NColours 00377 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00378 { 00379 Mat A = (Mat)prob; 00380 int v = 0; PetscBool 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 Compute an array that has the number of points of each colour, 00403 on this processor. If computing the colouring is done in parallel, 00404 then the multicolouring is stored in parallel: it is not gathered onto 00405 any processor. 00406 */ 00407 static PetscErrorCode ColourSizes 00408 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00409 { 00410 Mat A = (Mat)prob; 00411 int llv=0,*v = NULL; PetscBool 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 Compute an array that has the offsets of the locations of the colours, 00446 on this processor. If computing the colouring is done in parallel, 00447 then the multicolouring is stored in parallel: it is not gathered onto 00448 any processor. 00449 */ 00450 static PetscErrorCode ColourOffsets 00451 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00452 { 00453 Mat A = (Mat)prob; 00454 int llv=0,*v = NULL; PetscBool 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 Compute an array that has the colour sets on this processor. In order 00489 to know where a set starts or ends, the result of ColourOffsets() is 00490 needed. 00491 00492 If computing the colouring is done in parallel, 00493 then the multicolouring is stored in parallel: it is not gathered onto 00494 any processor. 00495 */ 00496 static PetscErrorCode Colours 00497 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00498 { 00499 Mat A = (Mat)prob; int N; 00500 int *v = NULL; PetscBool 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 }