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