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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
!
!-----------------------------------------------------------------------
!
!Purpose:
! GMRES(m) linear equation solver.
!
!Copyright:
! Copyright (c) 2019 Ecole Polytechnique de Montreal
! this library is free software; you can redistribute it and/or
! modify it under the terms of the gnu lesser general public
! license as published by the free software foundation; either
! version 2.1 of the license, or (at your option) any later version
!
!Author(s): A. Hebert
!
!Reference:
! based on a Matlab script by C. T. Kelley, July 10, 1994.
!
!Parameters: input
! B fixed source
! atv function pointer for the matrix-vector product returning
! X+M*(B-A*X) where X is the unknown and B is the source.
! The format for atv is "Y=atv(X,B,n,...)"
! n order of matrix A
! ertol iteration convergence criterion
! nstart restarts the GMRES method every nstart iterations
! maxit maximum number of GMRES iterations.
! impx print parameter: =0: no print; =1: minimum printing.
! iptrk L_TRACK pointer to the tracking information
! ipsys L_SYSTEM pointer to system matrices
! ipflux L_FLUX pointer to the solution
!
!Parameters: input/output
! X initial estimate / solution of the linear system.
!
!Parameters: output
! iter actual number of iterations
!
!----------------------------------------------------------------------------
!
subroutine FLDMRA(B,atv,n,ertol,nstart,maxit,impx,iptrk,ipsys,ipflux,X,iter)
use GANLIB
implicit real(kind=8) (a-h,o-z)
!----
! subroutine arguments
!----
real(kind=8), dimension(n), intent(in) :: B
integer, intent(in) :: nstart,maxit,impx
real(kind=8), intent(in) :: ertol
interface
function atv(X,B,n,iptrk,ipsys,ipflux) result(Y)
use GANLIB
integer, intent(in) :: n
real(kind=8), dimension(n), intent(in) :: X, B
real(kind=8), dimension(n) :: Y
type(c_ptr) iptrk,ipsys,ipflux
end function atv
end interface
real(kind=8), dimension(n), intent(inout) :: X
integer, intent(out) :: iter
type(c_ptr) iptrk,ipsys,ipflux
!----
! local variables
!----
integer, parameter :: iunout=6
!----
! allocatable arays
!----
real(kind=8), allocatable, dimension(:) :: r,qq,g,c,s
real(kind=8), allocatable, dimension(:,:) :: v,h
!----
! scratch storage allocation
!----
allocate(v(n,nstart+1),g(nstart+1),h(nstart+1,nstart+1), &
c(nstart+1),s(nstart+1))
!----
! global GMRES(m) iteration.
!----
allocate(r(n),qq(n))
eps1=ertol*sqrt(dot_product(B(:n),B(:n)))
rho=1.0d10
iter=0
do while((rho > eps1).and.(iter < maxit))
r(:)=atv(X,B,n,iptrk,ipsys,ipflux)-X(:)
rho=sqrt(dot_product(r(:n),r(:n)))
!----
! test for termination on entry
!----
if(rho < eps1) then
deallocate(qq,r)
go to 100
endif
!
g(:nstart+1)=0.0d0
h(:nstart,:nstart)=0.0d0
v(:n,:nstart+1)=0.0d0
c(:nstart+1)=0.0d0
s(:nstart+1)=0.0d0
g(1)=rho
v(:n,1)=r(:n)/rho
!----
! gmres(1) iteration
!----
k=0
do while((rho > eps1).and.(k < nstart).and.(iter < maxit))
k=k+1
iter=iter+1
if(impx > 2) write(iunout,200) iter,rho,eps1
qq(:n)=0.0d0
r(:)=atv(v(:,k),qq,n,iptrk,ipsys,ipflux)
v(:n,k+1)=v(:n,k)-r(:n)
!----
! modified Gram-Schmidt
!----
do j=1,k
hr=dot_product(v(:n,j),v(:n,k+1))
h(j,k)=hr
v(:n,k+1)=v(:n,k+1)-hr*v(:n,j)
enddo
h(k+1,k)=sqrt(dot_product(v(:n,k+1),v(:n,k+1)))
!----
! reorthogonalize
!----
do j=1,k
hr=dot_product(v(:n,j),v(:n,k+1))
h(j,k)=h(j,k)+hr
v(:n,k+1)=v(:n,k+1)-hr*v(:n,j)
enddo
h(k+1,k)=sqrt(dot_product(v(:n,k+1),v(:n,k+1)))
!----
! watch out for happy breakdown
!----
if(h(k+1,k) /= 0.0) then
v(:n,k+1)=v(:n,k+1)/h(k+1,k)
endif
!----
! form and store the information for the new Givens rotation
!----
do i=1,k-1
w1=c(i)*h(i,k)-s(i)*h(i+1,k)
w2=s(i)*h(i,k)+c(i)*h(i+1,k)
h(i,k)=w1
h(i+1,k)=w2
enddo
znu=sqrt(h(k,k)**2+h(k+1,k)**2)
if(znu /= 0.0) then
c(k)=h(k,k)/znu
s(k)=-h(k+1,k)/znu
h(k,k)=c(k)*h(k,k)-s(k)*h(k+1,k)
h(k+1,k)=0.0d0
w1=c(k)*g(k)-s(k)*g(k+1)
w2=s(k)*g(k)+c(k)*g(k+1)
g(k)=w1
g(k+1)=w2
endif
!----
! update the residual norm
!----
rho=abs(g(k+1))
enddo
!----
! at this point either k > nstart or rho < eps1.
! it's time to compute x and cycle.
!----
h(:k,k+1)=g(:k)
call ALSBD(k,1,h,ier,nstart+1)
if(ier /= 0) call XABORT('FLDMRA: singular matrix.')
do i=1,n
X(i)=X(i)+dot_product(v(i,:k),h(:k,k+1))
enddo
enddo
deallocate(qq,r)
!----
! scratch storage deallocation
!----
100 deallocate(s,c,h,g,v)
return
!
200 format(24h FLDMRA: outer iteration,i4,10h L2 norm=,1p,e11.4, &
6h eps1=,e11.4)
end subroutine FLDMRA
|
