Actual source code: ex50.c
2: static char help[] = "Tests point block Jacobi and ILU for different block sizes\n\n";
4: #include <petscksp.h>
6: int main(int argc, char **args)
7: {
8: Vec x, b, u;
9: Mat A; /* linear system matrix */
10: KSP ksp; /* linear solver context */
11: PetscRandom rctx; /* random number generator context */
12: PetscReal norm; /* norm of solution error */
13: PetscInt i, j, k, l, n = 27, its, bs = 2, Ii, J;
14: PetscScalar v;
17: PetscInitialize(&argc, &args, (char *)0, help);
18: PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL);
19: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
21: MatCreate(PETSC_COMM_WORLD, &A);
22: MatSetSizes(A, n * bs, n * bs, PETSC_DETERMINE, PETSC_DETERMINE);
23: MatSetBlockSize(A, bs);
24: MatSetFromOptions(A);
25: MatSetUp(A);
27: /*
28: Don't bother to preallocate matrix
29: */
30: PetscRandomCreate(PETSC_COMM_SELF, &rctx);
31: for (i = 0; i < n; i++) {
32: for (j = 0; j < n; j++) {
33: PetscRandomGetValue(rctx, &v);
34: if (PetscRealPart(v) < .25 || i == j) {
35: for (k = 0; k < bs; k++) {
36: for (l = 0; l < bs; l++) {
37: PetscRandomGetValue(rctx, &v);
38: Ii = i * bs + k;
39: J = j * bs + l;
40: if (Ii == J) v += 10.;
41: MatSetValue(A, Ii, J, v, INSERT_VALUES);
42: }
43: }
44: }
45: }
46: }
48: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
49: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
50: MatCreateVecs(A, &u, &b);
51: VecDuplicate(u, &x);
52: VecSet(u, 1.0);
53: MatMult(A, u, b);
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Create the linear solver and set various options
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: /*
60: Create linear solver context
61: */
62: KSPCreate(PETSC_COMM_WORLD, &ksp);
64: /*
65: Set operators. Here the matrix that defines the linear system
66: also serves as the preconditioning matrix.
67: */
68: KSPSetOperators(ksp, A, A);
70: KSPSetFromOptions(ksp);
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Solve the linear system
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: KSPSolve(ksp, b, x);
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Check solution and clean up
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82: /*
83: Check the error
84: */
85: VecAXPY(x, -1.0, u);
86: VecNorm(x, NORM_2, &norm);
87: KSPGetIterationNumber(ksp, &its);
89: /*
90: Print convergence information. PetscPrintf() produces a single
91: print statement from all processes that share a communicator.
92: An alternative is PetscFPrintf(), which prints to a file.
93: */
94: if (norm > .1) PetscPrintf(PETSC_COMM_WORLD, "Norm of residual %g iterations %" PetscInt_FMT " bs %" PetscInt_FMT "\n", (double)norm, its, bs);
96: /*
97: Free work space. All PETSc objects should be destroyed when they
98: are no longer needed.
99: */
100: KSPDestroy(&ksp);
101: VecDestroy(&u);
102: VecDestroy(&x);
103: VecDestroy(&b);
104: MatDestroy(&A);
105: PetscRandomDestroy(&rctx);
107: /*
108: Always call PetscFinalize() before exiting a program. This routine
109: - finalizes the PETSc libraries as well as MPI
110: - provides summary and diagnostic information if certain runtime
111: options are chosen (e.g., -log_view).
112: */
113: PetscFinalize();
114: return 0;
115: }
117: /*TEST
118: testset:
119: args: -mat_type {{aij baij}}
121: test:
122: args: -bs {{1 2 3 4 5 6 7 8}} -pc_type pbjacobi
124: test:
125: suffix: 2
126: args: -bs {{8 9 10 11 12 13 14 15}} -pc_type ilu
128: TEST*/