Actual source code: ex126f.F

  1: !
  2: ! This program is modified from a user's contribution.
  3: ! It illustrates how to call MUMPS's LU solver
  4: !

  6:       program main
  7:       implicit none

 9:  #include finclude/petscsys.h
 10:  #include finclude/petscvec.h
 11:  #include finclude/petscmat.h

 13:       Vec            x,b,u
 14:       Mat            A, fact
 15:       PetscInt       i,j,II,JJ,m
 16:       PetscInt       Istart, Iend
 17:       PetscInt       ione
 18:       PetscTruth     wmumps
 19:       PetscTruth     flg
 20:       PetscScalar    one, v
 21:       IS             perm,iperm
 22:       PetscErrorCode ierr

 24:       call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
 25:       m    = 10
 26:       one  = 1.0
 27:       ione = 1.0

 29:       wmumps = PETSC_FALSE

 31:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg, ierr)
 32:       call PetscOptionsGetTruth(PETSC_NULL_CHARACTER,'-use_mumps',
 33:      > wmumps,flg,ierr)

 35:       call MatCreate(PETSC_COMM_WORLD, A, ierr)
 36:       call MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*m, m*m, ierr)
 37:       call MatSetType(A, MATAIJ, ierr)
 38:       call MatSetFromOptions(A, ierr)
 39:       call MatSeqAIJSetPreallocation(A, 1, PETSC_NULL_INTEGER, ierr)
 40:       call MatMPIAIJSetPreallocation(A,5,PETSC_NULL_INTEGER,5,
 41:      >     PETSC_NULL_INTEGER,ierr)

 43:       call MatGetOwnershipRange(A,Istart,Iend,ierr)

 45:       do 10, II=Istart,Iend - 1
 46:         v = -1.0
 47:         i = II/m
 48:         j = II - i*m
 49:         if (i.gt.0) then
 50:           JJ = II - m
 51:           call MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr)
 52:         endif
 53:         if (i.lt.m-1) then
 54:           JJ = II + m
 55:           call MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr)
 56:         endif
 57:         if (j.gt.0) then
 58:           JJ = II - 1
 59:           call MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr)
 60:         endif
 61:         if (j.lt.m-1) then
 62:           JJ = II + 1
 63:           call MatSetValues(A,ione,II,ione,JJ,v,INSERT_VALUES,ierr)
 64:         endif
 65:         v = 4.0
 66:         call  MatSetValues(A,ione,II,ione,II,v,INSERT_VALUES,ierr)
 67:  10   continue
 68: 
 69:       call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)
 70:       call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)

 72:       call VecCreate(PETSC_COMM_WORLD, u, ierr)
 73:       call VecSetSizes(u, PETSC_DECIDE, m*m, ierr)
 74:       call VecSetFromOptions(u, ierr)
 75:       call VecDuplicate(u,b,ierr)
 76:       call VecDuplicate(b,x,ierr)
 77:       call VecSet(u, one, ierr)
 78:       call MatMult(A, u, b, ierr)

 80:       if (wmumps) then
 81:           write(*,*) 'use MUMPS LU...'
 82:           call MatGetFactor(A, MAT_SOLVER_MUMPS, MAT_FACTOR_LU,
 83:      >         fact, ierr)
 84:           call MatLUFactorSymbolic(fact, A, PETSC_NULL, PETSC_NULL,
 85:      >         PETSC_NULL, ierr)
 86:       else
 87:          write(*,*) 'use PETSc LU...'
 88:          call MatGetOrdering(A,MATORDERING_NATURAL,perm,iperm,ierr)
 89:          call MatGetFactor(A, MAT_SOLVER_PETSC, MAT_FACTOR_LU,
 90:      >         fact, ierr)
 91: 
 92:          call MatLUFactorSymbolic(fact, A, perm, iperm,
 93:      >         PETSC_NULL, ierr)
 94:          call ISDestroy(perm,ierr)
 95:          call ISDestroy(iperm,ierr)
 96:       endif

 98:       call MatLUFactorNumeric(fact, A, PETSC_NULL, ierr)
 99:       call MatSolve(fact, b, x,ierr)
100:       call MatDestroy(fact, ierr)
101: 
102:       call MatDestroy(A, ierr)
103:       call VecDestroy(u, ierr)
104:       call VecDestroy(x, ierr)
105:       call VecDestroy(b, ierr)

107:       call PetscFinalize(ierr)
108:       end