1! 2! This program is modified from a user's contribution. 3! It illustrates how to PetscCallA(MUMPS's LU solver 4! 5 6 program 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.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 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 endif 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