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: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
42: if (ierr .ne. 0) then
43: print*,'Unable to initialize PETSc'
44: stop
45: endif
46: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
47: if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
48: none = -1.0
49: one = 1.0
50: n = 10
51: i1 = 1
52: i2 = 2
53: i3 = 3
54: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
56: call PetscLogStageRegister("MatVec Assembly",stages(1),ierr)
57: call PetscLogStageRegister("KSP Solve",stages(2),ierr)
58: call PetscLogStagePush(stages(1),ierr)
59: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: ! Compute the matrix and right-hand-side vector that define
61: ! the linear system, Ax = b.
62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: ! Create matrix. When using MatCreate(), the matrix format can
65: ! be specified at runtime.
67: call MatCreate(PETSC_COMM_WORLD,A,ierr)
68: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
69: call MatSetFromOptions(A,ierr)
70: call MatSetUp(A,ierr)
72: ! Assemble matrix.
73: ! - Note that MatSetValues() uses 0-based row and column numbers
74: ! in Fortran as well as in C (as set here in the array "col").
76: value(1) = -1.0
77: value(2) = 2.0
78: value(3) = -1.0
79: do 50 i=1,n-2
80: col(1) = i-1
81: col(2) = i
82: col(3) = i+1
83: call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr)
84: 50 continue
85: i = n - 1
86: col(1) = n - 2
87: col(2) = n - 1
88: call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
89: i = 0
90: col(1) = 0
91: col(2) = 1
92: value(1) = 2.0
93: value(2) = -1.0
94: call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
95: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
96: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
98: ! Create vectors. Note that we form 1 vector from scratch and
99: ! then duplicate as needed.
101: call VecCreate(PETSC_COMM_WORLD,x,ierr)
102: call VecSetSizes(x,PETSC_DECIDE,n,ierr)
103: call VecSetFromOptions(x,ierr)
104: call VecDuplicate(x,b,ierr)
105: call VecDuplicate(x,u,ierr)
107: ! Set exact solution; then compute right-hand-side vector.
109: call VecSet(u,one,ierr)
110: call MatMult(A,u,b,ierr)
111: call PetscLogStagePop(ierr)
112: call PetscLogStagePush(stages(2),ierr)
114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: ! Create the linear solver and set various options
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: ! Create linear solver context
120: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
122: ! Set operators. Here the matrix that defines the linear system
123: ! also serves as the preconditioning matrix.
125: call KSPSetOperators(ksp,A,A,ierr)
127: ! Set linear solver defaults for this problem (optional).
128: ! - By extracting the KSP and PC contexts from the KSP context,
129: ! we can then directly directly call any KSP and PC routines
130: ! to set various options.
131: ! - The following four statements are optional; all of these
132: ! parameters could alternatively be specified at runtime via
133: ! KSPSetFromOptions();
135: call KSPGetPC(ksp,pc,ierr)
136: call PCSetType(pc,PCJACOBI,ierr)
137: tol = .0000001
138: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, &
139: & PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
141: ! Set runtime options, e.g.,
142: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
143: ! These options will override those specified above as long as
144: ! KSPSetFromOptions() is called _after_ any other customization
145: ! routines.
147: call KSPSetFromOptions(ksp,ierr)
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: ! Solve the linear system
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: call KSPSolve(ksp,b,x,ierr)
154: call PetscLogStagePop(ierr)
156: ! View solver converged reason; we could instead use the option -ksp_converged_reason
157: call KSPConvergedReasonView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)
159: ! View solver info; we could instead use the option -ksp_view
161: call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: ! Check solution and clean up
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: ! Check the error
169: call VecAXPY(x,none,u,ierr)
170: call VecNorm(x,NORM_2,norm,ierr)
171: call KSPGetIterationNumber(ksp,its,ierr)
172: if (norm .gt. 1.e-12) then
173: write(6,100) norm,its
174: else
175: write(6,200) its
176: endif
177: 100 format('Norm of error ',e11.4,', Iterations = ',i5)
178: 200 format('Norm of error < 1.e-12, Iterations = ',i5)
180: ! Free work space. All PETSc objects should be destroyed when they
181: ! are no longer needed.
183: call VecDestroy(x,ierr)
184: call VecDestroy(u,ierr)
185: call VecDestroy(b,ierr)
186: call MatDestroy(A,ierr)
187: call KSPDestroy(ksp,ierr)
188: call PetscFinalize(ierr)
190: end
192: !/*TEST
193: !
194: ! test:
195: ! args: -ksp_monitor_short
196: !
197: !TEST*/