Actual source code: ex126f.F
1: !
2: ! This program is modified from a user's contribution.
3: ! It illustrates how to call MUMPS's LU solver
4: !
6: program main
7: implicit none
9: #include finclude/petscsys.h
10: #include finclude/petscvec.h
11: #include finclude/petscmat.h
13: Vec x,b,u
14: Mat A, fact
15: PetscInt i,j,II,JJ,m
16: PetscInt Istart, Iend
17: PetscInt ione
18: PetscTruth wmumps
19: PetscTruth flg
20: PetscScalar one, v
21: IS perm,iperm
22: PetscErrorCode ierr
24: call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
25: m = 10
26: one = 1.0
27: ione = 1.0
29: wmumps = PETSC_FALSE
31: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg, ierr)
32: call PetscOptionsGetTruth(PETSC_NULL_CHARACTER,'-use_mumps',
33: > wmumps,flg,ierr)
35: call MatCreate(PETSC_COMM_WORLD, A, ierr)
36: call MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*m, m*m, ierr)
37: call MatSetType(A, MATAIJ, ierr)
38: call MatSetFromOptions(A, ierr)
39: call MatSeqAIJSetPreallocation(A, 1, PETSC_NULL_INTEGER, ierr)
40: call MatMPIAIJSetPreallocation(A,5,PETSC_NULL_INTEGER,5,
41: > PETSC_NULL_INTEGER,ierr)
43: call MatGetOwnershipRange(A,Istart,Iend,ierr)
45: do 10, II=Istart,Iend - 1
46: v = -1.0
47: i = II/m
48: j = II - i*m
49: if (i.gt.0) then
50: JJ = II - m
51: call MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr)
52: endif
53: if (i.lt.m-1) then
54: JJ = II + m
55: call MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr)
56: endif
57: if (j.gt.0) then
58: JJ = II - 1
59: call MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr)
60: endif
61: if (j.lt.m-1) then
62: JJ = II + 1
63: call MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr)
64: endif
65: v = 4.0
66: call MatSetValues(A,ione,II,ione,II,v,INSERT_VALUES,ierr)
67: 10 continue
68:
69: call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)
70: call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)
72: call VecCreate(PETSC_COMM_WORLD, u, ierr)
73: call VecSetSizes(u, PETSC_DECIDE, m*m, ierr)
74: call VecSetFromOptions(u, ierr)
75: call VecDuplicate(u,b,ierr)
76: call VecDuplicate(b,x,ierr)
77: call VecSet(u, one, ierr)
78: call MatMult(A, u, b, ierr)
80: if (wmumps) then
81: write(*,*) 'use MUMPS LU...'
82: call MatGetFactor(A, MAT_SOLVER_MUMPS, MAT_FACTOR_LU,
83: > fact, ierr)
84: call MatLUFactorSymbolic(fact, A, PETSC_NULL, PETSC_NULL,
85: > PETSC_NULL, ierr)
86: else
87: write(*,*) 'use PETSc LU...'
88: call MatGetOrdering(A,MATORDERING_NATURAL,perm,iperm,ierr)
89: call MatGetFactor(A, MAT_SOLVER_PETSC, MAT_FACTOR_LU,
90: > fact, ierr)
91:
92: call MatLUFactorSymbolic(fact, A, perm, iperm,
93: > PETSC_NULL, ierr)
94: call ISDestroy(perm,ierr)
95: call ISDestroy(iperm,ierr)
96: endif
98: call MatLUFactorNumeric(fact, A, PETSC_NULL, ierr)
99: call MatSolve(fact, b, x,ierr)
100: call MatDestroy(fact, ierr)
101:
102: call MatDestroy(A, ierr)
103: call VecDestroy(u, ierr)
104: call VecDestroy(x, ierr)
105: call VecDestroy(b, ierr)
107: call PetscFinalize(ierr)
108: end