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