Actual source code: ex52f.F90

  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

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

 31: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 32: !                 Beginning of program
 33: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 34:       PetscCallA(PetscInitialize(ierr))
 35:       one     = 1.0
 36:       neg_one = -1.0
 37:       i1      = 1
 38:       m       = 8
 39:       n       = 7
 40:       five    = 5
 41:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr))
 42:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
 43:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))

 45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 46: !      Compute the matrix and right-hand-side vector that define
 47: !      the linear system, Ax = b.
 48: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 49:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
 50:       PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr))
 51:       PetscCallA(MatSetType(A, MATAIJ,ierr))
 52:       PetscCallA(MatSetFromOptions(A,ierr))
 53:       PetscCallA(MatMPIAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,five,PETSC_NULL_INTEGER,ierr))
 54:       PetscCallA(MatSeqAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,ierr))

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

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

 89: !  Assemble matrix, using the 2-step process:
 90:       PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
 91:       PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))

 93: !  Create parallel vectors.
 94:       PetscCallA(VecCreateFromOptions(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,i1,PETSC_DECIDE,m*n,u,ierr))
 95:       PetscCallA(VecDuplicate(u,b,ierr))
 96:       PetscCallA(VecDuplicate(b,x,ierr))

 98: !  Set exact solution; then compute right-hand-side vector.
 99:       PetscCallA(VecSet(u,one,ierr))
100:       PetscCallA(MatMult(A,u,b,ierr))

102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: !         Create the linear solver and set various options
104: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105:       PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp,ierr))
106:       PetscCallA(KSPSetOperators(ksp,A,A,ierr))
107:       tol = 1.e-7
108:       PetscCallA(KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr))

110: !  Test MUMPS
111: #if defined(PETSC_HAVE_MUMPS)
112:       PetscCallA(KSPSetType(ksp,KSPPREONLY,ierr))
113:       PetscCallA(KSPGetPC(ksp,pc,ierr))
114:       PetscCallA(PCSetType(pc,PCLU,ierr))
115:       PetscCallA(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS,ierr))
116:       PetscCallA(PCFactorSetUpMatSolverType(pc,ierr))
117:       PetscCallA(PCFactorGetMatrix(pc,F,ierr))

119: !     sequential ordering
120:       icntl = 7
121:       ival  = 2
122:       PetscCallA(MatMumpsSetIcntl(F,icntl,ival,ierr))

124: !     threshold for row pivot detection
125:       icntl = 24
126:       ival  = 1
127:       PetscCallA(MatMumpsSetIcntl(F,icntl,ival,ierr))
128:       icntl = 3
129:       val = 1.e-6
130:       PetscCallA(MatMumpsSetCntl(F,icntl,val,ierr))

132: !     compute determinant of A
133:       icntl = 33
134:       ival  = 1
135:       PetscCallA(MatMumpsSetIcntl(F,icntl,ival,ierr))
136: #endif

138:       PetscCallA(KSPSetFromOptions(ksp,ierr))
139:       PetscCallA(KSPSetUp(ksp,ierr))
140: #if defined(PETSC_HAVE_MUMPS)
141:       icntl = 3;
142:       PetscCallA(MatMumpsGetCntl(F,icntl,cntl,ierr))
143:       icntl = 34
144:       PetscCallA(MatMumpsGetInfog(F,icntl,infog34,ierr))
145:       icntl = 12
146:       PetscCallA(MatMumpsGetRinfog(F,icntl,rinfo12,ierr))
147:       icntl = 13
148:       PetscCallA(MatMumpsGetRinfog(F,icntl,rinfo13,ierr))
149:       if (rank .eq. 0) then
150:          write(6,98) cntl
151:          write(6,99) rinfo12,rinfo13,infog34
152:       endif
153:  98   format('Mumps row pivot threshold = ',1pe11.2)
154:  99   format('Mumps determinant=(',1pe11.2,1pe11.2,')*2^',i3)
155: #endif

157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: !                      Solve the linear system
159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160:       PetscCallA(KSPSolve(ksp,b,x,ierr))

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

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

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

182:       PetscCallA(KSPDestroy(ksp,ierr))
183:       PetscCallA(VecDestroy(u,ierr))
184:       PetscCallA(VecDestroy(x,ierr))
185:       PetscCallA(VecDestroy(b,ierr))
186:       PetscCallA(MatDestroy(A,ierr))

188: !  Always call PetscFinalize() before exiting a program.
189:       PetscCallA(PetscFinalize(ierr))
190:       end

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