1: static char help[] = "Losing nullspaces in PCFIELDSPLIT after zeroing rows.\n";

  3: // Contributed by Jeremy Theler

  5: #include <petscksp.h>

  7: int main(int argc, char **args)
  8: {
  9:   KSP             ksp, *sub_ksp;
 10:   Vec             x, b, rigid_mode[6];
 11:   PetscViewer     viewer;
 12:   PetscInt        rows, cols, size, bs, n_splits = 0;
 13:   PetscBool       has_columns = PETSC_FALSE;
 14:   Mat             A, K;
 15:   MatNullSpace    nullsp, near_null_space;
 16:   IS              is_thermal, is_mech;
 17:   PC              pc, pc_thermal, pc_mech;
 18:   const PetscInt *bc_thermal_indexes, *bc_mech_indexes;
 19:   char            datafilespath[PETSC_MAX_PATH_LEN], datafile[PETSC_MAX_PATH_LEN];

 21:   PetscFunctionBeginUser;
 22:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 23:   PetscCall(PetscOptionsGetString(NULL, NULL, "-datafilespath", datafilespath, sizeof(datafilespath), NULL));

 25:   PetscCall(PetscStrcpy(datafile, datafilespath));
 26:   PetscCall(PetscStrcat(datafile, "/lostnullspace/"));
 27:   PetscCall(PetscStrcat(datafile, "A.bin"));
 28:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, datafile, FILE_MODE_READ, &viewer));
 29:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 30:   PetscCall(MatLoad(A, viewer));
 31:   PetscCall(PetscViewerDestroy(&viewer));

 33:   PetscCall(MatGetNearNullSpace(A, &nullsp));
 34:   PetscCall(MatGetSize(A, &rows, &cols));
 35:   PetscCall(MatGetBlockSize(A, &bs));
 36:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "A has rows = %" PetscInt_FMT ", cols = %" PetscInt_FMT ", bs = %" PetscInt_FMT ", nearnullsp = %p\n", rows, cols, bs, (void *)nullsp));

 38:   PetscCall(PetscStrcpy(datafile, datafilespath));
 39:   PetscCall(PetscStrcat(datafile, "/lostnullspace/"));
 40:   PetscCall(PetscStrcat(datafile, "is.bin"));
 41:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, datafile, FILE_MODE_READ, &viewer));

 43:   PetscCall(ISCreate(PETSC_COMM_WORLD, &is_thermal));
 44:   PetscCall(ISLoad(is_thermal, viewer));
 45:   PetscCall(ISGetSize(is_thermal, &size));
 46:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "thermal field size = %" PetscInt_FMT " \n", size));

 48:   PetscCall(ISCreate(PETSC_COMM_WORLD, &is_mech));
 49:   PetscCall(ISLoad(is_mech, viewer));
 50:   PetscCall(ISGetSize(is_mech, &size));
 51:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mechanical field size = %" PetscInt_FMT " \n", size));
 52:   PetscCall(PetscViewerDestroy(&viewer));

 54:   PetscCall(MatCreateVecs(A, &x, &b));

 56:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 57:   PetscCall(KSPSetOperators(ksp, A, A));

 59:   PetscCall(KSPSetType(ksp, KSPPREONLY));

 61:   PetscCall(KSPGetPC(ksp, &pc));
 62:   PetscCall(PCSetType(pc, PCFIELDSPLIT));
 63:   PetscCall(PCFieldSplitSetIS(pc, "thermal", is_thermal));
 64:   PetscCall(PCFieldSplitSetIS(pc, "mechanical", is_mech));
 65:   PetscCall(PCSetUp(pc));
 66:   PetscCall(PCFieldSplitGetSubKSP(pc, &n_splits, &sub_ksp));
 67:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "n_splits = %" PetscInt_FMT " \n", n_splits));
 68:   PetscCall(KSPSetType(sub_ksp[0], KSPGMRES));

 70:   PetscCall(KSPGetPC(sub_ksp[0], &pc_thermal));
 71:   PetscCall(PCSetType(pc_thermal, PCJACOBI));
 72:   PetscCall(KSPSetFromOptions(sub_ksp[0]));

 74:   PetscCall(KSPSetType(sub_ksp[1], KSPCG));
 75:   PetscCall(KSPGetPC(sub_ksp[1], &pc_mech));
 76:   PetscCall(PCSetType(pc_mech, PCGAMG));
 77:   PetscCall(KSPSetFromOptions(sub_ksp[1]));

 79:   PetscCall(PetscStrcpy(datafile, datafilespath));
 80:   PetscCall(PetscStrcat(datafile, "/lostnullspace/"));
 81:   PetscCall(PetscStrcat(datafile, "rigid-modes.bin"));
 82:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, datafile, FILE_MODE_READ, &viewer));
 83:   for (PetscInt i = 0; i < 6; i++) {
 84:     PetscCall(VecCreate(PETSC_COMM_WORLD, &rigid_mode[i]));
 85:     PetscCall(VecLoad(rigid_mode[i], viewer));
 86:   }
 87:   PetscCall(PetscViewerDestroy(&viewer));

 89:   PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_FALSE, 6, rigid_mode, &near_null_space));
 90:   PetscCall(KSPGetOperators(sub_ksp[1], &K, PETSC_NULLPTR));
 91:   PetscCall(MatSetNearNullSpace(K, near_null_space));
 92:   PetscCall(MatSetBlockSize(K, 3));

 94:   PetscCall(MatGetSize(K, &rows, &cols));
 95:   PetscCall(MatGetBlockSize(K, &bs));
 96:   PetscCall(MatGetNearNullSpace(K, &nullsp));
 97:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "K has rows = %" PetscInt_FMT ", cols = %" PetscInt_FMT ", bs = %" PetscInt_FMT ", nearnullsp = %p\n", rows, cols, bs, (void *)nullsp));

 99:   PetscCall(ISGetIndices(is_thermal, &bc_thermal_indexes));
100:   PetscCall(ISGetIndices(is_mech, &bc_mech_indexes));

102:   // check that MatZeroRows() without MAT_KEEP_NONZERO_PATTERN does not remove the near null spaces attached to the submatrices
103:   PetscCall(PetscOptionsHasName(PETSC_NULLPTR, PETSC_NULLPTR, "-columns", &has_columns));
104:   if (has_columns == PETSC_TRUE) {
105:     PetscCall(MatZeroRowsColumns(A, 3, bc_mech_indexes, 1, PETSC_NULLPTR, PETSC_NULLPTR));
106:     PetscCall(MatZeroRowsColumns(A, 1, bc_thermal_indexes, 1, PETSC_NULLPTR, PETSC_NULLPTR));
107:   } else {
108:     PetscCall(MatZeroRows(A, 3, bc_mech_indexes, 1, PETSC_NULLPTR, PETSC_NULLPTR));
109:     PetscCall(MatZeroRows(A, 1, bc_thermal_indexes, 1, PETSC_NULLPTR, PETSC_NULLPTR));
110:   }
111:   PetscCall(ISRestoreIndices(is_mech, &bc_mech_indexes));
112:   PetscCall(ISRestoreIndices(is_thermal, &bc_thermal_indexes));

114:   PetscCall(KSPSetFromOptions(ksp));
115:   PetscCall(KSPSetUp(ksp));
116:   PetscCall(KSPSolve(ksp, b, x));

118:   PetscCall(MatDestroy(&A));
119:   PetscCall(KSPDestroy(&ksp));
120:   PetscCall(ISDestroy(&is_mech));
121:   PetscCall(ISDestroy(&is_thermal));
122:   PetscCall(VecDestroy(&x));
123:   PetscCall(VecDestroy(&b));
124:   PetscCall(MatNullSpaceDestroy(&near_null_space));
125:   for (PetscInt i = 0; i < 6; i++) PetscCall(VecDestroy(&rigid_mode[i]));
126:   PetscCall(PetscFree(sub_ksp));
127:   PetscCall(PetscFinalize());
128:   return 0;
129: }

131: /*TEST

133:     test:
134:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
135:       args: -datafilespath ${DATAFILESPATH} -ksp_view
136:       filter: grep "near null"

138: TEST*/