1: static char help[] = "Tests PCMGGetRScale with a square interpolation matrix and different row and column layouts.\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **argv)
  6: {
  7:   const PetscScalar *values;
  8:   Mat                P;
  9:   PC                 pc;
 10:   Vec                rscale;
 11:   PetscInt           localRows, localColumns, start, end;
 12:   PetscMPIInt        rank;

 14:   PetscFunctionBeginUser;
 15:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 16:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 17:   localRows    = rank ? 3 : 1;
 18:   localColumns = rank ? 1 : 3;
 19:   PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, localRows, localColumns, 4, 4, 4, NULL, 4, NULL, &P));
 20:   PetscCall(MatGetOwnershipRange(P, &start, &end));
 21:   for (PetscInt row = start; row < end; ++row) {
 22:     for (PetscInt col = 0; col < 4; ++col) PetscCall(MatSetValue(P, row, col, row + col + 1, INSERT_VALUES));
 23:   }
 24:   PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
 25:   PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));

 27:   PetscCall(PCCreate(PETSC_COMM_WORLD, &pc));
 28:   PetscCall(PCSetType(pc, PCMG));
 29:   PetscCall(PCMGSetLevels(pc, 2, NULL));
 30:   PetscCall(PCMGSetInterpolation(pc, 1, P));
 31:   PetscCall(PCMGGetRScale(pc, 1, &rscale));
 32:   PetscCall(VecGetOwnershipRange(rscale, &start, &end));
 33:   PetscCall(VecGetArrayRead(rscale, &values));
 34:   for (PetscInt i = 0; i < end - start; ++i) {
 35:     const PetscInt    col      = start + i;
 36:     const PetscScalar expected = 1.0 / (PetscReal)(10 + 4 * col);

 38:     PetscCheck(PetscAbsScalar(values[i] - expected) < 100.0 * PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect restriction scaling at column %" PetscInt_FMT, col);
 39:   }
 40:   PetscCall(VecRestoreArrayRead(rscale, &values));

 42:   PetscCall(PCDestroy(&pc));
 43:   PetscCall(MatDestroy(&P));
 44:   PetscCall(PetscFinalize());
 45:   return 0;
 46: }

 48: /*TEST

 50:   test:
 51:     nsize: 2
 52:     output_file: output/empty.out

 54: TEST*/