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 { 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