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