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:   MatFactorInfo info

 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_ARRAY, ierr))
 39:   PetscCallA(MatMPIAIJSetPreallocation(A, ifive, PETSC_NULL_INTEGER_ARRAY, ifive, PETSC_NULL_INTEGER_ARRAY, 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 > 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))
 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:     end if
 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*/