Actual source code: ex126f.F

petsc-3.12.0 2019-09-29
Report Typos and Errors
  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:  #include <petsc/finclude/petscmat.h>
  8:       use petscmat
  9:       implicit none


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

 24:       call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
 25:       if (ierr .ne. 0) then
 26:         print*,'Unable to initialize PETSc'
 27:         stop
 28:       endif
 29:       m    = 10
 30:       one  = 1.0
 31:       ione = 1
 32:       ifive = 5

 34:       wmumps = PETSC_FALSE

 36:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,                   &
 37:      &                        '-m',m,flg, ierr)
 38:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,                  &
 39:      &                         '-use_mumps',wmumps,flg,ierr)

 41:       call MatCreate(PETSC_COMM_WORLD, A, ierr)
 42:       call MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*m, m*m, ierr)
 43:       call MatSetType(A, MATAIJ, ierr)
 44:       call MatSetFromOptions(A, ierr)
 45:       call MatSeqAIJSetPreallocation(A,ifive, PETSC_NULL_INTEGER, ierr)
 46:       call MatMPIAIJSetPreallocation(A,ifive,PETSC_NULL_INTEGER,ifive,  &
 47:      &     PETSC_NULL_INTEGER,ierr)

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

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

 75:       call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)
 76:       call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)

 78:       call VecCreate(PETSC_COMM_WORLD, u, ierr)
 79:       call VecSetSizes(u, PETSC_DECIDE, m*m, ierr)
 80:       call VecSetFromOptions(u, ierr)
 81:       call VecDuplicate(u,b,ierr)
 82:       call VecDuplicate(b,x,ierr)
 83:       call VecSet(u, one, ierr)
 84:       call MatMult(A, u, b, ierr)

 86:       call MatFactorInfoInitialize(info,ierr)
 87:       call MatGetOrdering(A,MATORDERINGNATURAL,perm,iperm,ierr)
 88:       if (wmumps) then
 89:           write(*,*) 'use MUMPS LU...'
 90:           call MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,fact,ierr)
 91:       else
 92:          write(*,*) 'use PETSc LU...'
 93:          call MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,fact,ierr)
 94:       endif
 95:       call MatLUFactorSymbolic(fact, A, perm, iperm,                                   &
 96:      &         info, ierr)
 97:       call ISDestroy(perm,ierr)
 98:       call ISDestroy(iperm,ierr)

100:       call MatLUFactorNumeric(fact, A, info, ierr)
101:       call MatSolve(fact, b, x,ierr)
102:       call MatDestroy(fact, ierr)

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

109:       call PetscFinalize(ierr)
110:       end

112: !/*TEST
113: !
114: !   test:
115: !
116: !   test:
117: !     suffix: 2
118: !     args: -use_mumps
119: !     requires: mumps
120: !
121: !TEST*/