Actual source code: ex28.c

  1: static char help[] = "Test procedural KSPSetFromOptions() or at runtime; Test PCREDUNDANT.\n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   Vec         x, b, u; /* approx solution, RHS, exact solution */
  8:   Mat         A;       /* linear system matrix */
  9:   KSP         ksp;     /* linear solver context */
 10:   PC          pc;      /* preconditioner context */
 11:   PetscReal   norm;    /* norm of solution error */
 12:   PetscInt    i, n = 10, col[3], its, rstart, rend, nlocal;
 13:   PetscScalar one             = 1.0, value[3];
 14:   PetscBool   TEST_PROCEDURAL = PETSC_FALSE;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 18:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 19:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-procedural", &TEST_PROCEDURAL, NULL));

 21:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 22:          Compute the matrix and right-hand-side vector that define
 23:          the linear system, Ax = b.
 24:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 26:   /*
 27:      Create vectors.  Note that we form 1 vector from scratch and
 28:      then duplicate as needed. For this simple case let PETSc decide how
 29:      many elements of the vector are stored on each processor. The second
 30:      argument to VecSetSizes() below causes PETSc to decide.
 31:   */
 32:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 33:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
 34:   PetscCall(VecSetFromOptions(x));
 35:   PetscCall(VecDuplicate(x, &b));
 36:   PetscCall(VecDuplicate(x, &u));

 38:   /* Identify the starting and ending mesh points on each
 39:      processor for the interior part of the mesh. We let PETSc decide
 40:      above. */

 42:   PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
 43:   PetscCall(VecGetLocalSize(x, &nlocal));

 45:   /* Create a tridiagonal matrix. See ../tutorials/ex23.c */
 46:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 47:   PetscCall(MatSetSizes(A, nlocal, nlocal, n, n));
 48:   PetscCall(MatSetFromOptions(A));
 49:   PetscCall(MatSetUp(A));
 50:   /* Assemble matrix */
 51:   if (!rstart) {
 52:     rstart   = 1;
 53:     i        = 0;
 54:     col[0]   = 0;
 55:     col[1]   = 1;
 56:     value[0] = 2.0;
 57:     value[1] = -1.0;
 58:     PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
 59:   }
 60:   if (rend == n) {
 61:     rend     = n - 1;
 62:     i        = n - 1;
 63:     col[0]   = n - 2;
 64:     col[1]   = n - 1;
 65:     value[0] = -1.0;
 66:     value[1] = 2.0;
 67:     PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
 68:   }

 70:   /* Set entries corresponding to the mesh interior */
 71:   value[0] = -1.0;
 72:   value[1] = 2.0;
 73:   value[2] = -1.0;
 74:   for (i = rstart; i < rend; i++) {
 75:     col[0] = i - 1;
 76:     col[1] = i;
 77:     col[2] = i + 1;
 78:     PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 79:   }
 80:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 81:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 83:   /* Set exact solution; then compute right-hand-side vector. */
 84:   PetscCall(VecSet(u, one));
 85:   PetscCall(MatMult(A, u, b));

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88:                 Create the linear solver and set various options
 89:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 90:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 91:   PetscCall(KSPSetOperators(ksp, A, A));

 93:   /*
 94:      Set linear solver defaults for this problem (optional).
 95:      - By extracting the KSP and PC contexts from the KSP context,
 96:        we can then directly call any KSP and PC routines to set
 97:        various options.
 98:      - The following statements are optional; all of these
 99:        parameters could alternatively be specified at runtime via
100:        KSPSetFromOptions();
101:   */
102:   if (TEST_PROCEDURAL) {
103:     /* Example of runtime options: '-pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi' */
104:     PetscMPIInt size, rank, subsize;
105:     Mat         A_redundant;
106:     KSP         innerksp;
107:     PC          innerpc;
108:     MPI_Comm    subcomm;

110:     PetscCall(KSPGetPC(ksp, &pc));
111:     PetscCall(PCSetType(pc, PCREDUNDANT));
112:     PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
113:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
114:     PetscCheck(size > 2, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "Num of processes %d must greater than 2", size);
115:     PetscCall(PCRedundantSetNumber(pc, size - 2));
116:     PetscCall(KSPSetFromOptions(ksp));

118:     /* Get subcommunicator and redundant matrix */
119:     PetscCall(KSPSetUp(ksp));
120:     PetscCall(PCRedundantGetKSP(pc, &innerksp));
121:     PetscCall(KSPGetPC(innerksp, &innerpc));
122:     PetscCall(PCGetOperators(innerpc, NULL, &A_redundant));
123:     PetscCall(PetscObjectGetComm((PetscObject)A_redundant, &subcomm));
124:     PetscCallMPI(MPI_Comm_size(subcomm, &subsize));
125:     if (subsize == 1 && rank == 0) {
126:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "A_redundant:\n"));
127:       PetscCall(MatView(A_redundant, PETSC_VIEWER_STDOUT_SELF));
128:     }
129:   } else {
130:     PetscCall(KSPSetFromOptions(ksp));
131:   }

133:   /*  Solve linear system */
134:   PetscCall(KSPSolve(ksp, b, x));

136:   /* Check the error */
137:   PetscCall(VecAXPY(x, -1.0, u));
138:   PetscCall(VecNorm(x, NORM_2, &norm));
139:   PetscCall(KSPGetIterationNumber(ksp, &its));
140:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));

142:   /* Free work space. */
143:   PetscCall(VecDestroy(&x));
144:   PetscCall(VecDestroy(&u));
145:   PetscCall(VecDestroy(&b));
146:   PetscCall(MatDestroy(&A));
147:   PetscCall(KSPDestroy(&ksp));
148:   PetscCall(PetscFinalize());
149:   return 0;
150: }

152: /*TEST

154:     test:
155:       nsize: 3
156:       output_file: output/ex28.out

158:     test:
159:       suffix: 2
160:       args: -procedural -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi
161:       nsize: 3

163:     test:
164:       suffix: 3
165:       args: -procedural -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi
166:       nsize: 5

168: TEST*/