Actual source code: ex52f.F90

petsc-3.13.4 2020-08-01
Report Typos and Errors
  1: !
  2: !   Modified from ex15f.F for testing MUMPS
  3: !   Solves a linear system in parallel with KSP.  
  4: !  -------------------------------------------------------------------------

  6:       program main
  7:  #include <petsc/finclude/petscksp.h>
  8:       use petscksp
  9:       implicit none


 12: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 13: !                   Variable declarations
 14: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 15:       Vec              x,b,u
 16:       Mat              A
 17:       KSP              ksp
 18:       PetscScalar      v,one,neg_one
 19:       PetscReal        norm,tol
 20:       PetscErrorCode   ierr
 21:       PetscInt         i,j,II,JJ,Istart
 22:       PetscInt         Iend,m,n,i1,its,five
 23:       PetscMPIInt      rank
 24:       PetscBool        flg
 25: #if defined(PETSC_HAVE_MUMPS)
 26:       PC               pc
 27:       Mat              F
 28:       PetscInt         ival,icntl,infog34
 29:       PetscReal        cntl,rinfo12,rinfo13,val
 30: #endif

 32: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 33: !                 Beginning of program
 34: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 35:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 36:       if (ierr .ne. 0) then
 37:         print*,'Unable to initialize PETSc'
 38:         stop
 39:       endif
 40:       one     = 1.0
 41:       neg_one = -1.0
 42:       i1      = 1
 43:       m       = 8
 44:       n       = 7
 45:       five    = 5
 46:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
 47:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
 48:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

 50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 51: !      Compute the matrix and right-hand-side vector that define
 52: !      the linear system, Ax = b.
 53: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 54:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
 55:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
 56:       call MatSetType(A, MATAIJ,ierr)
 57:       call MatSetFromOptions(A,ierr)
 58:       call MatMPIAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,five,PETSC_NULL_INTEGER,ierr)
 59:       call MatSeqAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,ierr)

 61:       call MatGetOwnershipRange(A,Istart,Iend,ierr)

 63: !  Set matrix elements for the 2-D, five-point stencil in parallel.
 64: !   - Each processor needs to insert only elements that it owns
 65: !     locally (but any non-local elements will be sent to the
 66: !     appropriate processor during matrix assembly).
 67: !   - Always specify global row and columns of matrix entries.
 68: !   - Note that MatSetValues() uses 0-based row and column numbers
 69: !     in Fortran as well as in C.
 70:       do 10, II=Istart,Iend-1
 71:         v = -1.0
 72:         i = II/n
 73:         j = II - i*n
 74:         if (i.gt.0) then
 75:           JJ = II - n
 76:           call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
 77:         endif
 78:         if (i.lt.m-1) then
 79:           JJ = II + n
 80:           call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
 81:         endif
 82:         if (j.gt.0) then
 83:           JJ = II - 1
 84:           call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
 85:         endif
 86:         if (j.lt.n-1) then
 87:           JJ = II + 1
 88:           call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
 89:         endif
 90:         v = 4.0
 91:         call  MatSetValues(A,i1,II,i1,II,v,ADD_VALUES,ierr)
 92:  10   continue

 94: !  Assemble matrix, using the 2-step process:
 95:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
 96:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

 98: !  Create parallel vectors.
 99:       call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,u,ierr)
100:       call VecDuplicate(u,b,ierr)
101:       call VecDuplicate(b,x,ierr)

103: !  Set exact solution; then compute right-hand-side vector.
104:       call VecSet(u,one,ierr)
105:       call MatMult(A,u,b,ierr)

107: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: !         Create the linear solver and set various options
109: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110:       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
111:       call KSPSetOperators(ksp,A,A,ierr)
112:       tol = 1.e-7
113:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)

115: !  Test MUMPS
116: #if defined(PETSC_HAVE_MUMPS)
117:       call KSPSetType(ksp,KSPPREONLY,ierr)
118:       call KSPGetPC(ksp,pc,ierr)
119:       call PCSetType(pc,PCLU,ierr)
120:       call PCFactorSetMatSolverType(pc,MATSOLVERMUMPS,ierr)
121:       call PCFactorSetUpMatSolverType(pc,ierr)
122:       call PCFactorGetMatrix(pc,F,ierr)

124: !     sequential ordering
125:       icntl = 7 
126:       ival  = 2
127:       call MatMumpsSetIcntl(F,icntl,ival,ierr)

129: !     threshold for row pivot detection
130:       icntl = 24
131:       ival  = 1
132:       call MatMumpsSetIcntl(F,icntl,ival,ierr)
133:       icntl = 3
134:       val = 1.e-6
135:       call MatMumpsSetCntl(F,icntl,val,ierr)

137: !     compute determinant of A
138:       icntl = 33
139:       ival  = 1
140:       call MatMumpsSetIcntl(F,icntl,ival,ierr)
141: #endif

143:       call KSPSetFromOptions(ksp,ierr)
144:       call KSPSetUp(ksp,ierr)
145: #if defined(PETSC_HAVE_MUMPS)
146:       icntl = 3;
147:       call MatMumpsGetCntl(F,icntl,cntl,ierr)
148:       icntl = 34
149:       call MatMumpsGetInfog(F,icntl,infog34,ierr)
150:       icntl = 12
151:       call MatMumpsGetRinfog(F,icntl,rinfo12,ierr)
152:       icntl = 13
153:       call MatMumpsGetRinfog(F,icntl,rinfo13,ierr)
154:       if (rank .eq. 0) then
155:          write(6,98) cntl
156:          write(6,99) rinfo12,rinfo13,infog34
157:       endif
158:  98   format('Mumps row pivot threshold = ',1pe11.2)
159:  99   format('Mumps determinant=(',1pe11.2,1pe11.2,')*2^',i3)
160: #endif
161:       
162: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: !                      Solve the linear system
164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165:       call KSPSolve(ksp,b,x,ierr)

167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: !                     Check solution and clean up
169: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170:       call VecAXPY(x,neg_one,u,ierr)
171:       call VecNorm(x,NORM_2,norm,ierr)
172:       call KSPGetIterationNumber(ksp,its,ierr)

174:       if (rank .eq. 0) then
175:         if (norm .gt. 1.e-12) then
176:            write(6,100) norm,its
177:         else
178:            write(6,110) its
179:         endif
180:       endif     
181:   100 format('Norm of error ',1pe11.4,' iterations ',i5)
182:   110 format('Norm of error < 1.e-12,iterations ',i5)

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

187:       call KSPDestroy(ksp,ierr)
188:       call VecDestroy(u,ierr)
189:       call VecDestroy(x,ierr)
190:       call VecDestroy(b,ierr)
191:       call MatDestroy(A,ierr)

193: !  Always call PetscFinalize() before exiting a program.
194:       call PetscFinalize(ierr)
195:       end

197: !
198: !/*TEST
199: !
200: !   test:
201: !      suffix: mumps
202: !      nsize: 3
203: !      requires: mumps double
204: !      output_file: output/ex52f_1.out
205: !
206: !TEST*/