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


2
#include "huti_fdefs.h"

3


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

5
!******************************************************************

6
!

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

8
!

9


10
program hutiexample

11


12
  use globals_module

13
  implicit none

14


15
  external huti_d_tfqmr, huti_d_cg

16
  external own_matvec

17


18
  ! local variables

19


20
  character (len=64) :: filename

21
  integer :: ndim, i, j

22
  double precision :: value

23


24
  ! Variables needed for HUTI

25


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

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

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

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

30


31
  integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar

32
  double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar

33


34
  integer :: workdim

35
  parameter ( workdim = 9 )

36


37
  !******************************************************************

38


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

40
  read (5, *) filename

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

42
  read (5, *) ndim

43


44
  allocate( A(ndim, ndim) )

45
  A_ptr => A

46


47
  allocate( b(ndim) )

48
  allocate( x(ndim) )

49
  allocate( work(ndim, workdim) )

50


51
  A = 0; b = 1

52


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

54


55
10 continue

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

57
  A(i, j) = value

58
  go to 10

59


60
20 close(11)

61


62
  HUTI_NDIM = ndim

63
  HUTI_WRKDIM = workdim

64
  HUTI_DBUGLVL = HUTI_ITEROUTPUT

65
  HUTI_MAXIT = HUTI_DFLTMAXIT

66
  HUTI_TOLERANCE = HUTI_DFLTTOLERANCE

67


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

69


70
  print *, 'HUTI returned ', HUTI_INFO

71


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

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

74
  close( 11 )

75


76
end program hutiexample

77


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

79
!******************************************************************

80
!

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

82
! solvers.

83
!

84


85
subroutine own_matvec ( u, v, ipar )

86


87
  use globals_module

88
  implicit none

89


90
  ! Parameters

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

92
  integer, dimension(*) :: ipar

93


94
  ! Local variables

95
  integer :: i,j

96
  double precision :: dsum

97


98
  do i = 1,HUTI_NDIM

99
     dsum = 0

100
     do j = 1,HUTI_NDIM

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

102
     end do

103
     v(i) = dsum

104
  end do

105


106
end subroutine own_matvec

107


108
!******************************************************************

109


110


111