1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2 #include <petsc/private/isimpl.h>
3
MatColoringApply_Natural(MatColoring mc,ISColoring * iscoloring)4 static PetscErrorCode MatColoringApply_Natural(MatColoring mc, ISColoring *iscoloring)
5 {
6 PetscInt start, end, i, bs = 1, n;
7 ISColoringValue *colors;
8 MPI_Comm comm;
9 PetscBool flg1, flg2;
10 Mat mat = mc->mat;
11 Mat mat_seq = mc->mat;
12 PetscMPIInt size;
13 ISColoring iscoloring_seq;
14 ISColoringValue *colors_loc;
15 PetscInt rstart, rend, N_loc, nc;
16
17 PetscFunctionBegin;
18 /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
19 PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1));
20 PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2));
21 if (flg1 || flg2) PetscCall(MatGetBlockSize(mat, &bs));
22
23 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
24 PetscCallMPI(MPI_Comm_size(comm, &size));
25 if (size > 1) {
26 /* create a sequential iscoloring on all processors */
27 PetscCall(MatGetSeqNonzeroStructure(mat, &mat_seq));
28 }
29
30 PetscCall(MatGetSize(mat_seq, &n, NULL));
31 PetscCall(MatGetOwnershipRange(mat_seq, &start, &end));
32 n = n / bs;
33 PetscCheck(n <= IS_COLORING_MAX - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Maximum color size exceeded");
34
35 start = start / bs;
36 end = end / bs;
37 PetscCall(PetscMalloc1(end - start + 1, &colors));
38 for (i = start; i < end; i++) colors[i - start] = (ISColoringValue)i;
39 PetscCall(ISColoringCreate(comm, n, end - start, colors, PETSC_OWN_POINTER, iscoloring));
40
41 if (size > 1) {
42 PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
43
44 /* convert iscoloring_seq to a parallel iscoloring */
45 iscoloring_seq = *iscoloring;
46 rstart = mat->rmap->rstart / bs;
47 rend = mat->rmap->rend / bs;
48 N_loc = rend - rstart; /* number of local nodes */
49
50 /* get local colors for each local node */
51 PetscCall(PetscMalloc1(N_loc + 1, &colors_loc));
52 for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
53 /* create a parallel iscoloring */
54 nc = iscoloring_seq->n;
55 PetscCall(ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring));
56 PetscCall(ISColoringDestroy(&iscoloring_seq));
57 }
58 PetscFunctionReturn(PETSC_SUCCESS);
59 }
60
61 /*MC
62 MATCOLORINGNATURAL - implements a trivial coloring routine with one color per column
63
64 Level: beginner
65
66 Note:
67 Using this coloring would be extremely inefficient but it is useful for testing
68
69 .seealso: `MatColoring`, `MatColoringType`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MatColoringType`
70 M*/
MatColoringCreate_Natural(MatColoring mc)71 PETSC_EXTERN PetscErrorCode MatColoringCreate_Natural(MatColoring mc)
72 {
73 PetscFunctionBegin;
74 mc->data = NULL;
75 mc->ops->apply = MatColoringApply_Natural;
76 mc->ops->view = NULL;
77 mc->ops->destroy = NULL;
78 mc->ops->setfromoptions = NULL;
79 PetscFunctionReturn(PETSC_SUCCESS);
80 }
81