Actual source code: ex53.c
2: static char help[] = "Solves a tridiagonal linear system with KSP. \n\
3: Modified from ex1.c to illustrate reuse of preconditioner \n\
4: Written as requested by [petsc-maint #63875] \n\n";
6: #include <petscksp.h>
8: int main(int argc, char **args)
9: {
10: Vec x, x2, b, u; /* approx solution, RHS, exact solution */
11: Mat A; /* linear system matrix */
12: KSP ksp; /* linear solver context */
13: PC pc; /* preconditioner context */
14: PetscReal norm, tol = 100. * PETSC_MACHINE_EPSILON; /* norm of solution error */
15: PetscInt i, n = 10, col[3], its;
16: PetscMPIInt rank, size;
17: PetscScalar one = 1.0, value[3];
20: PetscInitialize(&argc, &args, (char *)0, help);
21: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
22: MPI_Comm_size(PETSC_COMM_WORLD, &size);
23: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
25: /* Create vectors.*/
26: VecCreate(PETSC_COMM_WORLD, &x);
27: PetscObjectSetName((PetscObject)x, "Solution");
28: VecSetSizes(x, PETSC_DECIDE, n);
29: VecSetFromOptions(x);
30: VecDuplicate(x, &b);
31: VecDuplicate(x, &u);
32: VecDuplicate(x, &x2);
34: /* Create matrix. Only proc[0] sets values - not efficient for parallel processing!
35: See ex23.c for efficient parallel assembly matrix */
36: MatCreate(PETSC_COMM_WORLD, &A);
37: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
38: MatSetFromOptions(A);
39: MatSetUp(A);
41: if (rank == 0) {
42: value[0] = -1.0;
43: value[1] = 2.0;
44: value[2] = -1.0;
45: for (i = 1; i < n - 1; i++) {
46: col[0] = i - 1;
47: col[1] = i;
48: col[2] = i + 1;
49: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
50: }
51: i = n - 1;
52: col[0] = n - 2;
53: col[1] = n - 1;
54: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
55: i = 0;
56: col[0] = 0;
57: col[1] = 1;
58: value[0] = 2.0;
59: value[1] = -1.0;
60: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
62: i = 0;
63: col[0] = n - 1;
64: value[0] = 0.5; /* make A non-symmetric */
65: MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES);
66: }
67: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
68: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
70: /* Set exact solution */
71: VecSet(u, one);
73: /* Create linear solver context */
74: KSPCreate(PETSC_COMM_WORLD, &ksp);
75: KSPSetOperators(ksp, A, A);
76: KSPGetPC(ksp, &pc);
77: PCSetType(pc, PCLU);
78: #if defined(PETSC_HAVE_MUMPS)
79: if (size > 1) PCFactorSetMatSolverType(pc, MATSOLVERMUMPS);
80: #endif
81: KSPSetFromOptions(ksp);
83: /* 1. Solve linear system A x = b */
84: MatMult(A, u, b);
85: KSPSolve(ksp, b, x);
87: /* Check the error */
88: VecAXPY(x, -1.0, u);
89: VecNorm(x, NORM_2, &norm);
90: KSPGetIterationNumber(ksp, &its);
91: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "1. Norm of error for Ax=b: %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
93: /* 2. Solve linear system A^T x = b*/
94: MatMultTranspose(A, u, b);
95: KSPSolveTranspose(ksp, b, x2);
97: /* Check the error */
98: VecAXPY(x2, -1.0, u);
99: VecNorm(x2, NORM_2, &norm);
100: KSPGetIterationNumber(ksp, &its);
101: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "2. Norm of error for A^T x=b: %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
103: /* 3. Change A and solve A x = b with an iterative solver using A=LU as a preconditioner*/
104: if (rank == 0) {
105: i = 0;
106: col[0] = n - 1;
107: value[0] = 1.e-2;
108: MatSetValues(A, 1, &i, 1, col, value, ADD_VALUES);
109: }
110: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
111: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
113: MatMult(A, u, b);
114: KSPSolve(ksp, b, x);
116: /* Check the error */
117: VecAXPY(x, -1.0, u);
118: VecNorm(x, NORM_2, &norm);
119: KSPGetIterationNumber(ksp, &its);
120: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "3. Norm of error for (A+Delta) x=b: %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
122: /* Free work space. */
123: VecDestroy(&x);
124: VecDestroy(&u);
125: VecDestroy(&x2);
126: VecDestroy(&b);
127: MatDestroy(&A);
128: KSPDestroy(&ksp);
130: PetscFinalize();
131: return 0;
132: }
134: /*TEST
136: test:
137: requires: mumps
139: test:
140: suffix: 2
141: nsize: 2
142: requires: mumps
143: output_file: output/ex53.out
145: TEST*/