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