Actual source code: ex12.c

  1: static char help[] = "Solves a linear system in parallel with KSP,\n\
  2: demonstrating how to register a new preconditioner (PC) type.\n\
  3: Input parameters include:\n\
  4:   -m <mesh_x>       : number of mesh points in x-direction\n\
  5:   -n <mesh_y>       : number of mesh points in y-direction\n\n";

  7: /*
  8:    Demonstrates registering a new preconditioner (PC) type.

 10:    To register a PC type whose code is linked into the executable,
 11:    use PCRegister(). To register a PC type in a dynamic library use PCRegister()

 13:    Also provide the prototype for your PCCreate_XXX() function. In
 14:    this example we use the PETSc implementation of the Jacobi method,
 15:    PCCreate_Jacobi() just as an example.

 17:    See the file src/ksp/pc/impls/jacobi/jacobi.c for details on how to
 18:    write a new PC component.

 20:    See the manual page PCRegister() for details on how to register a method.
 21: */

 23: /*
 24:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 25:   automatically includes:
 26:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 27:      petscmat.h - matrices
 28:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 29:      petscviewer.h - viewers               petscpc.h  - preconditioners
 30: */
 31: #include <petscksp.h>

 33: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC);

 35: int main(int argc, char **args)
 36: {
 37:   Vec         x, b, u; /* approx solution, RHS, exact solution */
 38:   Mat         A;       /* linear system matrix */
 39:   KSP         ksp;     /* linear solver context */
 40:   PetscReal   norm;    /* norm of solution error */
 41:   PetscInt    i, j, Ii, J, Istart, Iend, m = 8, n = 7, its;
 42:   PetscScalar v, one = 1.0;
 43:   PC          pc; /* preconditioner context */

 45:   PetscFunctionBeginUser;
 46:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 47:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 48:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

 50:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 51:          Compute the matrix and right-hand-side vector that define
 52:          the linear system, Ax = b.
 53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 54:   /*
 55:      Create parallel matrix, specifying only its global dimensions.
 56:      When using MatCreate(), the matrix format can be specified at
 57:      runtime. Also, the parallel partitioning of the matrix can be
 58:      determined by PETSc at runtime.
 59:   */
 60:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 61:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
 62:   PetscCall(MatSetFromOptions(A));
 63:   PetscCall(MatSetUp(A));

 65:   /*
 66:      Currently, all PETSc parallel matrix formats are partitioned by
 67:      contiguous chunks of rows across the processors.  Determine which
 68:      rows of the matrix are locally owned.
 69:   */
 70:   PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));

 72:   /*
 73:      Set matrix elements for the 2-D, five-point stencil in parallel.
 74:       - Each processor needs to insert only elements that it owns
 75:         locally (but any non-local elements will be sent to the
 76:         appropriate processor during matrix assembly).
 77:       - Always specify global rows and columns of matrix entries.
 78:    */
 79:   for (Ii = Istart; Ii < Iend; Ii++) {
 80:     v = -1.0;
 81:     i = Ii / n;
 82:     j = Ii - i * n;
 83:     if (i > 0) {
 84:       J = Ii - n;
 85:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 86:     }
 87:     if (i < m - 1) {
 88:       J = Ii + n;
 89:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 90:     }
 91:     if (j > 0) {
 92:       J = Ii - 1;
 93:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 94:     }
 95:     if (j < n - 1) {
 96:       J = Ii + 1;
 97:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 98:     }
 99:     v = 4.0;
100:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
101:   }

103:   /*
104:      Assemble matrix, using the 2-step process:
105:        MatAssemblyBegin(), MatAssemblyEnd()
106:      Computations can be done while messages are in transition
107:      by placing code between these two statements.
108:   */
109:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
110:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

112:   /*
113:      Create parallel vectors.
114:       - When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
115:       we specify only the vector's global
116:         dimension; the parallel partitioning is determined at runtime.
117:       - When solving a linear system, the vectors and matrices MUST
118:         be partitioned accordingly.  PETSc automatically generates
119:         appropriately partitioned matrices and vectors when MatCreate()
120:         and VecCreate() are used with the same communicator.
121:       - Note: We form 1 vector from scratch and then duplicate as needed.
122:   */
123:   PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
124:   PetscCall(VecSetSizes(u, PETSC_DECIDE, m * n));
125:   PetscCall(VecSetFromOptions(u));
126:   PetscCall(VecDuplicate(u, &b));
127:   PetscCall(VecDuplicate(b, &x));

129:   /*
130:      Set exact solution; then compute right-hand-side vector.
131:      Use an exact solution of a vector with all elements of 1.0;
132:   */
133:   PetscCall(VecSet(u, one));
134:   PetscCall(MatMult(A, u, b));

136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137:                 Create the linear solver and set various options
138:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

140:   /*
141:      Create linear solver context
142:   */
143:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

145:   /*
146:      Set operators. Here the matrix that defines the linear system
147:      also serves as the preconditioning matrix.
148:   */
149:   PetscCall(KSPSetOperators(ksp, A, A));

151:   /*
152:        First register a new PC type with the command PCRegister()
153:   */
154:   PetscCall(PCRegister("ourjacobi", PCCreate_Jacobi));

156:   /*
157:      Set the PC type to be the new method
158:   */
159:   PetscCall(KSPGetPC(ksp, &pc));
160:   PetscCall(PCSetType(pc, "ourjacobi"));

162:   /*
163:     Set runtime options, e.g.,
164:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
165:     These options will override those specified above as long as
166:     KSPSetFromOptions() is called _after_ any other customization
167:     routines.
168:   */
169:   PetscCall(KSPSetFromOptions(ksp));

171:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172:                       Solve the linear system
173:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

175:   PetscCall(KSPSolve(ksp, b, x));

177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:                       Check solution and clean up
179:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

181:   /*
182:      Check the error
183:   */
184:   PetscCall(VecAXPY(x, -1.0, u));
185:   PetscCall(VecNorm(x, NORM_2, &norm));
186:   PetscCall(KSPGetIterationNumber(ksp, &its));

188:   /*
189:      Print convergence information.  PetscPrintf() produces a single
190:      print statement from all processes that share a communicator.
191:   */
192:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));

194:   /*
195:      Free work space.  All PETSc objects should be destroyed when they
196:      are no longer needed.
197:   */
198:   PetscCall(KSPDestroy(&ksp));
199:   PetscCall(VecDestroy(&u));
200:   PetscCall(VecDestroy(&x));
201:   PetscCall(VecDestroy(&b));
202:   PetscCall(MatDestroy(&A));

204:   /*
205:      Always call PetscFinalize() before exiting a program.  This routine
206:        - finalizes the PETSc libraries as well as MPI
207:        - provides summary and diagnostic information if certain runtime
208:          options are chosen (e.g., -log_view).
209:   */
210:   PetscCall(PetscFinalize());
211:   return 0;
212: }

214: /*TEST

216:    test:
217:       args: -ksp_gmres_cgs_refinement_type refine_always

219: TEST*/