Actual source code: ex85.c

  1: static char help[] = "Estimate eigenvalues with KSP.\n\n";

  3: /*
  4:     Test example that demonstrates how KSP can estimate eigenvalues.

  6:     Contributed by: Pablo Brubeck <brubeck@protonmail.com>
  7: */
  8: #include <petscksp.h>

 10: int main(int argc, char **args)
 11: {
 12:   Vec         x, b; /* approx solution, RHS */
 13:   Mat         A;    /* linear system matrix */
 14:   KSP         ksp;  /* linear solver context */
 15:   PC          pc;   /* preconditioner context */
 16:   PetscInt    n = 6, i, j, col[2];
 17:   PetscScalar value[4];
 18:   PetscMPIInt size;

 20:   PetscFunctionBeginUser;
 21:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 22:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 23:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

 25:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 26:          Compute the matrix and right-hand-side vector that define
 27:          the linear system, Ax = b.
 28:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 30:   /*
 31:      Create vectors.  Note that we form 1 vector from scratch and
 32:      then duplicate as needed.
 33:   */
 34:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 35:   PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
 36:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
 37:   PetscCall(VecSetFromOptions(x));
 38:   PetscCall(VecDuplicate(x, &b));

 40:   /*
 41:      Create matrix.  When using MatCreate(), the matrix format can
 42:      be specified at runtime.

 44:      Performance tuning note:  For problems of substantial size,
 45:      preallocation of matrix memory is crucial for attaining good
 46:      performance. See the matrix chapter of the users manual for details.
 47:   */
 48:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 49:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 50:   PetscCall(MatSetFromOptions(A));
 51:   PetscCall(MatSetUp(A));

 53:   /*
 54:      Assemble matrix
 55:   */
 56:   value[0] = 2.0;
 57:   value[1] = -1.0;
 58:   value[2] = -1.0;
 59:   value[3] = 2.0;
 60:   for (i = 0; 2 * i < n; i++) {
 61:     col[0] = 2 * i;
 62:     col[1] = 2 * i + 1;
 63:     PetscCall(MatSetValues(A, 2, col, 2, col, value, INSERT_VALUES));
 64:     for (j = 0; j < 4; j++) value[j] *= 3.0;
 65:   }
 66:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 67:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 69:   /*
 70:      Set random right-hand-side vector.
 71:   */
 72:   PetscCall(VecSetRandom(b, NULL));

 74:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 75:                 Create the linear solver and set various options
 76:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 77:   /*
 78:      Create linear solver context
 79:   */
 80:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

 82:   /*
 83:      Set operators. Here the matrix that defines the linear system
 84:      also serves as the preconditioning matrix.
 85:   */
 86:   PetscCall(KSPSetOperators(ksp, A, A));

 88:   /*
 89:      Set linear solver defaults for this problem (optional).
 90:      - By extracting the KSP and PC contexts from the KSP context,
 91:        we can then directly call any KSP and PC routines to set
 92:        various options.
 93:      - The following four statements are optional; all of these
 94:        parameters could alternatively be specified at runtime via
 95:        KSPSetFromOptions();
 96:   */
 97:   PetscCall(KSPGetPC(ksp, &pc));
 98:   PetscCall(PCSetType(pc, PCJACOBI));
 99:   PetscCall(KSPSetTolerances(ksp, 1.e-5, 1.e-5, PETSC_CURRENT, PETSC_CURRENT));

101:   /*
102:     Set runtime options, e.g.,
103:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
104:     These options will override those specified above as long as
105:     KSPSetFromOptions() is called _after_ any other customization
106:     routines.
107:   */
108:   PetscCall(KSPSetFromOptions(ksp));

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:                       Solve the linear system
112:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113:   /*
114:      Solve linear system
115:   */
116:   PetscCall(KSPSolve(ksp, b, x));

118:   /*
119:      View solver info; we could instead use the option -ksp_view to
120:      print this info to the screen at the conclusion of KSPSolve().
121:   */
122:   PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));

124:   /*
125:      Free work space.  All PETSc objects should be destroyed when they
126:      are no longer needed.
127:   */
128:   PetscCall(VecDestroy(&x));
129:   PetscCall(VecDestroy(&b));
130:   PetscCall(MatDestroy(&A));
131:   PetscCall(KSPDestroy(&ksp));

133:   /*
134:      Always call PetscFinalize() before exiting a program.  This routine
135:        - finalizes the PETSc libraries as well as MPI
136:        - provides summary and diagnostic information if certain runtime
137:          options are chosen (e.g., -log_view).
138:   */
139:   PetscCall(PetscFinalize());
140:   return 0;
141: }

143: /*TEST

145:    test:
146:       suffix: cg_none
147:       filter: grep -v "variant HERMITIAN"
148:       args: -ksp_type cg -pc_type none -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

150:    test:
151:       suffix: cg_jacobi
152:       filter: grep -v "variant HERMITIAN"
153:       args: -ksp_type cg -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

155:    test:
156:       suffix: fcg_none
157:       args: -ksp_type fcg -pc_type none -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

159:    test:
160:       suffix: fcg_jacobi
161:       args: -ksp_type fcg -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

163:    test:
164:       suffix: minres_none
165:       args: -ksp_type minres -pc_type none -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

167:    test:
168:       suffix: minres_jacobi
169:       args: -ksp_type minres -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

171:    test:
172:       suffix: gmres_none
173:       args: -ksp_type gmres -pc_type none -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

175:    test:
176:       suffix: gmres_jacobi_left
177:       args: -ksp_type gmres -ksp_pc_side left -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

179:    test:
180:       suffix: gmres_jacobi_right
181:       args: -ksp_type gmres -ksp_pc_side right -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

183:    test:
184:       suffix: fgmres_none
185:       args: -ksp_type fgmres -pc_type none -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

187:    test:
188:       suffix: fgmres_jacobi
189:       args: -ksp_type fgmres -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason

191: TEST*/