Actual source code: ex1f.F
1: !
2: ! "$Id: ex1f.F,v 1.30 2001/08/07 03:04:00 balay Exp $";
3: !
4: ! Description: Solves a tridiagonal linear system with SLES.
5: !
6: !/*T
7: ! Concepts: SLES^solving a system of linear equations
8: ! Processors: 1
9: !T*/
10: ! -----------------------------------------------------------------------
12: program main
13: implicit none
15: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
16: ! Include files
17: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
18: !
19: ! This program uses CPP for preprocessing, as indicated by the use of
20: ! PETSc include files in the directory petsc/include/finclude. This
21: ! convention enables use of the CPP preprocessor, which allows the use
22: ! of the #include statements that define PETSc objects and variables.
23: !
24: ! Use of the conventional Fortran include statements is also supported
25: ! In this case, the PETsc include files are located in the directory
26: ! petsc/include/foldinclude.
27: !
28: ! Since one must be very careful to include each file no more than once
29: ! in a Fortran routine, application programmers must exlicitly list
30: ! each file needed for the various PETSc components within their
31: ! program (unlike the C/C++ interface).
32: !
33: ! See the Fortran section of the PETSc users manual for details.
34: !
35: ! The following include statements are required for SLES Fortran programs:
36: ! petsc.h - base PETSc routines
37: ! petscvec.h - vectors
38: ! petscmat.h - matrices
39: ! petscksp.h - Krylov subspace methods
40: ! petscpc.h - preconditioners
41: ! petscsles.h - SLES interface
42: ! Other include statements may be needed if using additional PETSc
43: ! routines in a Fortran program, e.g.,
44: ! petscviewer.h - viewers
45: ! petscis.h - index sets
46: !
47: #include include/finclude/petsc.h
48: #include include/finclude/petscvec.h
49: #include include/finclude/petscmat.h
50: #include include/finclude/petscksp.h
51: #include include/finclude/petscpc.h
52: #include include/finclude/petscsles.h
53: !
54: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: ! Variable declarations
56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: !
58: ! Variables:
59: ! sles - linear solver context
60: ! ksp - Krylov subspace method context
61: ! pc - preconditioner context
62: ! x, b, u - approx solution, right-hand-side, exact solution vectors
63: ! A - matrix that defines linear system
64: ! its - iterations for convergence
65: ! norm - norm of error in solution
66: !
67: Vec x,b,u
68: Mat A
69: KSP ksp
70: PC pc
71: SLES sles
72: double precision norm,tol
73: integer ierr,i,n,col(3),its,flg,size,rank
74: PetscScalar none,one,value(3)
76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: ! Beginning of program
78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
81: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
82: if (size .ne. 1) then
83: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
84: if (rank .eq. 0) then
85: write(6,*) 'This is a uniprocessor example only!'
86: endif
87: SETERRQ(1,' ',ierr)
88: endif
89: none = -1.0
90: one = 1.0
91: n = 10
92: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: ! Compute the matrix and right-hand-side vector that define
96: ! the linear system, Ax = b.
97: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: ! Create matrix. When using MatCreate(), the matrix format can
100: ! be specified at runtime.
102: call MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,A, &
103: & ierr)
104: call MatSetFromOptions(A,ierr)
106: ! Assemble matrix.
107: ! - Note that MatSetValues() uses 0-based row and column numbers
108: ! in Fortran as well as in C (as set here in the array "col").
110: value(1) = -1.0
111: value(2) = 2.0
112: value(3) = -1.0
113: do 50 i=1,n-2
114: col(1) = i-1
115: col(2) = i
116: col(3) = i+1
117: call MatSetValues(A,1,i,3,col,value,INSERT_VALUES,ierr)
118: 50 continue
119: i = n - 1
120: col(1) = n - 2
121: col(2) = n - 1
122: call MatSetValues(A,1,i,2,col,value,INSERT_VALUES,ierr)
123: i = 0
124: col(1) = 0
125: col(2) = 1
126: value(1) = 2.0
127: value(2) = -1.0
128: call MatSetValues(A,1,i,2,col,value,INSERT_VALUES,ierr)
129: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
130: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
132: ! Create vectors. Note that we form 1 vector from scratch and
133: ! then duplicate as needed.
135: call VecCreate(PETSC_COMM_WORLD,x,ierr)
136: call VecSetSizes(x,PETSC_DECIDE,n,ierr)
137: call VecSetFromOptions(x,ierr)
138: call VecDuplicate(x,b,ierr)
139: call VecDuplicate(x,u,ierr)
141: ! Set exact solution; then compute right-hand-side vector.
143: call VecSet(one,u,ierr)
144: call MatMult(A,u,b,ierr)
146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: ! Create the linear solver and set various options
148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: ! Create linear solver context
152: call SLESCreate(PETSC_COMM_WORLD,sles,ierr)
154: ! Set operators. Here the matrix that defines the linear system
155: ! also serves as the preconditioning matrix.
157: call SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN,ierr)
159: ! Set linear solver defaults for this problem (optional).
160: ! - By extracting the KSP and PC contexts from the SLES context,
161: ! we can then directly directly call any KSP and PC routines
162: ! to set various options.
163: ! - The following four statements are optional; all of these
164: ! parameters could alternatively be specified at runtime via
165: ! SLESSetFromOptions();
167: call SLESGetKSP(sles,ksp,ierr)
168: call SLESGetPC(sles,pc,ierr)
169: call PCSetType(pc,PCJACOBI,ierr)
170: tol = 1.d-7
171: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION, &
172: & PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr)
174: ! Set runtime options, e.g.,
175: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
176: ! These options will override those specified above as long as
177: ! SLESSetFromOptions() is called _after_ any other customization
178: ! routines.
180: call SLESSetFromOptions(sles,ierr)
182: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: ! Solve the linear system
184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: call SLESSolve(sles,b,x,its,ierr)
188: ! View solver info; we could instead use the option -sles_view
190: call SLESView(sles,PETSC_VIEWER_STDOUT_WORLD,ierr)
192: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193: ! Check solution and clean up
194: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: ! Check the error
198: call VecAXPY(none,u,x,ierr)
199: call VecNorm(x,NORM_2,norm,ierr)
200: if (norm .gt. 1.e-12) then
201: write(6,100) norm,its
202: else
203: write(6,200) its
204: endif
205: 100 format('Norm of error = ',e10.4,', Iterations = ',i5)
206: 200 format('Norm of error < 1.e-12,Iterations = ',i5)
208: ! Free work space. All PETSc objects should be destroyed when they
209: ! are no longer needed.
211: call VecDestroy(x,ierr)
212: call VecDestroy(u,ierr)
213: call VecDestroy(b,ierr)
214: call MatDestroy(A,ierr)
215: call SLESDestroy(sles,ierr)
216: call PetscFinalize(ierr)
218: end