Actual source code: ex52f.F

petsc-dev 2014-07-25
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:       implicit none

  9: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 10: !                    Include files
 11: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 12: #include <finclude/petscsys.h>
 13: #include <finclude/petscvec.h>
 14: #include <finclude/petscmat.h>
 15: #include <finclude/petscpc.h>
 16: #include <finclude/petscksp.h>

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

 38: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 39: !                 Beginning of program
 40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 41:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 42:       one     = 1.0
 43:       neg_one = -1.0
 44:       i1      = 1
 45:       m       = 8
 46:       n       = 7
 47:       five    = 5
 48:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
 49:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
 50:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

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

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

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

 97: !  Assemble matrix, using the 2-step process:
 98:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
 99:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

101: !  Create parallel vectors.
102:       call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,u,ierr)
103:       call VecDuplicate(u,b,ierr)
104:       call VecDuplicate(b,x,ierr)

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

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

119: !  Test MUMPS
120: #if defined(PETSC_HAVE_MUMPS)
121:       call KSPSetType(ksp,KSPPREONLY,ierr)
122:       call KSPGetPC(ksp,pc,ierr)
123:       call PCSetType(pc,PCLU,ierr)
124:       call PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS,ierr)
125:       call PCFactorSetUpMatSolverPackage(pc,ierr)
126:       call PCFactorGetMatrix(pc,F,ierr)

128: !     sequential ordering
129:       icntl = 7
130:       ival  = 2
131:       call MatMumpsSetIcntl(F,icntl,ival,ierr)

133: !     threshhold for row pivot detection
134:       call MatMumpsSetIcntl(F,24,1,ierr)
135:       icntl = 3
136:       val = 1.e-6
137:       call MatMumpsSetCntl(F,icntl,val,ierr)

139: !     compute determinant of A
140:       call MatMumpsSetIcntl(F,33,1,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:       call MatMumpsGetInfog(F,34,infog34,ierr)
149:       call MatMumpsGetRinfog(F,12,rinfo12,ierr)
150:       call MatMumpsGetRinfog(F,13,rinfo13,ierr)
151:       if (rank .eq. 0) then
152:          write(6,98) cntl
153:          write(6,99) rinfo12,rinfo13,infog34
154:       endif
155:  98   format('Mumps row pivot threshhold = ',1pe11.2)
156:  99   format('Mumps determinant=(',1pe11.2,1pe11.2,')*2^',i3)
157: #endif
158: 
159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: !                      Solve the linear system
161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162:       call KSPSolve(ksp,b,x,ierr)

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

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

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

184:       call KSPDestroy(ksp,ierr)
185:       call VecDestroy(u,ierr)
186:       call VecDestroy(x,ierr)
187:       call VecDestroy(b,ierr)
188:       call MatDestroy(A,ierr)

190: !  Always call PetscFinalize() before exiting a program.
191:       call PetscFinalize(ierr)
192:       end