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