Actual source code: ex87.c
1: static char help[] = "Block-structured Nest matrix involving a HermitianTranspose block.\n\n"
2: "The command line options are:\n"
3: " -n <n>, where <n> = dimension of the blocks.\n\n";
5: #include <petscksp.h>
7: /*
8: Solves a linear system with coefficient matrix
10: H = [ R C
11: -C^H -R^T ],
13: where R is Hermitian and C is complex symmetric. In particular, R and C have the
14: following Toeplitz structure:
16: R = pentadiag{a,b,c,conj(b),conj(a)}
17: C = tridiag{b,d,b}
19: where a,b,d are complex scalars, and c is real.
20: */
22: int main(int argc, char **argv)
23: {
24: Mat H, R, C, block[4];
25: Vec rhs, x, r;
26: KSP ksp;
27: PC pc;
28: PCType type;
29: PetscReal norm[2], tol = 100.0 * PETSC_MACHINE_EPSILON;
30: PetscScalar a, b, c, d;
31: PetscInt n = 24, Istart, Iend, i;
32: PetscBool flg;
34: PetscFunctionBeginUser;
35: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
37: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
38: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nBlock-structured matrix, n=%" PetscInt_FMT "\n\n", n));
40: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: Compute the blocks R and C, and the Nest H
42: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
44: #if defined(PETSC_USE_COMPLEX)
45: a = PetscCMPLX(-0.1, 0.2);
46: b = PetscCMPLX(1.0, 0.5);
47: d = PetscCMPLX(2.0, 0.2);
48: #else
49: a = -0.1;
50: b = 1.0;
51: d = 2.0;
52: #endif
53: c = 4.5;
55: PetscCall(MatCreate(PETSC_COMM_WORLD, &R));
56: PetscCall(MatSetSizes(R, PETSC_DECIDE, PETSC_DECIDE, n, n));
57: PetscCall(MatSetFromOptions(R));
59: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
60: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, n, n));
61: PetscCall(MatSetFromOptions(C));
63: PetscCall(MatGetOwnershipRange(R, &Istart, &Iend));
64: for (i = Istart; i < Iend; i++) {
65: if (i > 1) PetscCall(MatSetValue(R, i, i - 2, a, INSERT_VALUES));
66: if (i > 0) PetscCall(MatSetValue(R, i, i - 1, b, INSERT_VALUES));
67: PetscCall(MatSetValue(R, i, i, c, INSERT_VALUES));
68: if (i < n - 1) PetscCall(MatSetValue(R, i, i + 1, PetscConj(b), INSERT_VALUES));
69: if (i < n - 2) PetscCall(MatSetValue(R, i, i + 2, PetscConj(a), INSERT_VALUES));
70: }
72: PetscCall(MatGetOwnershipRange(C, &Istart, &Iend));
73: for (i = Istart; i < Iend; i++) {
74: if (i > 0) PetscCall(MatSetValue(C, i, i - 1, b, INSERT_VALUES));
75: PetscCall(MatSetValue(C, i, i, d, INSERT_VALUES));
76: if (i < n - 1) PetscCall(MatSetValue(C, i, i + 1, b, INSERT_VALUES));
77: }
79: PetscCall(MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY));
80: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
81: PetscCall(MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY));
82: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
84: /* top block row */
85: block[0] = R;
86: block[1] = C;
88: /* bottom block row */
89: PetscCall(MatCreateHermitianTranspose(C, &block[2]));
90: PetscCall(MatScale(block[2], -1.0));
91: PetscCall(MatCreateTranspose(R, &block[3]));
92: PetscCall(MatScale(block[3], -1.0));
94: /* create nest matrix */
95: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)R), 2, NULL, 2, NULL, block, &H));
96: PetscCall(MatDestroy(&block[2]));
97: PetscCall(MatDestroy(&block[3]));
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Create linear system and solve
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
104: PetscCall(KSPSetOperators(ksp, H, H));
105: PetscCall(KSPSetTolerances(ksp, tol, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
106: PetscCall(KSPSetFromOptions(ksp));
108: PetscCall(MatCreateVecs(H, &x, &rhs));
109: PetscCall(VecSet(rhs, 1.0));
110: PetscCall(KSPSolve(ksp, rhs, x));
112: /* check final residual */
113: PetscCall(VecDuplicate(rhs, &r));
114: PetscCall(MatMult(H, x, r));
115: PetscCall(VecAXPY(r, -1.0, rhs));
116: PetscCall(VecNorm(r, NORM_2, norm));
117: PetscCall(VecNorm(rhs, NORM_2, norm + 1));
118: PetscCheck(norm[0] / norm[1] < 10.0 * PETSC_SMALL, PetscObjectComm((PetscObject)H), PETSC_ERR_PLIB, "Error ||H x-rhs||_2 / ||rhs||_2: %1.6e", (double)(norm[0] / norm[1]));
120: /* check PetscMemType */
121: PetscCall(KSPGetPC(ksp, &pc));
122: PetscCall(PCGetType(pc, &type));
123: PetscCall(PetscStrcmp(type, PCFIELDSPLIT, &flg));
124: if (flg) {
125: KSP *subksp;
126: Mat pmat;
127: const PetscScalar *array;
128: PetscInt n;
129: PetscMemType type[2];
131: PetscCall(PCFieldSplitGetSubKSP(pc, &n, &subksp));
132: PetscCheck(n == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Wrong number of fields");
133: PetscCall(KSPGetOperators(subksp[1], NULL, &pmat));
134: PetscCall(PetscFree(subksp));
135: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)pmat, &flg, MATSEQDENSE, MATMPIDENSE, ""));
136: if (flg) {
137: PetscCall(VecGetArrayReadAndMemType(x, &array, type));
138: PetscCall(VecRestoreArrayReadAndMemType(x, &array));
139: PetscCall(MatDenseGetArrayReadAndMemType(pmat, &array, type + 1));
140: PetscCall(MatDenseRestoreArrayReadAndMemType(pmat, &array));
141: PetscCheck(type[0] == type[1], PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Failed PetscMemType comparison");
142: }
143: }
145: PetscCall(KSPDestroy(&ksp));
146: PetscCall(MatDestroy(&R));
147: PetscCall(MatDestroy(&C));
148: PetscCall(MatDestroy(&H));
149: PetscCall(VecDestroy(&rhs));
150: PetscCall(VecDestroy(&x));
151: PetscCall(VecDestroy(&r));
152: PetscCall(PetscFinalize());
153: return 0;
154: }
156: /*TEST
158: testset:
159: requires: !complex
160: output_file: output/ex87.out
161: test:
162: suffix: real
163: args: -ksp_pc_side right
164: test:
165: suffix: real_fieldsplit
166: args: -ksp_type preonly -pc_type fieldsplit -fieldsplit_ksp_type preonly -fieldsplit_0_pc_type lu -fieldsplit_1_pc_type cholesky -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full
167: test:
168: requires: cuda
169: nsize: {{1 4}}
170: suffix: real_fieldsplit_cuda
171: args: -ksp_type preonly -pc_type fieldsplit -fieldsplit_ksp_type preonly -fieldsplit_pc_type redundant -fieldsplit_redundant_ksp_type preonly -fieldsplit_redundant_pc_type lu -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full -mat_type aijcusparse
172: test:
173: requires: hip
174: nsize: 4 # this is broken with a single process, see https://gitlab.com/petsc/petsc/-/issues/1529
175: suffix: real_fieldsplit_hip
176: args: -ksp_type preonly -pc_type fieldsplit -fieldsplit_ksp_type preonly -fieldsplit_pc_type redundant -fieldsplit_redundant_ksp_type preonly -fieldsplit_redundant_pc_type lu -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full -mat_type aijhipsparse
178: testset:
179: requires: complex
180: output_file: output/ex87.out
181: test:
182: suffix: complex
183: args: -ksp_pc_side right
184: test:
185: requires: elemental
186: nsize: 4
187: suffix: complex_fieldsplit_elemental
188: args: -ksp_type preonly -pc_type fieldsplit -fieldsplit_ksp_type preonly -fieldsplit_0_pc_type redundant -fieldsplit_0_redundant_pc_type lu -fieldsplit_1_pc_type lu -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full -fieldsplit_1_explicit_operator_mat_type elemental
189: test:
190: requires: scalapack
191: nsize: 4
192: suffix: complex_fieldsplit_scalapack
193: args: -ksp_type preonly -pc_type fieldsplit -fieldsplit_ksp_type preonly -fieldsplit_pc_type lu -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full -mat_type scalapack -fieldsplit_1_explicit_operator_mat_type scalapack
194: test:
195: suffix: complex_fieldsplit
196: args: -ksp_type preonly -pc_type fieldsplit -fieldsplit_ksp_type preonly -fieldsplit_0_pc_type lu -fieldsplit_1_pc_type cholesky -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full -fieldsplit_1_explicit_operator_mat_hermitian
198: TEST*/