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];
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
21: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
22: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
23: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
25: /* Create vectors.*/
26: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
27: PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
28: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
29: PetscCall(VecSetFromOptions(x));
30: PetscCall(VecDuplicate(x, &b));
31: PetscCall(VecDuplicate(x, &u));
32: PetscCall(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: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
37: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
38: PetscCall(MatSetFromOptions(A));
39: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
66: }
67: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
68: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
70: /* Set exact solution */
71: PetscCall(VecSet(u, one));
73: /* Create linear solver context */
74: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
75: PetscCall(KSPSetOperators(ksp, A, A));
76: PetscCall(KSPGetPC(ksp, &pc));
77: PetscCall(PCSetType(pc, PCLU));
78: #if defined(PETSC_HAVE_MUMPS)
79: if (size > 1) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
80: #endif
81: PetscCall(KSPSetFromOptions(ksp));
83: /* 1. Solve linear system A x = b */
84: PetscCall(MatMult(A, u, b));
85: PetscCall(KSPSolve(ksp, b, x));
87: /* Check the error */
88: PetscCall(VecAXPY(x, -1.0, u));
89: PetscCall(VecNorm(x, NORM_2, &norm));
90: PetscCall(KSPGetIterationNumber(ksp, &its));
91: if (norm > tol) PetscCall(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: PetscCall(MatMultTranspose(A, u, b));
95: PetscCall(KSPSolveTranspose(ksp, b, x2));
97: /* Check the error */
98: PetscCall(VecAXPY(x2, -1.0, u));
99: PetscCall(VecNorm(x2, NORM_2, &norm));
100: PetscCall(KSPGetIterationNumber(ksp, &its));
101: if (norm > tol) PetscCall(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: PetscCall(MatSetValues(A, 1, &i, 1, col, value, ADD_VALUES));
109: }
110: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
111: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
113: PetscCall(MatMult(A, u, b));
114: PetscCall(KSPSolve(ksp, b, x));
116: /* Check the error */
117: PetscCall(VecAXPY(x, -1.0, u));
118: PetscCall(VecNorm(x, NORM_2, &norm));
119: PetscCall(KSPGetIterationNumber(ksp, &its));
120: if (norm > tol) PetscCall(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: PetscCall(VecDestroy(&x));
124: PetscCall(VecDestroy(&u));
125: PetscCall(VecDestroy(&x2));
126: PetscCall(VecDestroy(&b));
127: PetscCall(MatDestroy(&A));
128: PetscCall(KSPDestroy(&ksp));
130: PetscCall(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*/