Actual source code: ex53.c
2: /* Program usage: ./ex53 [-help] [all PETSc options] */
4: static char help[] = "Solves a tridiagonal linear system with KSP. \n\
5: Modified from ex1.c to illustrate reuse of preconditioner \n\
6: Written as requested by [petsc-maint #63875] \n\n";
8: #include <petscksp.h>
12: int main(int argc,char **args)
13: {
14: Vec x,x2,b,u; /* approx solution, RHS, exact solution */
15: Mat A; /* linear system matrix */
16: KSP ksp; /* linear solver context */
17: PC pc; /* preconditioner context */
18: PetscReal norm; /* norm of solution error */
20: PetscInt i,n = 10,col[3],its;
21: PetscMPIInt rank;
22: PetscScalar neg_one = -1.0,one = 1.0,value[3];
24: PetscInitialize(&argc,&args,(char *)0,help);
25: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
26: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
28: /* Create vectors.*/
29: VecCreate(PETSC_COMM_WORLD,&x);
30: PetscObjectSetName((PetscObject) x, "Solution");
31: VecSetSizes(x,PETSC_DECIDE,n);
32: VecSetFromOptions(x);
33: VecDuplicate(x,&b);
34: VecDuplicate(x,&u);
35: VecDuplicate(x,&x2);
37: /* Create matrix. Only proc[0] sets values - not efficient for parallel processing!
38: See ex23.c for efficient parallel assembly matrix */
39: MatCreate(PETSC_COMM_WORLD,&A);
40: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
41: MatSetFromOptions(A);
43: if (!rank){
44: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
45: for (i=1; i<n-1; i++) {
46: col[0] = i-1; col[1] = i; col[2] = i+1;
47: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
48: }
49: i = n - 1; col[0] = n - 2; col[1] = n - 1;
50: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
51: i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
52: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
54: i = 0; col[0] = n-1; value[0] = 0.5; /* make A non-symmetric */
55: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
56: }
57: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
58: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
60: /* Set exact solution */
61: VecSet(u,one);
63: /* Create linear solver context */
64: KSPCreate(PETSC_COMM_WORLD,&ksp);
65: KSPSetOperators(ksp,A,A,SAME_PRECONDITIONER);
66: KSPGetPC(ksp,&pc);
67: PCSetType(pc,PCLU);
68: #ifdef PETSC_HAVE_MUMPS
69: PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
70: #endif
71: KSPSetFromOptions(ksp);
72:
73: /* 1. Solve linear system A x = b */
74: MatMult(A,u,b);
75: KSPSolve(ksp,b,x);
77: /* Check the error */
78: VecAXPY(x,neg_one,u);
79: VecNorm(x,NORM_2,&norm);
80: KSPGetIterationNumber(ksp,&its);
81: PetscPrintf(PETSC_COMM_WORLD,"1. Norm of error for Ax=b: %A, Iterations %D\n",
82: norm,its);
83:
84: /* 2. Solve linear system A^T x = b*/
85: MatMultTranspose(A,u,b);
86: KSPSolveTranspose(ksp,b,x2);
87:
88: /* Check the error */
89: VecAXPY(x2,neg_one,u);
90: VecNorm(x2,NORM_2,&norm);
91: KSPGetIterationNumber(ksp,&its);
92: PetscPrintf(PETSC_COMM_WORLD,"2. Norm of error for A^T x=b: %A, Iterations %D\n",
93: norm,its);
95: /* 3. Change A and solve A x = b with an iterative solver using A=LU as a preconditioner*/
96: if (!rank){
97: i = 0; col[0] = n-1; value[0] = 1.e-2;
98: MatSetValues(A,1,&i,1,col,value,ADD_VALUES);
99: }
100: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
103: MatMult(A,u,b);
104: KSPSolve(ksp,b,x);
106: /* Check the error */
107: VecAXPY(x,neg_one,u);
108: VecNorm(x,NORM_2,&norm);
109: KSPGetIterationNumber(ksp,&its);
110: PetscPrintf(PETSC_COMM_WORLD,"3. Norm of error for (A+Delta) x=b: %A, Iterations %D\n",
111: norm,its);
113: /* Free work space. */
114: VecDestroy(&x);
115: VecDestroy(&u);
116: VecDestroy(&x2);
117: VecDestroy(&b);
118: MatDestroy(&A);
119: KSPDestroy(&ksp);
121: PetscFinalize();
122: return 0;
123: }