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*/