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