Actual source code: ex15f.F90

  1: !
  2: !   Solves a linear system in parallel with KSP.  Also indicates
  3: !   use of a user-provided preconditioner.  Input parameters include:
  4: !      -user_defined_pc : Activate a user-defined preconditioner
  5: !
  6: !

  8: !
  9: !     -------------------------------------------------------------------------
 10: !
 11: !     Module contains diag needed by shell preconditioner
 12: !
 13:       module ex15fmodule
 14: #include <petsc/finclude/petscksp.h>
 15:       use petscksp
 16:       Vec    diag
 17:       end module

 19:       program main
 20:       use ex15fmodule
 21:       implicit none

 23: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 24: !                   Variable declarations
 25: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 26: !
 27: !  Variables:
 28: !     ksp     - linear solver context
 29: !     ksp      - Krylov subspace method context
 30: !     pc       - preconditioner context
 31: !     x, b, u  - approx solution, right-hand-side, exact solution vectors
 32: !     A        - matrix that defines linear system
 33: !     its      - iterations for convergence
 34: !     norm     - norm of solution error

 36:       Vec              x,b,u
 37:       Mat              A
 38:       PC               pc
 39:       KSP              ksp
 40:       PetscScalar      v,one,neg_one
 41:       PetscReal norm,tol
 42:       PetscErrorCode ierr
 43:       PetscInt   i,j,II,JJ,Istart
 44:       PetscInt   Iend,m,n,i1,its,five
 45:       PetscMPIInt rank
 46:       PetscBool  user_defined_pc,flg

 48: !  Note: Any user-defined Fortran routines MUST be declared as external.

 50:       external SampleShellPCSetUp, SampleShellPCApply
 51:       external  SampleShellPCDestroy

 53: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 54: !                 Beginning of program
 55: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 57:       PetscCallA(PetscInitialize(ierr))
 58:       one     = 1.0
 59:       neg_one = -1.0
 60:       i1 = 1
 61:       m       = 8
 62:       n       = 7
 63:       five    = 5
 64:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr))
 65:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
 66:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))

 68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69: !      Compute the matrix and right-hand-side vector that define
 70: !      the linear system, Ax = b.
 71: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 73: !  Create parallel matrix, specifying only its global dimensions.
 74: !  When using MatCreate(), the matrix format can be specified at
 75: !  runtime. Also, the parallel partitioning of the matrix is
 76: !  determined by PETSc at runtime.

 78:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
 79:       PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr))
 80:       PetscCallA(MatSetType(A, MATAIJ,ierr))
 81:       PetscCallA(MatSetFromOptions(A,ierr))
 82:       PetscCallA(MatMPIAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,five,PETSC_NULL_INTEGER,ierr))
 83:       PetscCallA(MatSeqAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,ierr))

 85: !  Currently, all PETSc parallel matrix formats are partitioned by
 86: !  contiguous chunks of rows across the processors.  Determine which
 87: !  rows of the matrix are locally owned.

 89:       PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))

 91: !  Set matrix elements for the 2-D, five-point stencil in parallel.
 92: !   - Each processor needs to insert only elements that it owns
 93: !     locally (but any non-local elements will be sent to the
 94: !     appropriate processor during matrix assembly).
 95: !   - Always specify global row and columns of matrix entries.
 96: !   - Note that MatSetValues() uses 0-based row and column numbers
 97: !     in Fortran as well as in C.

 99:       do 10, II=Istart,Iend-1
100:         v = -1.0
101:         i = II/n
102:         j = II - i*n
103:         if (i.gt.0) then
104:           JJ = II - n
105:           PetscCallA(MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr))
106:         endif
107:         if (i.lt.m-1) then
108:           JJ = II + n
109:           PetscCallA(MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr))
110:         endif
111:         if (j.gt.0) then
112:           JJ = II - 1
113:           PetscCallA(MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr))
114:         endif
115:         if (j.lt.n-1) then
116:           JJ = II + 1
117:           PetscCallA(MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr))
118:         endif
119:         v = 4.0
120:         PetscCallA( MatSetValues(A,i1,II,i1,II,v,ADD_VALUES,ierr))
121:  10   continue

123: !  Assemble matrix, using the 2-step process:
124: !       MatAssemblyBegin(), MatAssemblyEnd()
125: !  Computations can be done while messages are in transition,
126: !  by placing code between these two statements.

128:       PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
129:       PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))

131: !  Create parallel vectors.
132: !   - Here, the parallel partitioning of the vector is determined by
133: !     PETSc at runtime.  We could also specify the local dimensions
134: !     if desired -- or use the more general routine VecCreate().
135: !   - When solving a linear system, the vectors and matrices MUST
136: !     be partitioned accordingly.  PETSc automatically generates
137: !     appropriately partitioned matrices and vectors when MatCreate()
138: !     and VecCreate() are used with the same communicator.
139: !   - Note: We form 1 vector from scratch and then duplicate as needed.

141:       PetscCallA(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,u,ierr))
142:       PetscCallA(VecDuplicate(u,b,ierr))
143:       PetscCallA(VecDuplicate(b,x,ierr))

145: !  Set exact solution; then compute right-hand-side vector.

147:       PetscCallA(VecSet(u,one,ierr))
148:       PetscCallA(MatMult(A,u,b,ierr))

150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: !         Create the linear solver and set various options
152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

154: !  Create linear solver context

156:       PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp,ierr))

158: !  Set operators. Here the matrix that defines the linear system
159: !  also serves as the preconditioning matrix.

161:       PetscCallA(KSPSetOperators(ksp,A,A,ierr))

163: !  Set linear solver defaults for this problem (optional).
164: !   - By extracting the KSP and PC contexts from the KSP context,
165: !     we can then directly call any KSP and PC routines
166: !     to set various options.

168:       PetscCallA(KSPGetPC(ksp,pc,ierr))
169:       tol = 1.e-7
170:       PetscCallA(KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr))

172: !
173: !  Set a user-defined shell preconditioner if desired
174: !
175:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-user_defined_pc',user_defined_pc,ierr))

177:       if (user_defined_pc) then

179: !  (Required) Indicate to PETSc that we are using a shell preconditioner
180:          PetscCallA(PCSetType(pc,PCSHELL,ierr))

182: !  (Required) Set the user-defined routine for applying the preconditioner
183:          PetscCallA(PCShellSetApply(pc,SampleShellPCApply,ierr))

185: !  (Optional) Do any setup required for the preconditioner
186:          PetscCallA(PCShellSetSetUp(pc,SampleShellPCSetUp,ierr))

188: !  (Optional) Frees any objects we created for the preconditioner
189:          PetscCallA(PCShellSetDestroy(pc,SampleShellPCDestroy,ierr))

191:       else
192:          PetscCallA(PCSetType(pc,PCJACOBI,ierr))
193:       endif

195: !  Set runtime options, e.g.,
196: !      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
197: !  These options will override those specified above as long as
198: !  KSPSetFromOptions() is called _after_ any other customization
199: !  routines.

201:       PetscCallA(KSPSetFromOptions(ksp,ierr))

203: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: !                      Solve the linear system
205: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

207:       PetscCallA(KSPSolve(ksp,b,x,ierr))

209: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210: !                     Check solution and clean up
211: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

213: !  Check the error

215:       PetscCallA(VecAXPY(x,neg_one,u,ierr))
216:       PetscCallA(VecNorm(x,NORM_2,norm,ierr))
217:       PetscCallA(KSPGetIterationNumber(ksp,its,ierr))

219:       if (rank .eq. 0) then
220:         if (norm .gt. 1.e-12) then
221:            write(6,100) norm,its
222:         else
223:            write(6,110) its
224:         endif
225:       endif
226:   100 format('Norm of error ',1pe11.4,' iterations ',i5)
227:   110 format('Norm of error < 1.e-12,iterations ',i5)

229: !  Free work space.  All PETSc objects should be destroyed when they
230: !  are no longer needed.

232:       PetscCallA(KSPDestroy(ksp,ierr))
233:       PetscCallA(VecDestroy(u,ierr))
234:       PetscCallA(VecDestroy(x,ierr))
235:       PetscCallA(VecDestroy(b,ierr))
236:       PetscCallA(MatDestroy(A,ierr))

238: !  Always call PetscFinalize() before exiting a program.

240:       PetscCallA(PetscFinalize(ierr))
241:       end

243: !/***********************************************************************/
244: !/*          Routines for a user-defined shell preconditioner           */
245: !/***********************************************************************/

247: !
248: !   SampleShellPCSetUp - This routine sets up a user-defined
249: !   preconditioner context.
250: !
251: !   Input Parameters:
252: !   pc - preconditioner object
253: !
254: !   Output Parameter:
255: !   ierr  - error code (nonzero if error has been detected)
256: !
257: !   Notes:
258: !   In this example, we define the shell preconditioner to be Jacobi
259: !   method.  Thus, here we create a work vector for storing the reciprocal
260: !   of the diagonal of the preconditioner matrix; this vector is then
261: !   used within the routine SampleShellPCApply().
262: !
263:       subroutine SampleShellPCSetUp(pc,ierr)
264:       use ex15fmodule
265:       use petscksp
266:       implicit none

268:       PC      pc
269:       Mat     pmat
270:       PetscErrorCode ierr

272:       PetscCallA(PCGetOperators(pc,PETSC_NULL_MAT,pmat,ierr))
273:       PetscCallA(MatCreateVecs(pmat,diag,PETSC_NULL_VEC,ierr))
274:       PetscCallA(MatGetDiagonal(pmat,diag,ierr))
275:       PetscCallA(VecReciprocal(diag,ierr))

277:       end

279: ! -------------------------------------------------------------------
280: !
281: !   SampleShellPCApply - This routine demonstrates the use of a
282: !   user-provided preconditioner.
283: !
284: !   Input Parameters:
285: !   pc - preconditioner object
286: !   x - input vector
287: !
288: !   Output Parameters:
289: !   y - preconditioned vector
290: !   ierr  - error code (nonzero if error has been detected)
291: !
292: !   Notes:
293: !   This code implements the Jacobi preconditioner, merely as an
294: !   example of working with a PCSHELL.  Note that the Jacobi method
295: !   is already provided within PETSc.
296: !
297:       subroutine SampleShellPCApply(pc,x,y,ierr)
298:       use ex15fmodule
299:       implicit none

301:       PC      pc
302:       Vec     x,y
303:       PetscErrorCode ierr

305:       PetscCallA(VecPointwiseMult(y,x,diag,ierr))

307:       end

309: !/***********************************************************************/
310: !/*          Routines for a user-defined shell preconditioner           */
311: !/***********************************************************************/

313: !
314: !   SampleShellPCDestroy - This routine destroys (frees the memory of) any
315: !      objects we made for the preconditioner
316: !
317: !   Input Parameters:
318: !   pc - for this example we use the actual PC as our shell context
319: !
320: !   Output Parameter:
321: !   ierr  - error code (nonzero if error has been detected)
322: !

324:       subroutine SampleShellPCDestroy(pc,ierr)
325:       use ex15fmodule
326:       implicit none

328:       PC      pc
329:       PetscErrorCode ierr

331: !  Normally we would recommend storing all the work data (like diag) in
332: !  the context set with PCShellSetContext()

334:       PetscCallA(VecDestroy(diag,ierr))

336:       end

338: !
339: !/*TEST
340: !
341: !   test:
342: !      nsize: 2
343: !      args: -ksp_view -user_defined_pc -ksp_gmres_cgs_refinement_type refine_always
344: !
345: !TEST*/