Actual source code: ex3.c
1: static char help[] = "Demonstrates the use of fast Richardson for SOR. And\n\
2: also tests the MatSOR() routines. Input parameters are:\n\
3: -n <n> : problem dimension\n\n";
5: #include <petscksp.h>
6: #include <petscpc.h>
8: int main(int argc, char **args)
9: {
10: Mat mat; /* matrix */
11: Vec b, ustar, u; /* vectors (RHS, exact solution, approx solution) */
12: PC pc; /* PC context */
13: KSP ksp; /* KSP context */
14: PetscInt n = 10, i, its, col[3];
15: PetscScalar value[3];
16: KSPType kspname;
17: PCType pcname;
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &args, NULL, help));
21: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
23: /* Create and initialize vectors */
24: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &b));
25: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &ustar));
26: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &u));
27: PetscCall(VecSet(ustar, 1.0));
28: PetscCall(VecSet(u, 0.0));
30: /* Create and assemble matrix */
31: PetscCall(MatCreate(PETSC_COMM_SELF, &mat));
32: PetscCall(MatSetType(mat, MATSEQAIJ));
33: PetscCall(MatSetSizes(mat, n, n, n, n));
34: PetscCall(MatSetFromOptions(mat));
35: PetscCall(MatSeqAIJSetPreallocation(mat, 3, NULL));
36: PetscCall(MatSeqBAIJSetPreallocation(mat, 1, 3, NULL));
37: PetscCall(MatSeqSBAIJSetPreallocation(mat, 1, 3, NULL));
38: PetscCall(MatSeqSELLSetPreallocation(mat, 3, NULL));
39: value[0] = -1.0;
40: value[1] = 2.0;
41: value[2] = -1.0;
42: for (i = 1; i < n - 1; i++) {
43: col[0] = i - 1;
44: col[1] = i;
45: col[2] = i + 1;
46: PetscCall(MatSetValues(mat, 1, &i, 3, col, value, INSERT_VALUES));
47: }
48: i = n - 1;
49: col[0] = n - 2;
50: col[1] = n - 1;
51: PetscCall(MatSetValues(mat, 1, &i, 2, col, value, INSERT_VALUES));
52: i = 0;
53: col[0] = 0;
54: col[1] = 1;
55: value[0] = 2.0;
56: value[1] = -1.0;
57: PetscCall(MatSetValues(mat, 1, &i, 2, col, value, INSERT_VALUES));
58: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
59: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
61: /* Compute right-hand-side vector */
62: PetscCall(MatMult(mat, ustar, b));
64: /* Create PC context and set up data structures */
65: PetscCall(PCCreate(PETSC_COMM_WORLD, &pc));
66: PetscCall(PCSetType(pc, PCNONE));
67: PetscCall(PCSetFromOptions(pc));
68: PetscCall(PCSetOperators(pc, mat, mat));
69: PetscCall(PCSetUp(pc));
71: /* Create KSP context and set up data structures */
72: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
73: PetscCall(KSPSetType(ksp, KSPRICHARDSON));
74: PetscCall(KSPSetFromOptions(ksp));
75: PetscCall(PCSetOperators(pc, mat, mat));
76: PetscCall(KSPSetPC(ksp, pc));
77: PetscCall(KSPSetUp(ksp));
79: /* Solve the problem */
80: PetscCall(KSPGetType(ksp, &kspname));
81: PetscCall(PCGetType(pc, &pcname));
82: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Running %s with %s preconditioning\n", kspname, pcname));
83: PetscCall(KSPSolve(ksp, b, u));
84: PetscCall(KSPGetIterationNumber(ksp, &its));
85: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations %" PetscInt_FMT "\n", its));
87: /* Free data structures */
88: PetscCall(KSPDestroy(&ksp));
89: PetscCall(VecDestroy(&u));
90: PetscCall(VecDestroy(&ustar));
91: PetscCall(VecDestroy(&b));
92: PetscCall(MatDestroy(&mat));
93: PetscCall(PCDestroy(&pc));
94: PetscCall(PetscFinalize());
95: return 0;
96: }
98: /*TEST
100: testset:
101: args: -ksp_type gmres -ksp_monitor_short -pc_type sor -pc_sor_symmetric
102: output_file: output/ex3_1.out
103: test:
104: suffix: sor_aij
105: test:
106: suffix: sor_seqbaij
107: args: -mat_type seqbaij
108: test:
109: suffix: sor_seqsbaij
110: args: -mat_type seqbaij
111: test:
112: suffix: sor_seqsell
113: args: -mat_type seqsell
115: TEST*/