Actual source code: ex48.c


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

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

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

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

 24:   PetscFunctionBeginUser;
 25:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 26:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 27:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 28:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 29:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-nonzero_guess", &nonzeroguess, 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));

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

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

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

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

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

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

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

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

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

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

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

181: /*TEST

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

186: TEST*/