Actual source code: ex240.c
1: static char help[] = "Tests MatFDColoringSetValues()\n\n";
3: #include <petscdm.h>
4: #include <petscdmda.h>
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;
23: PetscInitialize(&argc, &argv, NULL, help);
24: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
25: 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: DMSetUp(da);
27: DMCreateMatrix(da, &A);
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: DMCreateColoring(da, IS_COLORING_LOCAL, &iscoloring);
32: ISColoringGetColors(iscoloring, &n, &nc, &icolors);
33: DMGetLocalToGlobalMapping(da, <og);
34: PetscMalloc1(n, &map);
35: for (i = 0; i < n; i++) map[i] = i;
36: ISLocalToGlobalMappingApply(ltog, n, map, map);
37: MatGetSize(A, &N, NULL);
38: 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: PetscFree(map);
42: PetscSynchronizedPrintf(MPI_COMM_WORLD, "[%d]Global colors \n", rank);
43: for (i = 0; i < N; i++) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " %" PetscInt_FMT "\n", i, colors[i]);
44: PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout);
46: /* second, compress the A matrix */
47: MatSetRandom(A, NULL);
48: MatView(A, NULL);
49: MatGetLocalSize(A, &nrow, NULL);
50: MatGetOwnershipRange(A, &rstart, NULL);
51: PetscCalloc1(nrow * nc, &cm);
52: for (i = 0; i < nrow; i++) {
53: MatGetRow(A, rstart + i, &ncols, &cols, &vals);
54: for (j = 0; j < ncols; j++) {
56: cm[i + nrow * colors[cols[j]]] = vals[j];
57: }
58: MatRestoreRow(A, rstart + i, &ncols, &cols, &vals);
59: }
61: /* print compressed matrix */
62: PetscSynchronizedPrintf(MPI_COMM_WORLD, "[%d] Compressed matrix \n", rank);
63: for (i = 0; i < nrow; i++) {
64: for (j = 0; j < nc; j++) PetscSynchronizedPrintf(MPI_COMM_WORLD, "%12.4e ", (double)PetscAbsScalar(cm[i + nrow * j]));
65: PetscSynchronizedPrintf(MPI_COMM_WORLD, "\n");
66: }
67: PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout);
69: /* put the compressed matrix into the standard matrix */
70: MatDuplicate(A, MAT_COPY_VALUES, &B);
71: MatZeroEntries(A);
72: MatView(B, 0);
73: MatFDColoringCreate(A, iscoloring, &fdcoloring);
74: PetscOptionsHasName(NULL, NULL, "-single_block", &single);
75: if (single) MatFDColoringSetBlockSize(fdcoloring, PETSC_DEFAULT, nc);
76: PetscOptionsHasName(NULL, NULL, "-two_block", &two);
77: if (two) MatFDColoringSetBlockSize(fdcoloring, PETSC_DEFAULT, 2);
78: MatFDColoringSetFromOptions(fdcoloring);
79: MatFDColoringSetUp(A, iscoloring, fdcoloring);
81: MatFDColoringSetValues(A, fdcoloring, cm);
82: MatView(A, NULL);
84: /* check the values were put in the correct locations */
85: MatAXPY(A, -1.0, B, SAME_NONZERO_PATTERN);
86: MatView(A, NULL);
87: MatNorm(A, NORM_FROBENIUS, &norm);
88: if (norm > PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Matrix is not identical, problem with MatFDColoringSetValues()\n");
89: PetscFree(colors);
90: PetscFree(cm);
91: ISColoringDestroy(&iscoloring);
92: MatFDColoringDestroy(&fdcoloring);
93: MatDestroy(&A);
94: MatDestroy(&B);
95: DMDestroy(&da);
96: PetscFinalize();
97: return 0;
98: }
100: /*TEST
102: test:
103: nsize: 2
104: requires: !complex
106: test:
107: suffix: single
108: requires: !complex
109: nsize: 2
110: args: -single_block
111: output_file: output/ex240_1.out
113: test:
114: suffix: two
115: requires: !complex
116: nsize: 2
117: args: -two_block
118: output_file: output/ex240_1.out
120: test:
121: suffix: 2
122: requires: !complex
123: nsize: 5
125: TEST*/