Actual source code: ex48.c

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

  3: /*
  4:     Test example that demonstrates how MINRES can produce a dp of zero
  5:     but still converge.

  7:     Provided by: Mark Filipiak <mjf@staffmail.ed.ac.uk>
  8: */
  9: #include <petscksp.h>

 11: int main(int argc, char **args)
 12: {
 13:   Vec         x, b, u; /* approx solution, RHS, exact solution */
 14:   Mat         A;       /* linear system matrix */
 15:   KSP         ksp;     /* linear solver context */
 16:   PC          pc;      /* preconditioner context */
 17:   PetscReal   norm;
 18:   PetscInt    i, n = 2, col[3], its;
 19:   PetscMPIInt size;
 20:   PetscScalar one = 1.0, value[3], shift = 0.0;
 21:   PetscBool   nonzeroguess = PETSC_FALSE;

 23:   PetscFunctionBeginUser;
 24:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 25:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 26:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 27:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 28:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-nonzero_guess", &nonzeroguess, NULL));
 29:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-shift", &shift, NULL));

 31:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 32:          Compute the matrix and right-hand-side vector that define
 33:          the linear system, Ax = b.
 34:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

 47:   /*
 48:      Create matrix.  When using MatCreate(), the matrix format can
 49:      be specified at runtime.

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

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

 86:   /*
 87:      Set constant right-hand-side vector.
 88:   */
 89:   PetscCall(VecSet(b, one));
 90:   /*
 91:      Solution = RHS for the matrix [[2 -1] [-1 2]] and RHS [1 1]
 92:   */
 93:   PetscCall(VecSet(u, one));

 95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:                 Create the linear solver and set various options
 97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 98:   /*
 99:      Create linear solver context
100:   */
101:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

103:   /*
104:      Set operators. Here the matrix that defines the linear system
105:      also serves as the preconditioning matrix.
106:   */
107:   PetscCall(KSPSetOperators(ksp, A, A));

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

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

131:   if (nonzeroguess) {
132:     PetscScalar p = .5;
133:     PetscCall(VecSet(x, p));
134:     PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
135:   }

137:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138:                       Solve the linear system
139:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140:   /*
141:      Solve linear system
142:   */
143:   PetscCall(KSPSolve(ksp, b, x));

145:   /*
146:      View solver info; we could instead use the option -ksp_view to
147:      print this info to the screen at the conclusion of KSPSolve().
148:   */
149:   PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));

151:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152:                       Check solution and clean up
153:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154:   /*
155:      Check the error
156:   */
157:   PetscCall(VecAXPY(x, -1.0, u));
158:   PetscCall(VecNorm(x, NORM_2, &norm));
159:   PetscCall(KSPGetIterationNumber(ksp, &its));
160:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));

162:   /*
163:      Free work space.  All PETSc objects should be destroyed when they
164:      are no longer needed.
165:   */
166:   PetscCall(VecDestroy(&x));
167:   PetscCall(VecDestroy(&u));
168:   PetscCall(VecDestroy(&b));
169:   PetscCall(MatDestroy(&A));
170:   PetscCall(KSPDestroy(&ksp));

172:   /*
173:      Always call PetscFinalize() before exiting a program.  This routine
174:        - finalizes the PETSc libraries as well as MPI
175:        - provides summary and diagnostic information if certain runtime
176:          options are chosen (e.g., -log_view).
177:   */
178:   PetscCall(PetscFinalize());
179:   return 0;
180: }

182: /*TEST

184:    test:
185:       args: -ksp_monitor_short -ksp_converged_reason -ksp_type minres -pc_type jacobi -ksp_error_if_not_converged

187:    test:
188:       suffix: 2
189:       args: -ksp_monitor_short -ksp_converged_reason -ksp_type minres -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -shift 1.5 -ksp_rtol 0 -ksp_atol 0 -n 5 -ksp_max_it 5

191: TEST*/