Actual source code: ex126f.F

petsc-3.4.2 2013-07-02
  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, ifive
 18:       PetscBool      wmumps
 19:       PetscBool      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
 28:       ifive = 5

 30:       wmumps = PETSC_FALSE

 32:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg, ierr)
 33:       call PetscOptionsGetBool(PETSC_NULL_CHARACTER,'-use_mumps',                       &
 34:      &                         wmumps,flg,ierr)

 36:       call MatCreate(PETSC_COMM_WORLD, A, ierr)
 37:       call MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*m, m*m, ierr)
 38:       call MatSetType(A, MATAIJ, ierr)
 39:       call MatSetFromOptions(A, ierr)
 40:       call MatSeqAIJSetPreallocation(A,ifive, PETSC_NULL_INTEGER, ierr)
 41:       call MatMPIAIJSetPreallocation(A,ifive,PETSC_NULL_INTEGER,ifive,  &
 42:      &     PETSC_NULL_INTEGER,ierr)

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

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

 70:       call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)
 71:       call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)

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

 81:       if (wmumps) then
 82:           write(*,*) 'use MUMPS LU...'
 83:           call MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU,                             &
 84:      &                      fact, ierr)
 85:           call MatLUFactorSymbolic(fact, A, PETSC_NULL_OBJECT,                            &
 86:      &                    PETSC_NULL_OBJECT,PETSC_NULL_INTEGER, ierr)
 87:       else
 88:          write(*,*) 'use PETSc LU...'
 89:          call MatGetOrdering(A,MATORDERINGNATURAL,perm,iperm,ierr)
 90:          call MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU,                              &
 91:      &         fact, ierr)

 93:          call MatLUFactorSymbolic(fact, A, perm, iperm,                                   &
 94:      &         PETSC_NULL_INTEGER, ierr)
 95:          call ISDestroy(perm,ierr)
 96:          call ISDestroy(iperm,ierr)
 97:       endif

 99:       call MatLUFactorNumeric(fact, A, PETSC_NULL_INTEGER, ierr)
100:       call MatSolve(fact, b, x,ierr)
101:       call MatDestroy(fact, ierr)

103:       call MatDestroy(A, ierr)
104:       call VecDestroy(u, ierr)
105:       call VecDestroy(x, ierr)
106:       call VecDestroy(b, ierr)

108:       call PetscFinalize(ierr)
109:       end