Actual source code: ex126f.F90
1: !
2: ! This program is modified from a user's contribution.
3: ! It illustrates how to PetscCallA(MUMPS's LU solver
4: !
6: program main
7: #include <petsc/finclude/petscmat.h>
8: use petscmat
9: implicit none
11: Vec x,b,u
12: Mat A, fact
13: PetscInt i,j,II,JJ,m
14: PetscInt Istart, Iend
15: PetscInt ione, ifive
16: PetscBool wmumps
17: PetscBool flg
18: PetscScalar one, v
19: IS perm,iperm
20: PetscErrorCode ierr
21: PetscReal info(MAT_FACTORINFO_SIZE)
23: PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
24: m = 10
25: one = 1.0
26: ione = 1
27: ifive = 5
29: wmumps = PETSC_FALSE
31: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg, ierr))
32: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-use_mumps',wmumps,flg,ierr))
34: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
35: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*m, m*m, ierr))
36: PetscCallA(MatSetType(A, MATAIJ, ierr))
37: PetscCallA(MatSetFromOptions(A, ierr))
38: PetscCallA(MatSeqAIJSetPreallocation(A,ifive, PETSC_NULL_INTEGER, ierr))
39: PetscCallA(MatMPIAIJSetPreallocation(A,ifive,PETSC_NULL_INTEGER,ifive,PETSC_NULL_INTEGER,ierr))
41: PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
43: do 10, II=Istart,Iend - 1
44: v = -1.0
45: i = II/m
46: j = II - i*m
47: if (i.gt.0) then
48: JJ = II - m
49: PetscCallA(MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr))
50: endif
51: if (i.lt.m-1) then
52: JJ = II + m
53: PetscCallA(MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr))
54: endif
55: if (j.gt.0) then
56: JJ = II - 1
57: PetscCallA(MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr))
58: endif
59: if (j.lt.m-1) then
60: JJ = II + 1
61: PetscCallA(MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr))
62: endif
63: v = 4.0
64: PetscCallA( MatSetValues(A,ione,II,ione,II,v,INSERT_VALUES,ierr))
65: 10 continue
67: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
68: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
70: PetscCallA(VecCreate(PETSC_COMM_WORLD, u, ierr))
71: PetscCallA(VecSetSizes(u, PETSC_DECIDE, m*m, ierr))
72: PetscCallA(VecSetFromOptions(u, ierr))
73: PetscCallA(VecDuplicate(u,b,ierr))
74: PetscCallA(VecDuplicate(b,x,ierr))
75: PetscCallA(VecSet(u, one, ierr))
76: PetscCallA(MatMult(A, u, b, ierr))
78: PetscCallA(MatFactorInfoInitialize(info,ierr))
79: PetscCallA(MatGetOrdering(A,MATORDERINGNATURAL,perm,iperm,ierr))
80: if (wmumps) then
81: write(*,*) 'use MUMPS LU...'
82: PetscCallA(MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,fact,ierr))
83: else
84: write(*,*) 'use PETSc LU...'
85: PetscCallA(MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,fact,ierr))
86: endif
87: PetscCallA(MatLUFactorSymbolic(fact, A, perm, iperm, info, ierr))
88: PetscCallA(ISDestroy(perm,ierr))
89: PetscCallA(ISDestroy(iperm,ierr))
91: PetscCallA(MatLUFactorNumeric(fact, A, info, ierr))
92: PetscCallA(MatSolve(fact, b, x,ierr))
93: PetscCallA(MatDestroy(fact, ierr))
95: PetscCallA(MatDestroy(A, ierr))
96: PetscCallA(VecDestroy(u, ierr))
97: PetscCallA(VecDestroy(x, ierr))
98: PetscCallA(VecDestroy(b, ierr))
100: PetscCallA(PetscFinalize(ierr))
101: end
103: !/*TEST
104: !
105: ! test:
106: !
107: ! test:
108: ! suffix: 2
109: ! args: -use_mumps
110: ! requires: mumps
111: !
112: !TEST*/