Actual source code: ex52f.F90

petsc-3.10.2 2018-10-09
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: !     threshhold for row pivot detection
130:       call MatMumpsSetIcntl(F,24,1,ierr)
131:       icntl = 3
132:       val = 1.e-6
133:       call MatMumpsSetCntl(F,icntl,val,ierr)

135: !     compute determinant of A
136:       call MatMumpsSetIcntl(F,33,1,ierr)
137: #endif

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

160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: !                     Check solution and clean up
162: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163:       call VecAXPY(x,neg_one,u,ierr)
164:       call VecNorm(x,NORM_2,norm,ierr)
165:       call KSPGetIterationNumber(ksp,its,ierr)

167:       if (rank .eq. 0) then
168:         if (norm .gt. 1.e-12) then
169:            write(6,100) norm,its
170:         else
171:            write(6,110) its
172:         endif
173:       endif
174:   100 format('Norm of error ',1pe11.4,' iterations ',i5)
175:   110 format('Norm of error < 1.e-12,iterations ',i5)

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

180:       call KSPDestroy(ksp,ierr)
181:       call VecDestroy(u,ierr)
182:       call VecDestroy(x,ierr)
183:       call VecDestroy(b,ierr)
184:       call MatDestroy(A,ierr)

186: !  Always call PetscFinalize() before exiting a program.
187:       call PetscFinalize(ierr)
188:       end

190: !
191: !/*TEST
192: !
193: !   test:
194: !      suffix: mumps
195: !      nsize: 3
196: !      requires: mumps
197: !      output_file: output/ex52f_1.out
198: !
199: !TEST*/