Actual source code: ex52f.F
petsc-3.7.5 2017-01-01
1: !
2: ! Modified from ex15f.F for testing MUMPS
3: ! Solves a linear system in parallel with KSP.
4: ! -------------------------------------------------------------------------
6: program main
7: implicit none
9: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
10: ! Include files
11: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
12: #include <petsc/finclude/petscsys.h>
13: #include <petsc/finclude/petscvec.h>
14: #include <petsc/finclude/petscmat.h>
15: #include <petsc/finclude/petscpc.h>
16: #include <petsc/finclude/petscksp.h>
18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
19: ! Variable declarations
20: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
21: Vec x,b,u
22: Mat A
23: PC pc
24: KSP ksp
25: PetscScalar v,one,neg_one
26: PetscReal norm,tol
27: PetscErrorCode ierr
28: PetscInt i,j,II,JJ,Istart
29: PetscInt Iend,m,n,i1,its,five
30: PetscMPIInt rank
31: PetscBool flg
32: #if defined(PETSC_HAVE_MUMPS)
33: Mat F
34: PetscInt ival,icntl,infog34
35: PetscReal cntl,rinfo12,rinfo13,val
36: #endif
38: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39: ! Beginning of program
40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
42: one = 1.0
43: neg_one = -1.0
44: i1 = 1
45: m = 8
46: n = 7
47: five = 5
48: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
49: & '-m',m,flg,ierr)
50: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
51: & '-n',n,flg,ierr)
52: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
54: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: ! Compute the matrix and right-hand-side vector that define
56: ! the linear system, Ax = b.
57: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: call MatCreate(PETSC_COMM_WORLD,A,ierr)
59: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
60: call MatSetType(A, MATAIJ,ierr)
61: call MatSetFromOptions(A,ierr)
62: call MatMPIAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,five, &
63: & PETSC_NULL_INTEGER,ierr)
64: call MatSeqAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,ierr)
66: call MatGetOwnershipRange(A,Istart,Iend,ierr)
68: ! Set matrix elements for the 2-D, five-point stencil in parallel.
69: ! - Each processor needs to insert only elements that it owns
70: ! locally (but any non-local elements will be sent to the
71: ! appropriate processor during matrix assembly).
72: ! - Always specify global row and columns of matrix entries.
73: ! - Note that MatSetValues() uses 0-based row and column numbers
74: ! in Fortran as well as in C.
75: do 10, II=Istart,Iend-1
76: v = -1.0
77: i = II/n
78: j = II - i*n
79: if (i.gt.0) then
80: JJ = II - n
81: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
82: endif
83: if (i.lt.m-1) then
84: JJ = II + n
85: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
86: endif
87: if (j.gt.0) then
88: JJ = II - 1
89: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
90: endif
91: if (j.lt.n-1) then
92: JJ = II + 1
93: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
94: endif
95: v = 4.0
96: call MatSetValues(A,i1,II,i1,II,v,ADD_VALUES,ierr)
97: 10 continue
99: ! Assemble matrix, using the 2-step process:
100: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
101: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
103: ! Create parallel vectors.
104: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,u,ierr)
105: call VecDuplicate(u,b,ierr)
106: call VecDuplicate(b,x,ierr)
108: ! Set exact solution; then compute right-hand-side vector.
109: call VecSet(u,one,ierr)
110: call MatMult(A,u,b,ierr)
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: ! Create the linear solver and set various options
114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
116: call KSPSetOperators(ksp,A,A,ierr)
117: tol = 1.e-7
118: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, &
119: & PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
121: ! Test MUMPS
122: #if defined(PETSC_HAVE_MUMPS)
123: call KSPSetType(ksp,KSPPREONLY,ierr)
124: call KSPGetPC(ksp,pc,ierr)
125: call PCSetType(pc,PCLU,ierr)
126: call PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS,ierr)
127: call PCFactorSetUpMatSolverPackage(pc,ierr)
128: call PCFactorGetMatrix(pc,F,ierr)
130: ! sequential ordering
131: icntl = 7
132: ival = 2
133: call MatMumpsSetIcntl(F,icntl,ival,ierr)
135: ! threshhold for row pivot detection
136: call MatMumpsSetIcntl(F,24,1,ierr)
137: icntl = 3
138: val = 1.e-6
139: call MatMumpsSetCntl(F,icntl,val,ierr)
141: ! compute determinant of A
142: call MatMumpsSetIcntl(F,33,1,ierr)
143: #endif
145: call KSPSetFromOptions(ksp,ierr)
146: call KSPSetUp(ksp,ierr)
147: #if defined(PETSC_HAVE_MUMPS)
148: icntl = 3;
149: call MatMumpsGetCntl(F,icntl,cntl,ierr)
150: call MatMumpsGetInfog(F,34,infog34,ierr)
151: call MatMumpsGetRinfog(F,12,rinfo12,ierr)
152: call MatMumpsGetRinfog(F,13,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 threshhold = ',1pe11.2)
158: 99 format('Mumps determinant=(',1pe11.2,1pe11.2,')*2^',i3)
159: #endif
160:
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