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;

 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));

 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, &ltog));
 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));

 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:   }

 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));

 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));

 81:   PetscCall(MatFDColoringSetValues(A, fdcoloring, cm));
 82:   PetscCall(MatView(A, NULL));

 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: }

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*/