Actual source code: ex50.c

  1: static char help[] = "Tests point block Jacobi and ILU for different block sizes\n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   Vec         x, b, u;
  8:   Mat         A;    /* linear system matrix */
  9:   KSP         ksp;  /* linear solver context */
 10:   PetscRandom rctx; /* random number generator context */
 11:   PetscReal   norm; /* norm of solution error */
 12:   PetscInt    i, j, k, l, n = 27, its, bs = 2, Ii, J;
 13:   PetscScalar v;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 17:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
 18:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

 20:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 21:   PetscCall(MatSetSizes(A, n * bs, n * bs, PETSC_DETERMINE, PETSC_DETERMINE));
 22:   PetscCall(MatSetBlockSize(A, bs));
 23:   PetscCall(MatSetFromOptions(A));
 24:   PetscCall(MatSetUp(A));

 26:   /*
 27:      Don't bother to preallocate matrix
 28:   */
 29:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rctx));
 30:   for (i = 0; i < n; i++) {
 31:     for (j = 0; j < n; j++) {
 32:       PetscCall(PetscRandomGetValue(rctx, &v));
 33:       if (PetscRealPart(v) < .25 || i == j) {
 34:         for (k = 0; k < bs; k++) {
 35:           for (l = 0; l < bs; l++) {
 36:             PetscCall(PetscRandomGetValue(rctx, &v));
 37:             Ii = i * bs + k;
 38:             J  = j * bs + l;
 39:             if (Ii == J) v += 10.;
 40:             PetscCall(MatSetValue(A, Ii, J, v, INSERT_VALUES));
 41:           }
 42:         }
 43:       }
 44:     }
 45:   }

 47:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 48:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 49:   PetscCall(MatCreateVecs(A, &u, &b));
 50:   PetscCall(VecDuplicate(u, &x));
 51:   PetscCall(VecSet(u, 1.0));
 52:   PetscCall(MatMult(A, u, b));

 54:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55:                 Create the linear solver and set various options
 56:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 58:   /*
 59:      Create linear solver context
 60:   */
 61:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

 63:   /*
 64:      Set operators. Here the matrix that defines the linear system
 65:      also serves as the preconditioning matrix.
 66:   */
 67:   PetscCall(KSPSetOperators(ksp, A, A));

 69:   PetscCall(KSPSetFromOptions(ksp));

 71:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 72:                       Solve the linear system
 73:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 75:   PetscCall(KSPSolve(ksp, b, x));

 77:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 78:                       Check solution and clean up
 79:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 81:   /*
 82:      Check the error
 83:   */
 84:   PetscCall(VecAXPY(x, -1.0, u));
 85:   PetscCall(VecNorm(x, NORM_2, &norm));
 86:   PetscCall(KSPGetIterationNumber(ksp, &its));

 88:   /*
 89:      Print convergence information.  PetscPrintf() produces a single
 90:      print statement from all processes that share a communicator.
 91:      An alternative is PetscFPrintf(), which prints to a file.
 92:   */
 93:   if (norm > .1) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of residual %g iterations %" PetscInt_FMT " bs %" PetscInt_FMT "\n", (double)norm, its, bs));

 95:   /*
 96:      Free work space.  All PETSc objects should be destroyed when they
 97:      are no longer needed.
 98:   */
 99:   PetscCall(KSPDestroy(&ksp));
100:   PetscCall(VecDestroy(&u));
101:   PetscCall(VecDestroy(&x));
102:   PetscCall(VecDestroy(&b));
103:   PetscCall(MatDestroy(&A));
104:   PetscCall(PetscRandomDestroy(&rctx));

106:   /*
107:      Always call PetscFinalize() before exiting a program.  This routine
108:        - finalizes the PETSc libraries as well as MPI
109:        - provides summary and diagnostic information if certain runtime
110:          options are chosen (e.g., -log_view).
111:   */
112:   PetscCall(PetscFinalize());
113:   return 0;
114: }

116: /*TEST

118:   testset:
119:     args: -bs {{1 2 3 4 5 6 7 8 11 15}} -pc_type {{pbjacobi ilu}}
120:     output_file: output/ex50_1.out

122:     test:
123:       args: -mat_type {{aij baij}}

125:     test:
126:       suffix: cuda
127:       requires: cuda
128:       args: -mat_type aijcusparse

130:     test:
131:       suffix: kok
132:       requires: kokkos_kernels
133:       args: -mat_type aijkokkos

135: TEST*/