Actual source code: ex126f.F90

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

  6:       program main
  7: #include <petsc/finclude/petscmat.h>
  8:       use petscmat
  9:       implicit none

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

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

 29:       wmumps = PETSC_FALSE

 31:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg, ierr))
 32:       PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-use_mumps',wmumps,flg,ierr))

 34:       PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
 35:       PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*m, m*m, ierr))
 36:       PetscCallA(MatSetType(A, MATAIJ, ierr))
 37:       PetscCallA(MatSetFromOptions(A, ierr))
 38:       PetscCallA(MatSeqAIJSetPreallocation(A,ifive, PETSC_NULL_INTEGER_ARRAY, ierr))
 39:       PetscCallA(MatMPIAIJSetPreallocation(A,ifive,PETSC_NULL_INTEGER_ARRAY,ifive,PETSC_NULL_INTEGER_ARRAY,ierr))

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

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

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

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

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

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

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

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

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