1! 2! This program is modified from a user's contribution. 3! It illustrates how to PetscCallA(MUMPS's LU solver 4! 5 6program main 7#include <petsc/finclude/petscmat.h> 8 use petscmat 9 implicit none 10 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 MatFactorInfo info 22 23 PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr)) 24 m = 10 25 one = 1.0 26 ione = 1 27 ifive = 5 28 29 wmumps = PETSC_FALSE 30 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)) 33 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_ARRAY, ierr)) 39 PetscCallA(MatMPIAIJSetPreallocation(A, ifive, PETSC_NULL_INTEGER_ARRAY, ifive, PETSC_NULL_INTEGER_ARRAY, ierr)) 40 41 PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr)) 42 43 do 10, II = Istart, Iend - 1 44 v = -1.0 45 i = II/m 46 j = II - i*m 47 if (i > 0) then 48 JJ = II - m 49 PetscCallA(MatSetValues(A, ione, [II], ione, [JJ], [v], INSERT_VALUES, ierr)) 50 end if 51 if (i < m - 1) then 52 JJ = II + m 53 PetscCallA(MatSetValues(A, ione, [II], ione, [JJ], [v], INSERT_VALUES, ierr)) 54 end if 55 if (j > 0) then 56 JJ = II - 1 57 PetscCallA(MatSetValues(A, ione, [II], ione, [JJ], [v], INSERT_VALUES, ierr)) 58 end if 59 if (j < m - 1) then 60 JJ = II + 1 61 PetscCallA(MatSetValues(A, ione, [II], ione, [JJ], [v], INSERT_VALUES, ierr)) 62 end if 63 v = 4.0 64 PetscCallA(MatSetValues(A, ione, [II], ione, [II], [v], INSERT_VALUES, ierr)) 6510 continue 66 67 PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) 68 PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) 69 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)) 77 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 end if 87 PetscCallA(MatLUFactorSymbolic(fact, A, perm, iperm, info, ierr)) 88 PetscCallA(ISDestroy(perm, ierr)) 89 PetscCallA(ISDestroy(iperm, ierr)) 90 91 PetscCallA(MatLUFactorNumeric(fact, A, info, ierr)) 92 PetscCallA(MatSolve(fact, b, x, ierr)) 93 PetscCallA(MatDestroy(fact, ierr)) 94 95 PetscCallA(MatDestroy(A, ierr)) 96 PetscCallA(VecDestroy(u, ierr)) 97 PetscCallA(VecDestroy(x, ierr)) 98 PetscCallA(VecDestroy(b, ierr)) 99 100 PetscCallA(PetscFinalize(ierr)) 101 end 102 103!/*TEST 104! 105! test: 106! 107! test: 108! suffix: 2 109! args: -use_mumps 110! requires: mumps 111! 112!TEST*/ 113