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>
  6: program main
  7:   use petscmat
  8:   implicit none

 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

 22:   PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
 23:   m = 10
 24:   one = 1.0
 25:   ione = 1
 26:   ifive = 5

 28:   wmumps = PETSC_FALSE

 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))

 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))

 40:   PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))

 42:   do 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))
 64:   end do

 66:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 67:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

 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))

 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))

 90:   PetscCallA(MatLUFactorNumeric(fact, A, info, ierr))
 91:   PetscCallA(MatSolve(fact, b, x, ierr))
 92:   PetscCallA(MatDestroy(fact, ierr))

 94:   PetscCallA(MatDestroy(A, ierr))
 95:   PetscCallA(VecDestroy(u, ierr))
 96:   PetscCallA(VecDestroy(x, ierr))
 97:   PetscCallA(VecDestroy(b, ierr))

 99:   PetscCallA(PetscFinalize(ierr))
100: end

102: !/*TEST
103: !
104: !   test:
105: !
106: !   test:
107: !     suffix: 2
108: !     args: -use_mumps
109: !     requires: mumps
110: !
111: !TEST*/