Actual source code: ex52f.F

petsc-master 2016-09-24
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 <petsc/finclude/petscsys.h>
 13:  #include <petsc/finclude/petscvec.h>
 14:  #include <petsc/finclude/petscmat.h>
 15:  #include <petsc/finclude/petscpc.h>
 16:  #include <petsc/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_OBJECT,PETSC_NULL_CHARACTER,        &
 49:      &                        '-m',m,flg,ierr)
 50:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,        &
 51:      &                        '-n',n,flg,ierr)
 52:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

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

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

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

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

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

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

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

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

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

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

141: !     compute determinant of A
142:       call MatMumpsSetIcntl(F,33,1,ierr)
143: #endif

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

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

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

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

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

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