Actual source code: ex4.c
1: static char help[] = "Tests BDDC reuse after PCReset() and KSPReset().\n\n";
3: #include <petsc/private/pcbddcimpl.h>
5: int main(int argc, char **args)
6: {
7: Mat A, local;
8: MatNullSpace nsp;
9: ISLocalToGlobalMapping map;
10: KSP ksp;
11: PC pc;
12: PC_BDDC *bddc;
13: Vec modes[2], x, b, exact, scaling;
14: PetscInt *indices, *xadj, *adjncy;
15: PetscInt n, nlocal, start, end, nconstraints;
16: PetscMPIInt rank, size, active;
17: PetscBool empty_rank = PETSC_FALSE, user_graph = PETSC_FALSE, use_nnsp, change, deluxe, stiffness;
18: PetscScalar factor, weight;
19: PetscScalar *values;
20: PetscReal error;
22: PetscFunctionBeginUser;
23: PetscCall(PetscInitialize(&argc, &args, NULL, help));
24: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
25: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
26: PetscCall(PetscOptionsGetBool(NULL, NULL, "-empty_rank", &empty_rank, NULL));
27: PetscCall(PetscOptionsGetBool(NULL, NULL, "-user_graph", &user_graph, NULL));
28: active = size - (empty_rank ? 1 : 0);
29: PetscCheck(active > 0, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "At least one nonempty subdomain is required");
30: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
31: PetscCall(KSPSetType(ksp, KSPCG));
32: PetscCall(KSPGetPC(ksp, &pc));
33: PetscCall(PCSetType(pc, PCBDDC));
34: PetscCall(PCBDDCSetCoarseningRatio(pc, 1));
35: PetscCall(KSPSetFromOptions(ksp));
36: PetscCall(PCISSetSubdomainScalingFactor(pc, rank + 1.0));
37: PetscCall(PCISSetUseStiffnessScaling(pc, PETSC_TRUE));
38: bddc = (PC_BDDC *)pc->data;
39: for (PetscInt step = 0; step < 3; step++) {
40: use_nnsp = bddc->use_nnsp;
41: change = bddc->use_change_of_basis;
42: deluxe = bddc->use_deluxe_scaling;
43: stiffness = bddc->pcis.use_stiffness_scaling;
44: factor = bddc->pcis.scaling_factor;
45: // Reset before the first setup and twice between systems of different sizes.
46: PetscCall(PCReset(pc));
47: PetscCall(KSPReset(ksp));
48: PetscCheck(bddc->use_nnsp == use_nnsp && bddc->use_change_of_basis == change && bddc->use_deluxe_scaling == deluxe && bddc->coarsening_ratio == 1, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "BDDC options changed during reset");
49: PetscCheck(bddc->pcis.use_stiffness_scaling == stiffness && bddc->pcis.scaling_factor == factor, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Scaling options changed during reset");
51: n = 6 + 2 * step;
52: nlocal = rank < active ? n : 0;
53: PetscCall(PetscMalloc1(nlocal, &indices));
54: for (PetscInt i = 0; i < nlocal; i++) indices[i] = i;
55: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, nlocal, indices, PETSC_OWN_POINTER, &map));
56: PetscCall(MatCreateIS(PETSC_COMM_WORLD, 1, PETSC_DECIDE, PETSC_DECIDE, n, n, map, map, &A));
57: PetscCall(ISLocalToGlobalMappingDestroy(&map));
58: PetscCall(MatISSetPreallocation(A, n, NULL, n, NULL));
59: PetscCall(MatISGetLocalMat(A, &local));
60: for (PetscInt i = 0; i < nlocal; i++)
61: for (PetscInt j = 0; j < nlocal; j++) PetscCall(MatSetValue(local, i, j, (rank + 1.0) * (i == j ? 2.0 : 0.1), INSERT_VALUES));
62: PetscCall(MatISRestoreLocalMat(A, &local));
63: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
64: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
65: PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE));
66: if (user_graph) {
67: PetscCall(PetscMalloc1(nlocal + 1, &xadj));
68: PetscCall(PetscMalloc1(nlocal * nlocal, &adjncy));
69: for (PetscInt i = 0; i <= nlocal; i++) xadj[i] = i * nlocal;
70: for (PetscInt i = 0; i < nlocal * nlocal; i++) adjncy[i] = i % nlocal;
71: PetscCall(PCBDDCSetLocalAdjacencyGraph(pc, nlocal, xadj, adjncy, PETSC_OWN_POINTER));
72: }
73: PetscCall(PCISSetUseStiffnessScaling(pc, (PetscBool)(step == 0)));
74: PetscCall(PCISSetSubdomainScalingFactor(pc, rank + step + 1.0));
75: if (step == 2) {
76: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nlocal, &scaling));
77: PetscCall(VecSet(scaling, rank + 5.0));
78: PetscCall(PCISSetSubdomainDiagonalScaling(pc, scaling));
79: PetscCall(VecDestroy(&scaling));
80: }
82: PetscCall(MatCreateVecs(A, &modes[0], &b));
83: PetscCall(VecDuplicate(modes[0], &modes[1]));
84: PetscCall(VecDuplicate(modes[0], &x));
85: PetscCall(VecDuplicate(modes[0], &exact));
86: PetscCall(VecGetOwnershipRange(modes[0], &start, &end));
87: for (PetscInt k = 0; k < 2; k++) {
88: PetscCall(VecGetArray(modes[k], &values));
89: for (PetscInt i = start; i < end; i++) values[i - start] = (i == 2 * k ? 1.0 : 0.0) - (i == 2 * k + 1 ? 1.0 : 0.0);
90: PetscCall(VecRestoreArray(modes[k], &values));
91: PetscCall(VecNormalize(modes[k], NULL));
92: }
93: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_FALSE, 2, modes, &nsp));
94: PetscCall(MatSetNearNullSpace(A, nsp));
95: PetscCall(MatNullSpaceDestroy(&nsp));
96: PetscCall(VecSet(exact, 1.0));
97: PetscCall(VecAXPY(exact, 0.5, modes[1]));
98: PetscCall(MatMult(A, exact, b));
99: PetscCall(KSPSetOperators(ksp, A, A));
100: PetscCall(KSPSolve(ksp, b, x));
101: PetscCall(VecAXPY(x, -1.0, exact));
102: PetscCall(VecNorm(x, NORM_INFINITY, &error));
103: PetscCheck(error < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Solution error after reset: %g", (double)error);
104: PetscCall(MatGetSize(bddc->ConstraintMatrix, &nconstraints, NULL));
105: PetscCheck(nconstraints == (nlocal && active > 1 ? (use_nnsp ? 2 : 1) : 0), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected number of constraints after reset: %" PetscInt_FMT, nconstraints);
106: factor = step == 2 ? 5.0 : step + 1.0;
107: weight = (rank + factor) / (active * (active + 2.0 * factor - 1.0) / 2.0);
108: PetscCall(VecDuplicate(bddc->pcis.D, &scaling));
109: PetscCall(VecCopy(bddc->pcis.D, scaling));
110: PetscCall(VecShift(scaling, -weight));
111: PetscCall(VecNorm(scaling, NORM_INFINITY, &error));
112: PetscCheck(error < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Scaling error after reset: %g", (double)error);
113: PetscCall(VecDestroy(&scaling));
114: PetscCall(VecDestroy(&modes[0]));
115: PetscCall(VecDestroy(&modes[1]));
116: PetscCall(VecDestroy(&x));
117: PetscCall(VecDestroy(&b));
118: PetscCall(VecDestroy(&exact));
119: PetscCall(MatDestroy(&A));
120: }
121: // Changing type must remove the PCIS callbacks that access the BDDC context.
122: PetscCall(PCSetType(pc, PCNONE));
123: PetscCall(PCISSetSubdomainScalingFactor(pc, 2.0));
124: PetscCall(PCISSetUseStiffnessScaling(pc, PETSC_TRUE));
125: PetscCall(KSPDestroy(&ksp));
126: PetscCall(PetscFinalize());
127: return 0;
128: }
130: /*TEST
132: testset:
133: requires: double
134: output_file: output/empty.out
135: args: -ksp_error_if_not_converged -ksp_rtol 1e-12
136: test:
137: suffix: reset
138: nsize: {{1 2 3}}
139: args: -pc_bddc_use_change_of_basis {{0 1}} -user_graph {{0 1}}
140: test:
141: suffix: no_nnsp
142: nsize: 2
143: args: -pc_bddc_use_nnsp 0
144: test:
145: suffix: deluxe
146: nsize: 2
147: args: -pc_bddc_use_deluxe_scaling -pc_bddc_use_change_of_basis {{0 1}}
148: test:
149: suffix: empty
150: nsize: 3
151: args: -empty_rank -user_graph
153: TEST*/