1 static char help[] = "Tests MatFDColoringSetValues()\n\n";
2
3 #include <petscdm.h>
4 #include <petscdmda.h>
5
main(int argc,char ** argv)6 int main(int argc, char **argv)
7 {
8 DM da;
9 PetscInt N, mx = 5, my = 4, i, j, nc, nrow, n, ncols, rstart, *colors, *map;
10 const PetscInt *cols;
11 const PetscScalar *vals;
12 Mat A, B;
13 PetscReal norm;
14 MatFDColoring fdcoloring;
15 ISColoring iscoloring;
16 PetscScalar *cm;
17 const ISColoringValue *icolors;
18 PetscMPIInt rank;
19 ISLocalToGlobalMapping ltog;
20 PetscBool single, two;
21
22 PetscFunctionBeginUser;
23 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
24 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
25 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, mx, my, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
26 PetscCall(DMSetUp(da));
27 PetscCall(DMCreateMatrix(da, &A));
28
29 /* as a test copy the matrices from the standard format to the compressed format; this is not scalable but is ok because just for testing */
30 /* first put the coloring in the global ordering */
31 PetscCall(DMCreateColoring(da, IS_COLORING_LOCAL, &iscoloring));
32 PetscCall(ISColoringGetColors(iscoloring, &n, &nc, &icolors));
33 PetscCall(DMGetLocalToGlobalMapping(da, <og));
34 PetscCall(PetscMalloc1(n, &map));
35 for (i = 0; i < n; i++) map[i] = i;
36 PetscCall(ISLocalToGlobalMappingApply(ltog, n, map, map));
37 PetscCall(MatGetSize(A, &N, NULL));
38 PetscCall(PetscMalloc1(N, &colors));
39 for (i = 0; i < N; i++) colors[i] = -1;
40 for (i = 0; i < n; i++) colors[map[i]] = icolors[i];
41 PetscCall(PetscFree(map));
42 PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "[%d]Global colors \n", rank));
43 for (i = 0; i < N; i++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " %" PetscInt_FMT "\n", i, colors[i]));
44 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
45
46 /* second, compress the A matrix */
47 PetscCall(MatSetRandom(A, NULL));
48 PetscCall(MatView(A, NULL));
49 PetscCall(MatGetLocalSize(A, &nrow, NULL));
50 PetscCall(MatGetOwnershipRange(A, &rstart, NULL));
51 PetscCall(PetscCalloc1(nrow * nc, &cm));
52 for (i = 0; i < nrow; i++) {
53 PetscCall(MatGetRow(A, rstart + i, &ncols, &cols, &vals));
54 for (j = 0; j < ncols; j++) {
55 PetscCheck(colors[cols[j]] >= 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Global column %" PetscInt_FMT " had no color", cols[j]);
56 cm[i + nrow * colors[cols[j]]] = vals[j];
57 }
58 PetscCall(MatRestoreRow(A, rstart + i, &ncols, &cols, &vals));
59 }
60
61 /* print compressed matrix */
62 PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "[%d] Compressed matrix \n", rank));
63 for (i = 0; i < nrow; i++) {
64 for (j = 0; j < nc; j++) PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "%12.4e ", (double)PetscAbsScalar(cm[i + nrow * j])));
65 PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "\n"));
66 }
67 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
68
69 /* put the compressed matrix into the standard matrix */
70 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
71 PetscCall(MatZeroEntries(A));
72 PetscCall(MatView(B, 0));
73 PetscCall(MatFDColoringCreate(A, iscoloring, &fdcoloring));
74 PetscCall(PetscOptionsHasName(NULL, NULL, "-single_block", &single));
75 if (single) PetscCall(MatFDColoringSetBlockSize(fdcoloring, PETSC_DEFAULT, nc));
76 PetscCall(PetscOptionsHasName(NULL, NULL, "-two_block", &two));
77 if (two) PetscCall(MatFDColoringSetBlockSize(fdcoloring, PETSC_DEFAULT, 2));
78 PetscCall(MatFDColoringSetFromOptions(fdcoloring));
79 PetscCall(MatFDColoringSetUp(A, iscoloring, fdcoloring));
80
81 PetscCall(MatFDColoringSetValues(A, fdcoloring, cm));
82 PetscCall(MatView(A, NULL));
83
84 /* check the values were put in the correct locations */
85 PetscCall(MatAXPY(A, -1.0, B, SAME_NONZERO_PATTERN));
86 PetscCall(MatView(A, NULL));
87 PetscCall(MatNorm(A, NORM_FROBENIUS, &norm));
88 if (norm > PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Matrix is not identical, problem with MatFDColoringSetValues()\n"));
89 PetscCall(PetscFree(colors));
90 PetscCall(PetscFree(cm));
91 PetscCall(ISColoringDestroy(&iscoloring));
92 PetscCall(MatFDColoringDestroy(&fdcoloring));
93 PetscCall(MatDestroy(&A));
94 PetscCall(MatDestroy(&B));
95 PetscCall(DMDestroy(&da));
96 PetscCall(PetscFinalize());
97 return 0;
98 }
99
100 /*TEST
101
102 test:
103 nsize: 2
104 requires: !complex
105
106 test:
107 suffix: single
108 requires: !complex
109 nsize: 2
110 args: -single_block
111 output_file: output/ex240_1.out
112
113 test:
114 suffix: two
115 requires: !complex
116 nsize: 2
117 args: -two_block
118 output_file: output/ex240_1.out
119
120 test:
121 suffix: 2
122 requires: !complex
123 nsize: 5
124
125 TEST*/
126