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