Actual source code: ex52f.F
petsc-3.5.4 2015-05-23
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 <finclude/petscsys.h>
13: #include <finclude/petscvec.h>
14: #include <finclude/petscmat.h>
15: #include <finclude/petscpc.h>
16: #include <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_CHARACTER,'-m',m,flg,ierr)
49: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
50: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: ! Compute the matrix and right-hand-side vector that define
54: ! the linear system, Ax = b.
55: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: call MatCreate(PETSC_COMM_WORLD,A,ierr)
57: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
58: call MatSetType(A, MATAIJ,ierr)
59: call MatSetFromOptions(A,ierr)
60: call MatMPIAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,five, &
61: & PETSC_NULL_INTEGER,ierr)
62: call MatSeqAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,ierr)
64: call MatGetOwnershipRange(A,Istart,Iend,ierr)
66: ! Set matrix elements for the 2-D, five-point stencil in parallel.
67: ! - Each processor needs to insert only elements that it owns
68: ! locally (but any non-local elements will be sent to the
69: ! appropriate processor during matrix assembly).
70: ! - Always specify global row and columns of matrix entries.
71: ! - Note that MatSetValues() uses 0-based row and column numbers
72: ! in Fortran as well as in C.
73: do 10, II=Istart,Iend-1
74: v = -1.0
75: i = II/n
76: j = II - i*n
77: if (i.gt.0) then
78: JJ = II - n
79: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
80: endif
81: if (i.lt.m-1) then
82: JJ = II + n
83: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
84: endif
85: if (j.gt.0) then
86: JJ = II - 1
87: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
88: endif
89: if (j.lt.n-1) then
90: JJ = II + 1
91: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
92: endif
93: v = 4.0
94: call MatSetValues(A,i1,II,i1,II,v,ADD_VALUES,ierr)
95: 10 continue
97: ! Assemble matrix, using the 2-step process:
98: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
99: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
101: ! Create parallel vectors.
102: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,u,ierr)
103: call VecDuplicate(u,b,ierr)
104: call VecDuplicate(b,x,ierr)
106: ! Set exact solution; then compute right-hand-side vector.
107: call VecSet(u,one,ierr)
108: call MatMult(A,u,b,ierr)
110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: ! Create the linear solver and set various options
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
114: call KSPSetOperators(ksp,A,A,ierr)
115: tol = 1.e-7
116: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, &
117: & PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
119: ! Test MUMPS
120: #if defined(PETSC_HAVE_MUMPS)
121: call KSPSetType(ksp,KSPPREONLY,ierr)
122: call KSPGetPC(ksp,pc,ierr)
123: call PCSetType(pc,PCLU,ierr)
124: call PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS,ierr)
125: call PCFactorSetUpMatSolverPackage(pc,ierr)
126: call PCFactorGetMatrix(pc,F,ierr)
128: ! sequential ordering
129: icntl = 7
130: ival = 2
131: call MatMumpsSetIcntl(F,icntl,ival,ierr)
133: ! threshhold for row pivot detection
134: call MatMumpsSetIcntl(F,24,1,ierr)
135: icntl = 3
136: val = 1.e-6
137: call MatMumpsSetCntl(F,icntl,val,ierr)
139: ! compute determinant of A
140: call MatMumpsSetIcntl(F,33,1,ierr)
141: #endif
143: call KSPSetFromOptions(ksp,ierr)
144: call KSPSetUp(ksp,ierr)
145: #if defined(PETSC_HAVE_MUMPS)
146: icntl = 3;
147: call MatMumpsGetCntl(F,icntl,cntl,ierr)
148: call MatMumpsGetInfog(F,34,infog34,ierr)
149: call MatMumpsGetRinfog(F,12,rinfo12,ierr)
150: call MatMumpsGetRinfog(F,13,rinfo13,ierr)
151: if (rank .eq. 0) then
152: write(6,98) cntl
153: write(6,99) rinfo12,rinfo13,infog34
154: endif
155: 98 format('Mumps row pivot threshhold = ',1pe11.2)
156: 99 format('Mumps determinant=(',1pe11.2,1pe11.2,')*2^',i3)
157: #endif
158:
159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: ! Solve the linear system
161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: call KSPSolve(ksp,b,x,ierr)
164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: ! Check solution and clean up
166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: call VecAXPY(x,neg_one,u,ierr)
168: call VecNorm(x,NORM_2,norm,ierr)
169: call KSPGetIterationNumber(ksp,its,ierr)
171: if (rank .eq. 0) then
172: if (norm .gt. 1.e-12) then
173: write(6,100) norm,its
174: else
175: write(6,110) its
176: endif
177: endif
178: 100 format('Norm of error ',1pe11.4,' iterations ',i5)
179: 110 format('Norm of error < 1.e-12,iterations ',i5)
181: ! Free work space. All PETSc objects should be destroyed when they
182: ! are no longer needed.
184: call KSPDestroy(ksp,ierr)
185: call VecDestroy(u,ierr)
186: call VecDestroy(x,ierr)
187: call VecDestroy(b,ierr)
188: call MatDestroy(A,ierr)
190: ! Always call PetscFinalize() before exiting a program.
191: call PetscFinalize(ierr)
192: end