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;

 25:   PetscInitialize(&argc, &args, (char *)0, help);
 26:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 28:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 29:   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:   VecCreate(PETSC_COMM_WORLD, &x);
 41:   PetscObjectSetName((PetscObject)x, "Solution");
 42:   VecSetSizes(x, PETSC_DECIDE, n);
 43:   VecSetFromOptions(x);
 44:   VecDuplicate(x, &b);
 45:   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:   MatCreate(PETSC_COMM_WORLD, &A);
 56:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
 57:   MatSetFromOptions(A);
 58:   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:     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:   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:   MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
 82:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 83:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

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

 94:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95:                 Create the linear solver and set various options
 96:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 97:   /*
 98:      Create linear solver context
 99:   */
100:   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:   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:   KSPGetPC(ksp, &pc);
118:   PCSetType(pc, PCJACOBI);
119:   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:   KSPSetFromOptions(ksp);

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

136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137:                       Solve the linear system
138:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139:   /*
140:      Solve linear system
141:   */
142:   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:   KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD);

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:                       Check solution and clean up
152:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153:   /*
154:      Check the error
155:   */
156:   VecAXPY(x, -1.0, u);
157:   VecNorm(x, NORM_2, &norm);
158:   KSPGetIterationNumber(ksp, &its);
159:   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:   VecDestroy(&x);
166:   VecDestroy(&u);
167:   VecDestroy(&b);
168:   MatDestroy(&A);
169:   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:   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*/