Actual source code: ex1f.F90
1: !
2: ! Description: Solves a tridiagonal linear system with KSP.
3: !
4: ! -----------------------------------------------------------------------
6: program main
7: #include <petsc/finclude/petscksp.h>
8: use petscksp
9: implicit none
11: !
12: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
13: ! Variable declarations
14: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
15: !
16: ! Variables:
17: ! ksp - linear solver context
18: ! ksp - Krylov subspace method context
19: ! pc - preconditioner context
20: ! x, b, u - approx solution, right-hand-side, exact solution vectors
21: ! A - matrix that defines linear system
22: ! its - iterations for convergence
23: ! norm - norm of error in solution
24: !
25: Vec x,b,u
26: Mat A
27: KSP ksp
28: PC pc
29: PetscReal norm,tol
30: PetscErrorCode ierr
31: PetscInt i,n,col(3),its,i1,i2,i3
32: PetscBool flg
33: PetscMPIInt size
34: PetscScalar none,one,value(3)
35: PetscLogStage stages(2);
37: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38: ! Beginning of program
39: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: PetscCallA(PetscInitialize(ierr))
42: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
43: if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
44: none = -1.0
45: one = 1.0
46: n = 10
47: i1 = 1
48: i2 = 2
49: i3 = 3
50: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
52: PetscCallA(PetscLogStageRegister("MatVec Assembly",stages(1),ierr))
53: PetscCallA(PetscLogStageRegister("KSP Solve",stages(2),ierr))
54: PetscCallA(PetscLogStagePush(stages(1),ierr))
55: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: ! Compute the matrix and right-hand-side vector that define
57: ! the linear system, Ax = b.
58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: ! Create matrix. When using MatCreate(), the matrix format can
61: ! be specified at runtime.
63: PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
64: PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
65: PetscCallA(MatSetFromOptions(A,ierr))
66: PetscCallA(MatSetUp(A,ierr))
68: ! Assemble matrix.
69: ! - Note that MatSetValues() uses 0-based row and column numbers
70: ! in Fortran as well as in C (as set here in the array "col").
72: value(1) = -1.0
73: value(2) = 2.0
74: value(3) = -1.0
75: do 50 i=1,n-2
76: col(1) = i-1
77: col(2) = i
78: col(3) = i+1
79: PetscCallA(MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr))
80: 50 continue
81: i = n - 1
82: col(1) = n - 2
83: col(2) = n - 1
84: PetscCallA(MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr))
85: i = 0
86: col(1) = 0
87: col(2) = 1
88: value(1) = 2.0
89: value(2) = -1.0
90: PetscCallA(MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr))
91: PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
92: PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
94: ! Create vectors. Note that we form 1 vector from scratch and
95: ! then duplicate as needed.
97: PetscCallA(VecCreate(PETSC_COMM_WORLD,x,ierr))
98: PetscCallA(VecSetSizes(x,PETSC_DECIDE,n,ierr))
99: PetscCallA(VecSetFromOptions(x,ierr))
100: PetscCallA(VecDuplicate(x,b,ierr))
101: PetscCallA(VecDuplicate(x,u,ierr))
103: ! Set exact solution; then compute right-hand-side vector.
105: PetscCallA(VecSet(u,one,ierr))
106: PetscCallA(MatMult(A,u,b,ierr))
107: PetscCallA(PetscLogStagePop(ierr))
108: PetscCallA(PetscLogStagePush(stages(2),ierr))
110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: ! Create the linear solver and set various options
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: ! Create linear solver context
116: PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp,ierr))
118: ! Set operators. Here the matrix that defines the linear system
119: ! also serves as the preconditioning matrix.
121: PetscCallA(KSPSetOperators(ksp,A,A,ierr))
123: ! Set linear solver defaults for this problem (optional).
124: ! - By extracting the KSP and PC contexts from the KSP context,
125: ! we can then directly call any KSP and PC routines
126: ! to set various options.
127: ! - The following four statements are optional; all of these
128: ! parameters could alternatively be specified at runtime via
129: ! KSPSetFromOptions();
131: PetscCallA(KSPGetPC(ksp,pc,ierr))
132: PetscCallA(PCSetType(pc,PCJACOBI,ierr))
133: tol = .0000001
134: PetscCallA(KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr))
136: ! Set runtime options, e.g.,
137: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
138: ! These options will override those specified above as long as
139: ! KSPSetFromOptions() is called _after_ any other customization
140: ! routines.
142: PetscCallA(KSPSetFromOptions(ksp,ierr))
144: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: ! Solve the linear system
146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: PetscCallA(KSPSolve(ksp,b,x,ierr))
149: PetscCallA(PetscLogStagePop(ierr))
151: ! View solver converged reason; we could instead use the option -ksp_converged_reason
152: PetscCallA(KSPConvergedReasonView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr))
154: ! View solver info; we could instead use the option -ksp_view
156: PetscCallA(KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr))
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: ! Check solution and clean up
160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: ! Check the error
164: PetscCallA(VecAXPY(x,none,u,ierr))
165: PetscCallA(VecNorm(x,NORM_2,norm,ierr))
166: PetscCallA(KSPGetIterationNumber(ksp,its,ierr))
167: if (norm .gt. 1.e-12) then
168: write(6,100) norm,its
169: else
170: write(6,200) its
171: endif
172: 100 format('Norm of error ',e11.4,', Iterations = ',i5)
173: 200 format('Norm of error < 1.e-12, Iterations = ',i5)
175: ! Free work space. All PETSc objects should be destroyed when they
176: ! are no longer needed.
178: PetscCallA(VecDestroy(x,ierr))
179: PetscCallA(VecDestroy(u,ierr))
180: PetscCallA(VecDestroy(b,ierr))
181: PetscCallA(MatDestroy(A,ierr))
182: PetscCallA(KSPDestroy(ksp,ierr))
183: PetscCallA(PetscFinalize(ierr))
185: end
187: !/*TEST
188: !
189: ! test:
190: ! args: -ksp_monitor_short
191: !
192: !TEST*/