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