Actual source code: ex28.c


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

  4: #include <petscksp.h>

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

 18:   PetscInitialize(&argc, &args, (char *)0, help);
 19:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 20:   PetscOptionsGetBool(NULL, NULL, "-procedural", &TEST_PROCEDURAL, NULL);

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

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

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

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

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

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

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

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

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

111:     KSPGetPC(ksp, &pc);
112:     PCSetType(pc, PCREDUNDANT);
113:     MPI_Comm_size(PETSC_COMM_WORLD, &size);
114:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
116:     PCRedundantSetNumber(pc, size - 2);
117:     KSPSetFromOptions(ksp);

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

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

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

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

153: /*TEST

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

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

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

169: TEST*/