Actual source code: ex52f.F90
1: !
2: ! Modified from ex15f.F for testing MUMPS
3: ! Solves a linear system in parallel with KSP.
4: ! -------------------------------------------------------------------------
5: #include <petsc/finclude/petscksp.h>
6: program main
7: use petscksp
8: implicit none
10: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
11: ! Variable declarations
12: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
13: Vec x, b, u
14: Mat A
15: KSP ksp
16: PetscScalar, parameter :: one = 1.0, neg_one = -1.0
17: PetscScalar v
18: PetscReal norm
19: PetscReal, parameter :: tol = 1.e-7
20: PetscErrorCode ierr
21: PetscInt i, j, II, JJ, Istart, Iend
22: PetscInt m, n, its
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: PetscCallA(PetscInitialize(ierr))
36: m = 8
37: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
38: n = 7
39: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
40: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: ! Compute the matrix and right-hand-side vector that define
44: ! the linear system, Ax = b.
45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
47: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*n, m*n, ierr))
48: PetscCallA(MatSetType(A, MATAIJ, ierr))
49: PetscCallA(MatSetFromOptions(A, ierr))
50: PetscCallA(MatMPIAIJSetPreallocation(A, 5_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 5_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
51: PetscCallA(MatSeqAIJSetPreallocation(A, 5_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
53: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
55: ! Set matrix elements for the 2-D, five-point stencil in parallel.
56: ! - Each processor needs to insert only elements that it owns
57: ! locally (but any non-local elements will be sent to the
58: ! appropriate processor during matrix assembly).
59: ! - Always specify global row and columns of matrix entries.
60: ! - Note that MatSetValues() uses 0-based row and column numbers
61: ! in Fortran as well as in C.
62: do II = Istart, Iend - 1
63: v = -1.0
64: i = II/n
65: j = II - i*n
66: if (i > 0) then
67: JJ = II - n
68: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
69: end if
70: if (i < m - 1) then
71: JJ = II + n
72: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
73: end if
74: if (j > 0) then
75: JJ = II - 1
76: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
77: end if
78: if (j < n - 1) then
79: JJ = II + 1
80: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
81: end if
82: v = 4.0
83: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [II], [v], ADD_VALUES, ierr))
84: end do
86: ! Assemble matrix, using the 2-step process:
87: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
88: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
90: ! Create parallel vectors.
91: PetscCallA(VecCreateFromOptions(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, 1_PETSC_INT_KIND, PETSC_DECIDE, m*n, u, ierr))
92: PetscCallA(VecDuplicate(u, b, ierr))
93: PetscCallA(VecDuplicate(b, x, ierr))
95: ! Set exact solution; then compute right-hand-side vector.
96: PetscCallA(VecSet(u, one, ierr))
97: PetscCallA(MatMult(A, u, b, ierr))
99: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: ! Create the linear solver and set various options
101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
103: PetscCallA(KSPSetOperators(ksp, A, A, ierr))
104: PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))
106: ! Test MUMPS
107: #if defined(PETSC_HAVE_MUMPS)
108: PetscCallA(KSPSetType(ksp, KSPPREONLY, ierr))
109: PetscCallA(KSPGetPC(ksp, pc, ierr))
110: PetscCallA(PCSetType(pc, PCLU, ierr))
111: PetscCallA(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS, ierr))
112: PetscCallA(PCFactorSetUpMatSolverType(pc, ierr))
113: PetscCallA(PCFactorGetMatrix(pc, F, ierr))
115: ! sequential ordering
116: icntl = 7
117: ival = 2
118: PetscCallA(MatMumpsSetIcntl(F, icntl, ival, ierr))
120: ! threshold for row pivot detection
121: icntl = 24
122: ival = 1
123: PetscCallA(MatMumpsSetIcntl(F, icntl, ival, ierr))
124: icntl = 3
125: val = 1.e-6
126: PetscCallA(MatMumpsSetCntl(F, icntl, val, ierr))
128: ! compute determinant of A
129: icntl = 33
130: ival = 1
131: PetscCallA(MatMumpsSetIcntl(F, icntl, ival, ierr))
132: #endif
134: PetscCallA(KSPSetFromOptions(ksp, ierr))
135: PetscCallA(KSPSetUp(ksp, ierr))
136: #if defined(PETSC_HAVE_MUMPS)
137: icntl = 3
138: PetscCallA(MatMumpsGetCntl(F, icntl, cntl, ierr))
139: icntl = 34
140: PetscCallA(MatMumpsGetInfog(F, icntl, infog34, ierr))
141: icntl = 12
142: PetscCallA(MatMumpsGetRinfog(F, icntl, rinfo12, ierr))
143: icntl = 13
144: PetscCallA(MatMumpsGetRinfog(F, icntl, rinfo13, ierr))
145: if (rank == 0) then
146: write (6, 98) cntl
147: write (6, 99) rinfo12, rinfo13, infog34
148: end if
149: 98 format('Mumps row pivot threshold = ', 1pe11.2)
150: 99 format('Mumps determinant=(', 1pe11.2, 1pe11.2, ')*2^', i3)
151: #endif
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: ! Solve the linear system
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: PetscCallA(KSPSolve(ksp, b, x, ierr))
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: ! Check solution and clean up
160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: PetscCallA(VecAXPY(x, neg_one, u, ierr))
162: PetscCallA(VecNorm(x, NORM_2, norm, ierr))
163: PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
165: if (rank == 0) then
166: if (norm > 1.e-12) then
167: write (6, 100) norm, its
168: else
169: write (6, 110) its
170: end if
171: end if
172: 100 format('Norm of error ', 1pe11.4, ' iterations ', i5)
173: 110 format('Norm of error < 1.e-12,iterations ', i5)
175: ! Free work space. All PETSc objects should be destroyed when they
176: ! are no longer needed.
178: PetscCallA(KSPDestroy(ksp, ierr))
179: PetscCallA(VecDestroy(u, ierr))
180: PetscCallA(VecDestroy(x, ierr))
181: PetscCallA(VecDestroy(b, ierr))
182: PetscCallA(MatDestroy(A, ierr))
184: ! Always call PetscFinalize() before exiting a program.
185: PetscCallA(PetscFinalize(ierr))
186: end
188: !
189: !/*TEST
190: !
191: ! test:
192: ! suffix: mumps
193: ! nsize: 3
194: ! requires: mumps double
195: ! output_file: output/ex52f_1.out
196: !
197: !TEST*/