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*/