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
 13:   PetscInt Istart, Iend
 14:   PetscInt m
 15:   PetscBool wmumps
 16:   PetscBool flg
 17:   PetscScalar, parameter :: one = 1.0
 18:   PetscScalar v
 19:   IS perm, iperm
 20:   PetscErrorCode ierr
 21:   MatFactorInfo info

 23:   PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))

 25:   wmumps = PETSC_FALSE

 27:   m = 10
 28:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
 29:   PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-use_mumps', wmumps, flg, ierr))

 31:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
 32:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m**2, m**2, ierr))
 33:   PetscCallA(MatSetType(A, MATAIJ, ierr))
 34:   PetscCallA(MatSetFromOptions(A, ierr))
 35:   PetscCallA(MatSeqAIJSetPreallocation(A, 5_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
 36:   PetscCallA(MatMPIAIJSetPreallocation(A, 5_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 5_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))

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

 40:   do II = Istart, Iend - 1
 41:     v = -1.0
 42:     i = II/m
 43:     j = II - i*m
 44:     if (i > 0) then
 45:       JJ = II - m
 46:       PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
 47:     end if
 48:     if (i < m - 1) then
 49:       JJ = II + m
 50:       PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
 51:     end if
 52:     if (j > 0) then
 53:       JJ = II - 1
 54:       PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
 55:     end if
 56:     if (j < m - 1) then
 57:       JJ = II + 1
 58:       PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
 59:     end if
 60:     v = 4.0
 61:     PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [II], [v], INSERT_VALUES, ierr))
 62:   end do

 64:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 65:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

 67:   PetscCallA(VecCreate(PETSC_COMM_WORLD, u, ierr))
 68:   PetscCallA(VecSetSizes(u, PETSC_DECIDE, m*m, ierr))
 69:   PetscCallA(VecSetFromOptions(u, ierr))
 70:   PetscCallA(VecDuplicate(u, b, ierr))
 71:   PetscCallA(VecDuplicate(b, x, ierr))
 72:   PetscCallA(VecSet(u, one, ierr))
 73:   PetscCallA(MatMult(A, u, b, ierr))

 75:   PetscCallA(MatFactorInfoInitialize(info, ierr))
 76:   PetscCallA(MatGetOrdering(A, MATORDERINGNATURAL, perm, iperm, ierr))
 77:   if (wmumps) then
 78:     write (*, *) 'use MUMPS LU...'
 79:     PetscCallA(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, fact, ierr))
 80:   else
 81:     write (*, *) 'use PETSc LU...'
 82:     PetscCallA(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, fact, ierr))
 83:   end if
 84:   PetscCallA(MatLUFactorSymbolic(fact, A, perm, iperm, info, ierr))
 85:   PetscCallA(ISDestroy(perm, ierr))
 86:   PetscCallA(ISDestroy(iperm, ierr))

 88:   PetscCallA(MatLUFactorNumeric(fact, A, info, ierr))
 89:   PetscCallA(MatSolve(fact, b, x, ierr))
 90:   PetscCallA(MatDestroy(fact, ierr))

 92:   PetscCallA(MatDestroy(A, ierr))
 93:   PetscCallA(VecDestroy(u, ierr))
 94:   PetscCallA(VecDestroy(x, ierr))
 95:   PetscCallA(VecDestroy(b, ierr))

 97:   PetscCallA(PetscFinalize(ierr))
 98: end

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