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: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
35: if (ierr .ne. 0) then
36: print*,'Unable to initialize PETSc'
37: stop
38: endif
39: one = 1.0
40: neg_one = -1.0
41: i1 = 1
42: m = 8
43: n = 7
44: five = 5
45: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
46: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
47: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: ! Compute the matrix and right-hand-side vector that define
51: ! the linear system, Ax = b.
52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: call MatCreate(PETSC_COMM_WORLD,A,ierr)
54: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
55: call MatSetType(A, MATAIJ,ierr)
56: call MatSetFromOptions(A,ierr)
57: call MatMPIAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,five,PETSC_NULL_INTEGER,ierr)
58: call MatSeqAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,ierr)
60: call MatGetOwnershipRange(A,Istart,Iend,ierr)
62: ! Set matrix elements for the 2-D, five-point stencil in parallel.
63: ! - Each processor needs to insert only elements that it owns
64: ! locally (but any non-local elements will be sent to the
65: ! appropriate processor during matrix assembly).
66: ! - Always specify global row and columns of matrix entries.
67: ! - Note that MatSetValues() uses 0-based row and column numbers
68: ! in Fortran as well as in C.
69: do 10, II=Istart,Iend-1
70: v = -1.0
71: i = II/n
72: j = II - i*n
73: if (i.gt.0) then
74: JJ = II - n
75: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
76: endif
77: if (i.lt.m-1) then
78: JJ = II + n
79: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
80: endif
81: if (j.gt.0) then
82: JJ = II - 1
83: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
84: endif
85: if (j.lt.n-1) then
86: JJ = II + 1
87: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
88: endif
89: v = 4.0
90: call MatSetValues(A,i1,II,i1,II,v,ADD_VALUES,ierr)
91: 10 continue
93: ! Assemble matrix, using the 2-step process:
94: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
95: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
97: ! Create parallel vectors.
98: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,u,ierr)
99: call VecDuplicate(u,b,ierr)
100: call VecDuplicate(b,x,ierr)
102: ! Set exact solution; then compute right-hand-side vector.
103: call VecSet(u,one,ierr)
104: call MatMult(A,u,b,ierr)
106: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: ! Create the linear solver and set various options
108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
110: call KSPSetOperators(ksp,A,A,ierr)
111: tol = 1.e-7
112: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
114: ! Test MUMPS
115: #if defined(PETSC_HAVE_MUMPS)
116: call KSPSetType(ksp,KSPPREONLY,ierr)
117: call KSPGetPC(ksp,pc,ierr)
118: call PCSetType(pc,PCLU,ierr)
119: call PCFactorSetMatSolverType(pc,MATSOLVERMUMPS,ierr)
120: call PCFactorSetUpMatSolverType(pc,ierr)
121: call PCFactorGetMatrix(pc,F,ierr)
123: ! sequential ordering
124: icntl = 7
125: ival = 2
126: call MatMumpsSetIcntl(F,icntl,ival,ierr)
128: ! threshold for row pivot detection
129: icntl = 24
130: ival = 1
131: call MatMumpsSetIcntl(F,icntl,ival,ierr)
132: icntl = 3
133: val = 1.e-6
134: call MatMumpsSetCntl(F,icntl,val,ierr)
136: ! compute determinant of A
137: icntl = 33
138: ival = 1
139: call MatMumpsSetIcntl(F,icntl,ival,ierr)
140: #endif
142: call KSPSetFromOptions(ksp,ierr)
143: call KSPSetUp(ksp,ierr)
144: #if defined(PETSC_HAVE_MUMPS)
145: icntl = 3;
146: call MatMumpsGetCntl(F,icntl,cntl,ierr)
147: icntl = 34
148: call MatMumpsGetInfog(F,icntl,infog34,ierr)
149: icntl = 12
150: call MatMumpsGetRinfog(F,icntl,rinfo12,ierr)
151: icntl = 13
152: call MatMumpsGetRinfog(F,icntl,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 threshold = ',1pe11.2)
158: 99 format('Mumps determinant=(',1pe11.2,1pe11.2,')*2^',i3)
159: #endif
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
196: !
197: !/*TEST
198: !
199: ! test:
200: ! suffix: mumps
201: ! nsize: 3
202: ! requires: mumps double
203: ! output_file: output/ex52f_1.out
204: !
205: !TEST*/