testgmres/testgmres.f90 at master · jayten/testgmres · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
!=========================================================================
subroutine noprecond(MinvX,X,n)

      implicit none

      integer,intent(in) :: n
      real*8, intent(in) :: X(n)
      real*8,intent(out) :: MinvX(n)

      MinvX=X

end subroutine noprecond
!=========================================================================
subroutine findAX(AX,X,n)

	implicit none

	integer,intent(in) :: n
        real*8, intent(in) :: X(n)
        real*8,intent(out) :: AX(n)

	real*8 :: sup,sub,diag
	integer :: i

	sup  =  1.0
	sub  =  1.0
	diag = -2.0

	AX = 0.d0

	do i=1,n

		if(i .eq. 1) then
			AX(i)=diag*x(i)+sup*x(i+1)
		else if(i .eq. n) then
			AX(i)=diag*x(i)+sub*x(i-1)
		else
			AX(i)=sub*x(i-1)+diag*x(i)+sup*x(i+1)
		endif
	enddo

end subroutine findAX
!=========================================================================
subroutine findbvec(b,n)

	integer, intent(in) :: n
	real*8, intent(out)  :: b(n)

	b=0.d0

	b(1) = 0.d0
	b(n) =-1.d0

end subroutine findbvec
!=========================================================================
program testgmres
      use solvergmres
      implicit none

      integer :: m,n,it,i
      real*8,allocatable :: x0(:),x(:),b(:)
      external :: findAX,noprecond

      open(unit=5,file="soln.dat")

      m=40
      n=50
      it=20

      allocate(x0(n))
      allocate(x(n))
      allocate(b(n))

      x0 = 0.d0
      x  = 0.d0

      call findbvec(b,n)

      call performgmres(b,x0,x,m,n,it,findAX,noprecond)

      do i=1,n
      	write(5,'(I5 F10.5)'),i,x(i)
      enddo

      close(5)

end program testgmres
!=========================================================================