Actual source code: ex1.c


  2: static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";

  4: /*
  5:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  6:   automatically includes:
  7:      petscsys.h    - base PETSc routines   petscvec.h - vectors
  8:      petscmat.h    - matrices              petscpc.h  - preconditioners
  9:      petscis.h     - index sets
 10:      petscviewer.h - viewers

 12:   Note:  The corresponding parallel example is ex23.c
 13: */
 14: #include <petscksp.h>

 16: int main(int argc, char **args)
 17: {
 18:   Vec         x, b, u; /* approx solution, RHS, exact solution */
 19:   Mat         A;       /* linear system matrix */
 20:   KSP         ksp;     /* linear solver context */
 21:   PC          pc;      /* preconditioner context */
 22:   PetscReal   norm;    /* norm of solution error */
 23:   PetscInt    i, n = 10, col[3], its;
 24:   PetscMPIInt size;
 25:   PetscScalar value[3];

 27:   PetscFunctionBeginUser;
 28:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 29:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 30:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

 32:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

 34:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 35:          Compute the matrix and right-hand-side vector that define
 36:          the linear system, Ax = b.
 37:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 39:   /*
 40:      Create vectors.  Note that we form 1 vector from scratch and
 41:      then duplicate as needed.
 42:   */
 43:   PetscCall(VecCreate(PETSC_COMM_SELF, &x));
 44:   PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
 45:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
 46:   PetscCall(VecSetFromOptions(x));
 47:   PetscCall(VecDuplicate(x, &b));
 48:   PetscCall(VecDuplicate(x, &u));

 50:   /*
 51:      Create matrix.  When using MatCreate(), the matrix format can
 52:      be specified at runtime.

 54:      Performance tuning note:  For problems of substantial size,
 55:      preallocation of matrix memory is crucial for attaining good
 56:      performance. See the matrix chapter of the users manual for details.
 57:   */
 58:   PetscCall(MatCreate(PETSC_COMM_SELF, &A));
 59:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 60:   PetscCall(MatSetFromOptions(A));
 61:   PetscCall(MatSetUp(A));

 63:   /*
 64:      Assemble matrix
 65:   */
 66:   value[0] = -1.0;
 67:   value[1] = 2.0;
 68:   value[2] = -1.0;
 69:   for (i = 1; i < n - 1; i++) {
 70:     col[0] = i - 1;
 71:     col[1] = i;
 72:     col[2] = i + 1;
 73:     PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 74:   }
 75:   i      = n - 1;
 76:   col[0] = n - 2;
 77:   col[1] = n - 1;
 78:   PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
 79:   i        = 0;
 80:   col[0]   = 0;
 81:   col[1]   = 1;
 82:   value[0] = 2.0;
 83:   value[1] = -1.0;
 84:   PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
 85:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 86:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 88:   /*
 89:      Set exact solution; then compute right-hand-side vector.
 90:   */
 91:   PetscCall(VecSet(u, 1.0));
 92:   PetscCall(MatMult(A, u, b));

 94:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95:                 Create the linear solver and set various options
 96:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 97:   PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));

 99:   /*
100:      Set operators. Here the matrix that defines the linear system
101:      also serves as the matrix that defines the preconditioner.
102:   */
103:   PetscCall(KSPSetOperators(ksp, A, A));

105:   /*
106:      Set linear solver defaults for this problem (optional).
107:      - By extracting the KSP and PC contexts from the KSP context,
108:        we can then directly call any KSP and PC routines to set
109:        various options.
110:      - The following four statements are optional; all of these
111:        parameters could alternatively be specified at runtime via
112:        KSPSetFromOptions();
113:   */
114:   PetscCall(KSPGetPC(ksp, &pc));
115:   PetscCall(PCSetType(pc, PCJACOBI));
116:   PetscCall(KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));

118:   /*
119:     Set runtime options, e.g.,
120:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
121:     These options will override those specified above as long as
122:     KSPSetFromOptions() is called _after_ any other customization
123:     routines.
124:   */
125:   PetscCall(KSPSetFromOptions(ksp));

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:                       Solve the linear system
129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130:   PetscCall(KSPSolve(ksp, b, x));

132:   /*
133:      View solver info; we could instead use the option -ksp_view to
134:      print this info to the screen at the conclusion of KSPSolve().
135:   */
136:   PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF));

138:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139:                       Check the solution and clean up
140:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141:   PetscCall(VecAXPY(x, -1.0, u));
142:   PetscCall(VecNorm(x, NORM_2, &norm));
143:   PetscCall(KSPGetIterationNumber(ksp, &its));
144:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));

146:   /* check that KSP automatically handles the fact that the the new non-zero values in the matrix are propagated to the KSP solver */
147:   PetscCall(MatShift(A, 2.0));
148:   PetscCall(KSPSolve(ksp, b, x));

150:   /*
151:      Free work space.  All PETSc objects should be destroyed when they
152:      are no longer needed.
153:   */
154:   PetscCall(VecDestroy(&x));
155:   PetscCall(VecDestroy(&u));
156:   PetscCall(VecDestroy(&b));
157:   PetscCall(MatDestroy(&A));
158:   PetscCall(KSPDestroy(&ksp));

160:   /*
161:      Always call PetscFinalize() before exiting a program.  This routine
162:        - finalizes the PETSc libraries as well as MPI
163:        - provides summary and diagnostic information if certain runtime
164:          options are chosen (e.g., -log_view).
165:   */
166:   PetscCall(PetscFinalize());
167:   return 0;
168: }

170: /*TEST

172:    test:
173:       args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

175:    test:
176:       suffix: 2
177:       args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

179:    test:
180:       suffix: 2_aijcusparse
181:       requires: cuda
182:       args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -mat_type aijcusparse -vec_type cuda

184:    test:
185:       suffix: 3
186:       args: -pc_type eisenstat -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

188:    test:
189:       suffix: 3_aijcusparse
190:       requires: cuda
191:       args: -pc_type eisenstat -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -mat_type aijcusparse -vec_type cuda

193:    test:
194:       suffix: aijcusparse
195:       requires: cuda
196:       args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -mat_type aijcusparse -vec_type cuda
197:       output_file: output/ex1_1_aijcusparse.out

199:    test:
200:       requires: defined(PETSC_USE_SINGLE_LIBRARY)
201:       suffix: mpi_linear_solver_server_1
202:       nsize: 3
203:       filter: sed 's?ATOL?RTOL?g'
204:       args: -mpi_linear_solver_server -mpi_linear_solver_server_view -pc_type mpi -ksp_type preonly -mpi_ksp_monitor -mpi_ksp_converged_reason -mat_view -mpi_pc_type none -mpi_ksp_view -mpi_mat_view -pc_mpi_minimum_count_per_rank 5

206:    test:
207:       requires: defined(PETSC_USE_SINGLE_LIBRARY)
208:       suffix: mpi_linear_solver_server_2
209:       nsize: 3
210:       filter: sed 's?ATOL?RTOL?g'
211:       args: -mpi_linear_solver_server  -mpi_linear_solver_server_view -pc_type mpi -ksp_type preonly -mpi_ksp_monitor -mpi_ksp_converged_reason -mat_view -mpi_pc_type none -mpi_ksp_view

213:    test:
214:       requires: defined(PETSC_USE_SINGLE_LIBRARY)
215:       suffix: mpi_linear_solver_server_3
216:       nsize: 3
217:       filter: sed 's?ATOL?RTOL?g'
218:       args: -mpi_linear_solver_server  -mpi_linear_solver_server_view -pc_type mpi -ksp_type preonly -mpi_ksp_monitor -mpi_ksp_converged_reason -mat_view -mpi_pc_type none -mpi_ksp_view -mpi_mat_view -pc_mpi_always_use_server

220: TEST*/