Actual source code: ex1f.F
petsc-3.7.5 2017-01-01
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: implicit none
13: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
14: ! Include files
15: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
16: !
17: ! This program uses CPP for preprocessing, as indicated by the use of
18: ! PETSc include files in the directory petsc/include/petsc/finclude. This
19: ! convention enables use of the CPP preprocessor, which allows the use
20: ! of the #include statements that define PETSc objects and variables.
21: !
22: ! Use of the conventional Fortran include statements is also supported
23: ! In this case, the PETsc include files are located in the directory
24: ! petsc/include/foldinclude.
25: !
26: ! Since one must be very careful to include each file no more than once
27: ! in a Fortran routine, application programmers must exlicitly list
28: ! each file needed for the various PETSc components within their
29: ! program (unlike the C/C++ interface).
30: !
31: ! See the Fortran section of the PETSc users manual for details.
32: !
33: ! The following include statements are required for KSP Fortran programs:
34: ! petscsys.h - base PETSc routines
35: ! petscvec.h - vectors
36: ! petscmat.h - matrices
37: ! petscksp.h - Krylov subspace methods
38: ! petscpc.h - preconditioners
39: ! Other include statements may be needed if using additional PETSc
40: ! routines in a Fortran program, e.g.,
41: ! petscviewer.h - viewers
42: ! petscis.h - index sets
43: !
44: #include <petsc/finclude/petscsys.h>
45: #include <petsc/finclude/petscvec.h>
46: #include <petsc/finclude/petscmat.h>
47: #include <petsc/finclude/petscksp.h>
48: #include <petsc/finclude/petscpc.h>
49: !
50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: ! Variable declarations
52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: !
54: ! Variables:
55: ! ksp - linear solver context
56: ! ksp - Krylov subspace method context
57: ! pc - preconditioner context
58: ! x, b, u - approx solution, right-hand-side, exact solution vectors
59: ! A - matrix that defines linear system
60: ! its - iterations for convergence
61: ! norm - norm of error in solution
62: !
63: Vec x,b,u
64: Mat A
65: KSP ksp
66: PC pc
67: PetscReal norm,tol
68: PetscErrorCode ierr
69: PetscInt i,n,col(3),its,i1,i2,i3
70: PetscBool flg
71: PetscMPIInt size,rank
72: PetscScalar none,one,value(3)
74: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: ! Beginning of program
76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
79: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
80: if (size .ne. 1) then
81: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
82: if (rank .eq. 0) then
83: write(6,*) 'This is a uniprocessor example only!'
84: endif
85: SETERRQ(PETSC_COMM_WORLD,1,' ',ierr)
86: endif
87: none = -1.0
88: one = 1.0
89: n = 10
90: i1 = 1
91: i2 = 2
92: i3 = 3
93: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
94: & '-n',n,flg,ierr)
96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: ! Compute the matrix and right-hand-side vector that define
98: ! the linear system, Ax = b.
99: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: ! Create matrix. When using MatCreate(), the matrix format can
102: ! be specified at runtime.
104: call MatCreate(PETSC_COMM_WORLD,A,ierr)
105: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
106: call MatSetFromOptions(A,ierr)
107: call MatSetUp(A,ierr)
109: ! Assemble matrix.
110: ! - Note that MatSetValues() uses 0-based row and column numbers
111: ! in Fortran as well as in C (as set here in the array "col").
113: value(1) = -1.0
114: value(2) = 2.0
115: value(3) = -1.0
116: do 50 i=1,n-2
117: col(1) = i-1
118: col(2) = i
119: col(3) = i+1
120: call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr)
121: 50 continue
122: i = n - 1
123: col(1) = n - 2
124: col(2) = n - 1
125: call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
126: i = 0
127: col(1) = 0
128: col(2) = 1
129: value(1) = 2.0
130: value(2) = -1.0
131: call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
132: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
133: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
135: ! Create vectors. Note that we form 1 vector from scratch and
136: ! then duplicate as needed.
138: call VecCreate(PETSC_COMM_WORLD,x,ierr)
139: call VecSetSizes(x,PETSC_DECIDE,n,ierr)
140: call VecSetFromOptions(x,ierr)
141: call VecDuplicate(x,b,ierr)
142: call VecDuplicate(x,u,ierr)
144: ! Set exact solution; then compute right-hand-side vector.
146: call VecSet(u,one,ierr)
147: call MatMult(A,u,b,ierr)
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: ! Create the linear solver and set various options
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: ! Create linear solver context
155: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
157: ! Set operators. Here the matrix that defines the linear system
158: ! also serves as the preconditioning matrix.
160: call KSPSetOperators(ksp,A,A,ierr)
162: ! Set linear solver defaults for this problem (optional).
163: ! - By extracting the KSP and PC contexts from the KSP context,
164: ! we can then directly directly call any KSP and PC routines
165: ! to set various options.
166: ! - The following four statements are optional; all of these
167: ! parameters could alternatively be specified at runtime via
168: ! KSPSetFromOptions();
170: call KSPGetPC(ksp,pc,ierr)
171: call PCSetType(pc,PCJACOBI,ierr)
172: tol = .0000001
173: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, &
174: & PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
176: ! Set runtime options, e.g.,
177: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
178: ! These options will override those specified above as long as
179: ! KSPSetFromOptions() is called _after_ any other customization
180: ! routines.
182: call KSPSetFromOptions(ksp,ierr)
184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: ! Solve the linear system
186: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: call KSPSolve(ksp,b,x,ierr)
190: ! View solver info; we could instead use the option -ksp_view
192: call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)
194: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: ! Check solution and clean up
196: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: ! Check the error
200: call VecAXPY(x,none,u,ierr)
201: call VecNorm(x,NORM_2,norm,ierr)
202: call KSPGetIterationNumber(ksp,its,ierr)
203: if (norm .gt. 1.e-12) then
204: write(6,100) norm,its
205: else
206: write(6,200) its
207: endif
208: 100 format('Norm of error ',e11.4,', Iterations = ',i5)
209: 200 format('Norm of error < 1.e-12, Iterations = ',i5)
211: ! Free work space. All PETSc objects should be destroyed when they
212: ! are no longer needed.
214: call VecDestroy(x,ierr)
215: call VecDestroy(u,ierr)
216: call VecDestroy(b,ierr)
217: call MatDestroy(A,ierr)
218: call KSPDestroy(ksp,ierr)
219: call PetscFinalize(ierr)
221: end