Actual source code: ex40.c


  2: static char help[] = "Solves a linear system in parallel with KSP.\n\
  3: Input parameters include:\n\
  4:   -random_exact_sol : use a random exact solution vector\n\
  5:   -view_exact_sol   : write exact solution vector to stdout\n\
  6:   -m <mesh_x>       : number of mesh points in x-direction\n\
  7:   -n <mesh_y>       : number of mesh points in y-direction\n\n";

  9: /*
 10:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 11:   automatically includes:
 12:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 13:      petscmat.h - matrices
 14:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 15:      petscviewer.h - viewers               petscpc.h  - preconditioners
 16: */
 17: #include <petscksp.h>

 19: int main(int argc, char **args)
 20: {
 21:   Vec         x, b, u; /* approx solution, RHS, exact solution */
 22:   Mat         A;       /* linear system matrix */
 23:   KSP         ksp;     /* linear solver context */
 24:   PetscRandom rctx;    /* random number generator context */
 25:   PetscReal   norm;    /* norm of solution error */
 26:   PetscInt    i, j, Ii, J, m = 8, n = 7, its;
 27:   PetscBool   flg = PETSC_FALSE;
 28:   PetscScalar v;
 29:   PetscMPIInt rank;

 32:   PetscInitialize(&argc, &args, (char *)0, help);
 33:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 34:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 35:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 36:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 37:          Compute the matrix and right-hand-side vector that define
 38:          the linear system, Ax = b.
 39:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 40:   /*
 41:      Create parallel matrix, specifying only its global dimensions.
 42:      When using MatCreate(), the matrix format can be specified at
 43:      runtime. Also, the parallel partitioning of the matrix is
 44:      determined by PETSc at runtime.

 46:      Performance tuning note:  For problems of substantial size,
 47:      preallocation of matrix memory is crucial for attaining good
 48:      performance. See the matrix chapter of the users manual for details.
 49:   */
 50:   MatCreate(PETSC_COMM_WORLD, &A);
 51:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n);
 52:   MatSetType(A, MATELEMENTAL);
 53:   MatSetFromOptions(A);
 54:   MatSetUp(A);
 55:   if (rank == 0) {
 56:     PetscInt M, N;
 57:     MatGetSize(A, &M, &N);
 58:     for (Ii = 0; Ii < M; Ii++) {
 59:       v = -1.0;
 60:       i = Ii / n;
 61:       j = Ii - i * n;
 62:       if (i > 0) {
 63:         J = Ii - n;
 64:         MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 65:       }
 66:       if (i < m - 1) {
 67:         J = Ii + n;
 68:         MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 69:       }
 70:       if (j > 0) {
 71:         J = Ii - 1;
 72:         MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 73:       }
 74:       if (j < n - 1) {
 75:         J = Ii + 1;
 76:         MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 77:       }
 78:       v = 4.0;
 79:       MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES);
 80:     }
 81:   }
 82:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 83:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 85:   /* MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE); */

 87:   /*
 88:      Create parallel vectors.
 89:       - We form 1 vector from scratch and then duplicate as needed.
 90:       - When using VecCreate(), VecSetSizes and VecSetFromOptions()
 91:         in this example, we specify only the
 92:         vector's global dimension; the parallel partitioning is determined
 93:         at runtime.
 94:       - When solving a linear system, the vectors and matrices MUST
 95:         be partitioned accordingly.  PETSc automatically generates
 96:         appropriately partitioned matrices and vectors when MatCreate()
 97:         and VecCreate() are used with the same communicator.
 98:       - The user can alternatively specify the local vector and matrix
 99:         dimensions when more sophisticated partitioning is needed
100:         (replacing the PETSC_DECIDE argument in the VecSetSizes() statement
101:         below).
102:   */
103:   VecCreate(PETSC_COMM_WORLD, &u);
104:   VecSetSizes(u, PETSC_DECIDE, m * n);
105:   VecSetFromOptions(u);
106:   VecDuplicate(u, &b);
107:   VecDuplicate(b, &x);

109:   /*
110:      Set exact solution; then compute right-hand-side vector.
111:      By default we use an exact solution of a vector with all
112:      elements of 1.0;  Alternatively, using the runtime option
113:      -random_sol forms a solution vector with random components.
114:   */
115:   PetscOptionsGetBool(NULL, NULL, "-random_exact_sol", &flg, NULL);
116:   if (flg) {
117:     PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
118:     PetscRandomSetFromOptions(rctx);
119:     VecSetRandom(u, rctx);
120:     PetscRandomDestroy(&rctx);
121:   } else {
122:     VecSet(u, 1.0);
123:   }
124:   MatMult(A, u, b);

126:   /*
127:      View the exact solution vector if desired
128:   */
129:   flg = PETSC_FALSE;
130:   PetscOptionsGetBool(NULL, NULL, "-view_exact_sol", &flg, NULL);
131:   if (flg) VecView(u, PETSC_VIEWER_STDOUT_WORLD);

133:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134:                 Create the linear solver and set various options
135:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

137:   /*
138:      Create linear solver context
139:   */
140:   KSPCreate(PETSC_COMM_WORLD, &ksp);

142:   /*
143:      Set operators. Here the matrix that defines the linear system
144:      also serves as the preconditioning matrix.
145:   */
146:   KSPSetOperators(ksp, A, A);

148:   /*
149:      Set linear solver defaults for this problem (optional).
150:      - By extracting the KSP and PC contexts from the KSP context,
151:        we can then directly call any KSP and PC routines to set
152:        various options.
153:      - The following two statements are optional; all of these
154:        parameters could alternatively be specified at runtime via
155:        KSPSetFromOptions().  All of these defaults can be
156:        overridden at runtime, as indicated below.
157:   */
158:   KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);

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

169:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170:                       Solve the linear system
171:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

173:   KSPSolve(ksp, b, x);

175:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176:                       Check solution and clean up
177:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

186:   /*
187:      Print convergence information.  PetscPrintf() produces a single
188:      print statement from all processes that share a communicator.
189:      An alternative is PetscFPrintf(), which prints to a file.
190:   */
191:   PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its);

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

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

213: /*TEST

215:    test:
216:       nsize: 6
217:       args: -pc_type none
218:       requires: elemental

220:    test:
221:       suffix: 2
222:       nsize: 6
223:       args: -pc_type lu -pc_factor_mat_solver_type elemental
224:       requires: elemental

226: TEST*/