47 subroutine gmres(x, iter, nbPrd, res, &
48 & b, A, tol, itmax, rst, verb)
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
58 real(RP),
dimension(rst+1, size(x,1) ) :: V
59 real(RP),
dimension(rst+1, rst ) :: H
61 real(RP),
dimension(size(x,1)) :: r, w
62 real(RP),
dimension(rst ) :: sn, cs, y
63 real(RP),
dimension(rst +1 ) :: s
65 real(RP) :: nb2, nr2, temp
68 if (verb>1)
write(*,*)
'gmres_mod : gmres' 76 if (nb2/
re(nn)<1
e-8_rp) nb2=sqrt(
re(nn))
85 if (verb>2)
write(*,*)
' iter', iter,
' residual', res
95 if (verb>2)
write(*,*)
' iter', iter,
' residual', res
109 w = w - h(kk,ii)*v(kk,:)
114 v(ii+1, : ) = w / h(ii+1,ii)
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)
124 call grotmat(cs(ii), sn(ii), h(ii,ii), h(ii+1,ii) )
127 s(ii+1) = -sn(ii)*s(ii)
130 h(ii,ii) = cs(ii)*h(ii,ii) + sn(ii)*h(ii+1,ii)
133 res = abs(s(ii+1)) * nb2
135 if (verb>2)
write(*,*)
' iter', iter,
' residual', res
137 if ( res < tol )
then 139 call invtrisup(y(1:ii), h(1:ii,1:ii), s(1:ii))
141 x = x + v(kk, : )*y(kk)
152 call invtrisup(y(1:rst), h(1:rst,1:rst), s(1:rst))
154 x = x + v(kk, : )*y(kk)
164 if (verb>2)
write(*,*)
' iter', iter,
' residual', res
166 if ( res < tol )
return 171 call warning(
"gmres_mod: gmres: not converged", verb-1)
177 subroutine grotmat(cs, sn, a, b)
179 real(RP),
intent(in) :: a,b
180 real(RP),
intent(out) :: cs,sn
188 else if (abs(b)>abs(a))
then 190 sn = 1._rp / sqrt( 1._rp + tmp**2)
195 cs = 1._rp / sqrt( 1._rp + tmp**2)
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
subroutine, public gmres(x, iter, nbPrd, res, b, A, tol, itmax, rst, verb)
GMRES (no preconditioning)
conversion integers or rational to real
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
subroutine, public invtrisup(res, M, vec)
Calcul de x=M^{-1}vec avec M triangulaire superieure.
real(rp) function e(x, v1, v2)
x = x + b*y // OR // z = x + b*y
subroutine grotmat(cs, sn, a, b)
Matrice de rotation de Givens.
real(rp), parameter, public real_tol