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
72:         endif
73:         if (i.lt.m-1) then
74:           JJ = II + n
76:         endif
77:         if (j.gt.0) then
78:           JJ = II - 1
80:         endif
81:         if (j.lt.n-1) then
82:           JJ = II + 1
84:         endif
85:         v = 4.0
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*/
```