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, M = 1;
32: PetscBool flg;
34: PetscFunctionBeginUser;
35: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
37: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
38: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
39: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nBlock-structured matrix, n=%" PetscInt_FMT "\n\n", n));
41: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: Compute the blocks R and C, and the Nest H
43: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
45: #if defined(PETSC_USE_COMPLEX)
46: a = PetscCMPLX(-0.1, 0.2);
47: b = PetscCMPLX(1.0, 0.5);
48: d = PetscCMPLX(2.0, 0.2);
49: #else
50: a = -0.1;
51: b = 1.0;
52: d = 2.0;
53: #endif
54: c = 4.5;
56: PetscCall(MatCreate(PETSC_COMM_WORLD, &R));
57: PetscCall(MatSetSizes(R, PETSC_DECIDE, PETSC_DECIDE, n, n));
58: PetscCall(MatSetFromOptions(R));
60: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
61: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, n, n));
62: PetscCall(MatSetFromOptions(C));
64: PetscCall(MatGetOwnershipRange(R, &Istart, &Iend));
65: for (i = Istart; i < Iend; i++) {
66: if (i > 1) PetscCall(MatSetValue(R, i, i - 2, a, INSERT_VALUES));
67: if (i > 0) PetscCall(MatSetValue(R, i, i - 1, b, INSERT_VALUES));
68: PetscCall(MatSetValue(R, i, i, c, INSERT_VALUES));
69: if (i < n - 1) PetscCall(MatSetValue(R, i, i + 1, PetscConj(b), INSERT_VALUES));
70: if (i < n - 2) PetscCall(MatSetValue(R, i, i + 2, PetscConj(a), INSERT_VALUES));
71: }
73: PetscCall(MatGetOwnershipRange(C, &Istart, &Iend));
74: for (i = Istart; i < Iend; i++) {
75: if (i > 0) PetscCall(MatSetValue(C, i, i - 1, b, INSERT_VALUES));
76: PetscCall(MatSetValue(C, i, i, d, INSERT_VALUES));
77: if (i < n - 1) PetscCall(MatSetValue(C, i, i + 1, b, INSERT_VALUES));
78: }
80: PetscCall(MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY));
81: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
82: PetscCall(MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY));
83: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
85: /* top block row */
86: block[0] = R;
87: block[1] = C;
89: /* bottom block row */
90: PetscCall(MatCreateHermitianTranspose(C, &block[2]));
91: PetscCall(MatScale(block[2], -1.0));
92: PetscCall(MatCreateTranspose(R, &block[3]));
93: PetscCall(MatScale(block[3], -1.0));
95: /* create nest matrix */
96: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)R), 2, NULL, 2, NULL, block, &H));
97: PetscCall(MatDestroy(&block[2]));
98: PetscCall(MatDestroy(&block[3]));
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create linear system and solve
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
105: PetscCall(KSPSetOperators(ksp, H, H));
106: PetscCall(KSPSetTolerances(ksp, tol, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
107: PetscCall(KSPSetFromOptions(ksp));
109: PetscCall(MatCreateVecs(H, &x, &rhs));
110: if (M == 1) {
111: PetscCall(VecSet(rhs, 1.0));
112: PetscCall(KSPSolve(ksp, rhs, x));
114: /* check final residual */
115: PetscCall(VecDuplicate(rhs, &r));
116: PetscCall(MatMult(H, x, r));
117: PetscCall(VecAXPY(r, -1.0, rhs));
118: PetscCall(VecNorm(r, NORM_2, norm));
119: PetscCall(VecNorm(rhs, NORM_2, norm + 1));
120: 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]));
121: PetscCall(VecDestroy(&r));
122: } else {
123: Mat B, X, R;
124: VecType type;
126: PetscCall(MatGetLocalSize(H, &n, NULL));
127: PetscCall(VecGetType(rhs, &type));
128: PetscCall(MatCreateDenseFromVecType(PETSC_COMM_WORLD, type, n, PETSC_DECIDE, PETSC_DECIDE, M, -1, NULL, &B));
129: PetscCall(MatCreateDenseFromVecType(PETSC_COMM_WORLD, type, n, PETSC_DECIDE, PETSC_DECIDE, M, -1, NULL, &X));
130: PetscCall(MatSetRandom(B, NULL));
131: PetscCall(KSPMatSolve(ksp, B, X));
132: /* check final residual */
133: PetscCall(MatMatMult(H, X, MAT_INITIAL_MATRIX, PETSC_DECIDE, &R));
134: PetscCall(MatAXPY(R, -1.0, B, SAME_NONZERO_PATTERN));
135: PetscCall(MatNorm(R, NORM_1, norm));
136: PetscCall(MatNorm(B, NORM_1, norm + 1));
137: PetscCheck(norm[0] / norm[1] < 10.0 * PETSC_SMALL, PetscObjectComm((PetscObject)H), PETSC_ERR_PLIB, "Error ||H X-B||_1 / ||B||_1: %1.6e", (double)(norm[0] / norm[1]));
138: PetscCall(MatDestroy(&R));
139: PetscCall(MatDestroy(&X));
140: PetscCall(MatDestroy(&B));
141: }
143: /* check PetscMemType */
144: PetscCall(KSPGetPC(ksp, &pc));
145: PetscCall(PCGetType(pc, &type));
146: PetscCall(PetscStrcmp(type, PCFIELDSPLIT, &flg));
147: if (flg) {
148: KSP *subksp;
149: Mat pmat;
150: const PetscScalar *array;
151: PetscInt n;
152: PetscMemType type[2];
154: PetscCall(PCFieldSplitGetSubKSP(pc, &n, &subksp));
155: PetscCheck(n == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Wrong number of fields");
156: PetscCall(KSPGetOperators(subksp[1], NULL, &pmat));
157: PetscCall(PetscFree(subksp));
158: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)pmat, &flg, MATSEQDENSE, MATMPIDENSE, ""));
159: if (flg) {
160: PetscCall(VecGetArrayReadAndMemType(x, &array, type));
161: PetscCall(VecRestoreArrayReadAndMemType(x, &array));
162: PetscCall(MatDenseGetArrayReadAndMemType(pmat, &array, type + 1));
163: PetscCall(MatDenseRestoreArrayReadAndMemType(pmat, &array));
164: PetscCheck(type[0] == type[1], PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Failed PetscMemType comparison");
165: }
166: }
168: PetscCall(KSPDestroy(&ksp));
169: PetscCall(MatDestroy(&R));
170: PetscCall(MatDestroy(&C));
171: PetscCall(MatDestroy(&H));
172: PetscCall(VecDestroy(&rhs));
173: PetscCall(VecDestroy(&x));
174: PetscCall(PetscFinalize());
175: return 0;
176: }
178: /*TEST
180: testset:
181: requires: !complex
182: output_file: output/ex87.out
183: test:
184: suffix: real
185: args: -ksp_pc_side right
186: test:
187: suffix: real_fieldsplit
188: 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 -M {{1 2}shared output}
189: test:
190: requires: cuda
191: nsize: {{1 4}}
192: suffix: real_fieldsplit_cuda
193: 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 -M {{1 2}shared output}
194: test:
195: requires: hip
196: nsize: 4 # this is broken with a single process, see https://gitlab.com/petsc/petsc/-/issues/1529
197: suffix: real_fieldsplit_hip
198: 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 -M {{1 2}shared output}
199: test:
200: requires: hpddm
201: nsize: 4
202: suffix: real_fieldsplit_hpddm
203: args: -ksp_type hpddm -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 {{diag lower upper}shared output} -pc_fieldsplit_schur_precondition full -M 2
205: testset:
206: requires: complex
207: output_file: output/ex87.out
208: test:
209: suffix: complex
210: args: -ksp_pc_side right
211: test:
212: requires: elemental
213: nsize: 4
214: suffix: complex_fieldsplit_elemental
215: 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
216: test:
217: requires: scalapack
218: nsize: 4
219: suffix: complex_fieldsplit_scalapack
220: 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
221: test:
222: suffix: complex_fieldsplit
223: 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 -M {{1 2}shared output}
224: test:
225: requires: hpddm
226: suffix: complex_fieldsplit_hpddm
227: args: -ksp_type hpddm -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 {{diag lower upper}shared output} -pc_fieldsplit_schur_precondition full -fieldsplit_1_explicit_operator_mat_hermitian -M 2
229: TEST*/