Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
 Path: blob/devel/fhutiter/examples/ex1/huti-ex.F90
3206 views
1
#include "huti_fdefs.h"

2


3
!******************************************************************

4
!******************************************************************

5
!

6
! Purpose of this example is to show how to use HUTIter library

7
!

8


9
program hutiexample

10


11
  use globals_module

12
  implicit none

13


14
  external huti_d_tfqmr, huti_d_cg

15
  external own_matvec

16


17
  ! local variables

18


19
  character (len=64) :: filename

20
  integer :: ndim, i, j

21
  double precision :: value

22


23
  ! Variables needed for HUTI

24


25
  double precision, dimension(:,:), allocatable, target :: A

26
  double precision, dimension(:), allocatable :: b

27
  double precision, dimension(:), allocatable :: x

28
  double precision, dimension(:,:), allocatable :: work

29


30
  integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar

31
  double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar

32


33
  integer :: workdim

34
  parameter ( workdim = HUTI_TFQMR_WORKSIZE )

35


36
  !******************************************************************

37


38
  write (6, *) 'Enter name of the matrix file :'

39
  read (5, *) filename

40
  write (6, *) 'Enter leading dimension of the matrix :'

41
  read (5, *) ndim

42


43
  allocate( A(ndim, ndim) )

44
  A_ptr => A

45


46
  allocate( b(ndim) )

47
  allocate( x(ndim) )

48
  allocate( work(ndim, workdim) )

49


50
  A = 0; b = 1

51


52
  open ( 11, file=filename, status='old', err=20 )

53


54
10 continue

55
  read ( 11, *, end=20 ) i, j, value

56
  A(i, j) = value

57
  go to 10

58


59
20 close(11)

60


61
  HUTI_NDIM = ndim

62
  HUTI_WRKDIM = workdim

63
  HUTI_DBUGLVL = HUTI_ITEROUTPUT

64
  HUTI_MAXIT = HUTI_DFLTMAXIT

65
  HUTI_TOLERANCE = HUTI_DFLTTOLERANCE

66


67
  call HUTI_D_TFQMR( x, b, ipar, dpar, work, own_matvec, 0, 0, 0, 0, 0 )

68


69
  print *, 'HUTI returned ', HUTI_INFO

70


71
  open( 11, file=trim( filename ) // '.out' )

72
  write( 11, '(I6,A,E25.17)') ( i,' ', x(i), i=1,ndim )

73
  close( 11 )

74


75
end program hutiexample

76


77
!******************************************************************

78
!******************************************************************

79
!

80
! This is an example of the user supplied matvec routine for HUTI

81
! solvers.

82
!

83


84
subroutine own_matvec ( u, v, ipar )

85


86
  use globals_module

87
  implicit none

88


89
  ! Parameters

90
  double precision, dimension(*) :: u, v

91
  integer, dimension(*) :: ipar

92


93
  ! Local variables

94
  integer :: i,j

95
  double precision :: dsum

96


97
  do i = 1,HUTI_NDIM

98
     dsum = 0

99
     do j = 1,HUTI_NDIM

100
	dsum = dsum + A_ptr(i,j) * u(j)

101
     end do

102
     v(i) = dsum

103
  end do

104


105
end subroutine own_matvec

106


107
!******************************************************************

108


109


110