Choral
gmres_mod.f90
Go to the documentation of this file.
1 
8 
9 module gmres_mod
10 
11  use real_type
12  use basic_tools
13  use abstract_interfaces, only: rntorn
14  use algebra_lin, only: invtrisup
15  use r1d_mod, only: norm_2, axpy, scalprod, xpay
16 
17  implicit none
18  private
19 
20  public :: gmres !! TESTED
21 
22 contains
23 
24 
47  subroutine gmres(x, iter, nbPrd, res, &
48  & b, A, tol, itmax, rst, verb)
49 
50  real(RP), dimension(:), intent(inout) :: x
51  integer , intent(out) :: iter, nbPrd
52  real(RP) , intent(out) :: res
53  procedure(rntorn) :: A
54  real(RP), dimension(:), intent(in) :: b
55  real(RP) , intent(in) :: tol
56  integer , intent(in) :: itmax, rst, verb
57 
58  real(RP), dimension(rst+1, size(x,1) ) :: V
59  real(RP), dimension(rst+1, rst ) :: H
60 
61  real(RP), dimension(size(x,1)) :: r, w
62  real(RP), dimension(rst ) :: sn, cs, y
63  real(RP), dimension(rst +1 ) :: s
64 
65  real(RP) :: nb2, nr2, temp
66  integer :: nn, ii, kk
67 
68  if (verb>1) write(*,*) 'gmres_mod : gmres'
69 
70  iter = 0
71  nbprd = 0
72 
73 
74  nb2 = norm_2(b)
75  nn = size(x,1)
76  if (nb2/re(nn)<1e-8_rp) nb2=sqrt(re(nn))
77  nb2 = 1._rp/nb2
78 
79  call a(r, x)
80  nbprd = 1
81  r = b-r
82  nr2 = norm_2(r)
83  res = nr2 * nb2
84 
85  if (verb>2) write(*,*)' iter', iter,' residual', res
86  if (res<tol) return
87 
88  v = 0._rp
89  h = 0._rp
90  cs = 0._rp
91  sn = 0._rp
92 
93  do iter=1, itmax
94 
95  if (verb>2) write(*,*)' iter', iter,' residual', res
96 
97  s = 0._rp
98  s(1) = nr2
99 
100  v(1,:) = r / nr2
101 
102  do ii = 1, rst ! orthonormal basis using Gram-Schmidt
103 
104  call a(w, v(ii,:))
105  nbprd = nbprd + 1
106 
107  do kk = 1, ii
108  h(kk,ii)= scalprod( w, v(kk,:))
109  w = w - h(kk,ii)*v(kk,:)
110  end do
111 
112  h(ii+1,ii) = norm_2(w)
113 
114  v(ii+1, : ) = w / h(ii+1,ii)
115 
116  ! apply Givens rotation
117  !
118  do kk = 1, ii-1
119  temp = cs(kk)*h(kk,ii) + sn(kk)*h(kk+1,ii)
120  h(kk+1,ii) = -sn(kk)*h(kk,ii) + cs(kk)*h(kk+1,ii)
121  h(kk,ii) = temp
122  end do
123 
124  call grotmat(cs(ii), sn(ii), h(ii,ii), h(ii+1,ii) )
125 
126  temp = cs(ii)*s(ii) ! approximate residual norm
127  s(ii+1) = -sn(ii)*s(ii)
128  s(ii) = temp
129 
130  h(ii,ii) = cs(ii)*h(ii,ii) + sn(ii)*h(ii+1,ii)
131  h(ii+1,ii) = 0._rp
132 
133  res = abs(s(ii+1)) * nb2
134 
135  if (verb>2) write(*,*)' iter', iter,' residual', res
136 
137  if ( res < tol ) then
138 
139  call invtrisup(y(1:ii), h(1:ii,1:ii), s(1:ii))
140  do kk=1, ii
141  x = x + v(kk, : )*y(kk)
142  end do
143 
144  return
145 
146  end if
147 
148  end do
149 
150  ! update approximation
151  !
152  call invtrisup(y(1:rst), h(1:rst,1:rst), s(1:rst))
153  do kk=1, rst
154  x = x + v(kk, : )*y(kk)
155  end do
156 
157  ! compute residual
158  !
159  call a(r, x)
160  nbprd = nbprd + 1
161  r = b-r
162  nr2 = norm_2(r)
163  res = nr2 * nb2
164  if (verb>2) write(*,*)' iter', iter,' residual', res
165 
166  if ( res < tol ) return
167 
168  end do
169 
170  ! Convergence failure flag
171  call warning("gmres_mod: gmres: not converged", verb-1)
172  iter = -1
173 
174  end subroutine gmres
175 
177  subroutine grotmat(cs, sn, a, b)
179  real(RP), intent(in) :: a,b
180  real(RP), intent(out) :: cs,sn
181 
182  real(RP) :: tmp
183 
184  if (abs(b)<real_tol) then
185  cs = 1._rp
186  sn = 0._rp
187 
188  else if (abs(b)>abs(a)) then
189  tmp = a/b
190  sn = 1._rp / sqrt( 1._rp + tmp**2)
191  cs = tmp*sn
192 
193  else
194  tmp = b/a
195  cs = 1._rp / sqrt( 1._rp + tmp**2)
196  sn = tmp*cs
197 
198  end if
199 
200  end subroutine grotmat
201 
202 end module gmres_mod
203 
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
LINEAR ALGEBRA TOOLS
Definition: algebra_lin.F90:9
BASIC TOOLS
Definition: basic_tools.F90:8
GMRES LINEAR SOLVER
Definition: gmres_mod.f90:9
L2 vector norm.
Definition: R1d_mod.F90:82
subroutine, public gmres(x, iter, nbPrd, res, b, A, tol, itmax, rst, verb)
GMRES (no preconditioning)
Definition: gmres_mod.f90:49
conversion integers or rational to real
Definition: real_type.F90:153
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
Definition: R1d_mod.F90:10
subroutine, public invtrisup(res, M, vec)
Calcul de x=M^{-1}vec avec M triangulaire superieure.
real(rp) function e(x, v1, v2)
x = a*x + y
Definition: R1d_mod.F90:55
scalar product
Definition: R1d_mod.F90:77
x = x + b*y // OR // z = x + b*y
Definition: R1d_mod.F90:60
subroutine grotmat(cs, sn, a, b)
Matrice de rotation de Givens.
Definition: gmres_mod.f90:178
real(rp), parameter, public real_tol
Definition: real_type.F90:91
subroutine, public warning(message, verb)
Warning message.