Actual source code: ex24.c

  1: static char help[] = "Tests CG, MINRES and SYMMLQ on symmetric matrices with SBAIJ format. The preconditioner ICC only works on sequential SBAIJ format. \n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat         C;
  8:   PetscScalar v, none = -1.0;
  9:   PetscInt    i, j, Ii, J, Istart, Iend, N, m = 4, n = 4, its, k;
 10:   PetscMPIInt size, rank;
 11:   PetscReal   err_norm, res_norm;
 12:   Vec         x, b, u, u_tmp;
 13:   PC          pc;
 14:   KSP         ksp;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 18:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 19:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 20:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 21:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 22:   N = m * n;

 24:   /* Generate matrix */
 25:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 26:   PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N));
 27:   PetscCall(MatSetFromOptions(C));
 28:   PetscCall(MatSetUp(C));
 29:   PetscCall(MatGetOwnershipRange(C, &Istart, &Iend));
 30:   for (Ii = Istart; Ii < Iend; Ii++) {
 31:     v = -1.0;
 32:     i = Ii / n;
 33:     j = Ii - i * n;
 34:     if (i > 0) {
 35:       J = Ii - n;
 36:       PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
 37:     }
 38:     if (i < m - 1) {
 39:       J = Ii + n;
 40:       PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
 41:     }
 42:     if (j > 0) {
 43:       J = Ii - 1;
 44:       PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
 45:     }
 46:     if (j < n - 1) {
 47:       J = Ii + 1;
 48:       PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
 49:     }
 50:     v = 4.0;
 51:     PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
 52:   }
 53:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 54:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

 56:   /* a shift can make C indefinite. Preconditioners LU, ILU (for BAIJ format) and ICC may fail */
 57:   /* PetscCall(MatShift(C,alpha)); */
 58:   /* PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); */

 60:   /* Setup and solve for system */
 61:   /* Create vectors.  */
 62:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 63:   PetscCall(VecSetSizes(x, PETSC_DECIDE, N));
 64:   PetscCall(VecSetFromOptions(x));
 65:   PetscCall(VecDuplicate(x, &b));
 66:   PetscCall(VecDuplicate(x, &u));
 67:   PetscCall(VecDuplicate(x, &u_tmp));
 68:   /* Set exact solution u; then compute right-hand-side vector b. */
 69:   PetscCall(VecSet(u, 1.0));
 70:   PetscCall(MatMult(C, u, b));

 72:   for (k = 0; k < 3; k++) {
 73:     if (k == 0) { /* CG  */
 74:       PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 75:       PetscCall(KSPSetOperators(ksp, C, C));
 76:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n CG: \n"));
 77:       PetscCall(KSPSetType(ksp, KSPCG));
 78:     } else if (k == 1) { /* MINRES */
 79:       PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 80:       PetscCall(KSPSetOperators(ksp, C, C));
 81:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MINRES: \n"));
 82:       PetscCall(KSPSetType(ksp, KSPMINRES));
 83:     } else { /* SYMMLQ */
 84:       PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 85:       PetscCall(KSPSetOperators(ksp, C, C));
 86:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n SYMMLQ: \n"));
 87:       PetscCall(KSPSetType(ksp, KSPSYMMLQ));
 88:     }
 89:     PetscCall(KSPGetPC(ksp, &pc));
 90:     /* PetscCall(PCSetType(pc,PCICC)); */
 91:     PetscCall(PCSetType(pc, PCJACOBI));
 92:     PetscCall(KSPSetTolerances(ksp, 1.e-7, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));

 94:     /*
 95:     Set runtime options, e.g.,
 96:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
 97:     These options will override those specified above as long as
 98:     KSPSetFromOptions() is called _after_ any other customization
 99:     routines.
100:     */
101:     PetscCall(KSPSetFromOptions(ksp));

103:     /* Solve linear system; */
104:     PetscCall(KSPSetUp(ksp));
105:     PetscCall(KSPSolve(ksp, b, x));

107:     PetscCall(KSPGetIterationNumber(ksp, &its));
108:     /* Check error */
109:     PetscCall(VecCopy(u, u_tmp));
110:     PetscCall(VecAXPY(u_tmp, none, x));
111:     PetscCall(VecNorm(u_tmp, NORM_2, &err_norm));
112:     PetscCall(MatMult(C, x, u_tmp));
113:     PetscCall(VecAXPY(u_tmp, none, b));
114:     PetscCall(VecNorm(u_tmp, NORM_2, &res_norm));

116:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
117:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g;", (double)res_norm));
118:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Error norm %g.\n", (double)err_norm));
119:     PetscCall(KSPDestroy(&ksp));
120:   }

122:   /*
123:        Free work space.  All PETSc objects should be destroyed when they
124:        are no longer needed.
125:   */
126:   PetscCall(VecDestroy(&b));
127:   PetscCall(VecDestroy(&u));
128:   PetscCall(VecDestroy(&x));
129:   PetscCall(VecDestroy(&u_tmp));
130:   PetscCall(MatDestroy(&C));

132:   PetscCall(PetscFinalize());
133:   return 0;
134: }

136: /*TEST

138:     test:
139:       args: -pc_type icc -mat_type seqsbaij -mat_ignore_lower_triangular

141:     test:
142:       suffix: 2
143:       args: -pc_type icc -pc_factor_levels 2 -mat_type seqsbaij -mat_ignore_lower_triangular

145:     test:
146:       suffix: 3
147:       nsize: 2
148:       args: -pc_type bjacobi -sub_pc_type icc -mat_type mpisbaij -mat_ignore_lower_triangular -ksp_max_it 8

150:     test:
151:       suffix: 4
152:       nsize: 2
153:       args: -pc_type bjacobi -sub_pc_type icc -sub_pc_factor_levels 1 -mat_type mpisbaij -mat_ignore_lower_triangular

155: TEST*/