V3FIT
dutfx.f
1 C**********************************************************************
2 C
3 C Copyright (C) 1992 Roland W. Freund and Noel M. Nachtigal
4 C All rights reserved.
5 C
6 C This code is part of a copyrighted package. For details, see the
7 C file "cpyrit.doc" in the top-level directory.
8 C
9 C *****************************************************************
10 C ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE
11 C COPYRIGHT NOTICE
12 C *****************************************************************
13 C
14 C**********************************************************************
15 C
16 C This file contains the routine for the TFQMR algorithm.
17 C
18 C**********************************************************************
19 C
20  SUBROUTINE dutfx (NDIM,NLEN,NLIM,VECS,TOL,INFO)
21  USE stel_kinds, ONLY: rprec
22  IMPLICIT NONE
23 C
24 C Purpose:
25 C This subroutine uses the TFQMR (Transpose Free, QMR) algorithm to
26 C solve linear systems.
27 C It runs the algorithm to convergence or until a user-specified
28 C limit on the number of iterations is reached.
29 C
30 C The code is set up to solve the system A x = b with initial
31 C guess x_0 = 0. Here A x = b denotes the preconditioned system,
32 C and it is connected with the original system as follows. Let
33 C B y = c be the original unpreconditioned system to be solved, and
34 C let y_0 be an arbitrary initial guess for its solution. Then:
35 C A x = b, where A = M_1^{-1} B M_2^{-1},
36 C x = M_2 (y - y_0), b = M_1^{-1} (c - B y_0).
37 C Here M = M_1 M_2 is the preconditioner.
38 C
39 C To recover the final iterate y_n for the original system B y = c
40 C from the final iterate x_n for the preconditioned system A x = b,
41 C set
42 C y_n = y_0 + M_2^{-1} x_n.
43 C
44 C The algorithm was first described in the RIACS Technical Report
45 C 91.18, `A Transpose-Free Quasi-Minimal Residual Algorithm for
46 C Non-Hermitian Linear Systems`, by Roland Freund, September 1991,
47 C which subsequently appeared in SIAM J. Sci. Comput., 14 (1993),
48 C pp. 470--482.
49 C
50 C Parameters:
51 C For a description of the parameters, see the file "dutfx.doc" in
52 C the current directory.
53 C
54 C External routines used:
55 C double precision dlamch(ch)
56 C LAPACK routine, computes machine-related constants.
57 C double precision dnrm2(n,x,incx)
58 C BLAS-1 routine, computes the 2-norm of x.
59 C subroutine daxpby(n,z,a,x,b,y)
60 C Library routine, computes z = a * x + b * y.
61 C double precision ddot(n,x,incx,y,incy)
62 C BLAS-1 routine, computes y^H * x.
63 C subroutine drandn(n,x,seed)
64 C Library routine, fills x with random numbers.
65 C
66 C Noel M. Nachtigal
67 C April 13, 1993
68 C
69 C**********************************************************************
70 C
71  EXTERNAL dlamch, daxpby, drandn, dnrm2, ddot
72  REAL(rprec) :: DLAMCH, DNRM2, DDOT
73 C
74  INTEGER :: INFO(4), NDIM, NLEN, NLIM
75  REAL(rprec) :: VECS(NDIM,9)
76  REAL(rprec) :: TOL
77 C
78 C Miscellaneous parameters.
79 C
80  REAL(rprec) DHUN, DONE, DTEN, DZERO
81  parameter(dhun = 1.0d2,done = 1.0d0,dten = 1.0d1,dzero = 0.0d0)
82 C
83 C Local variables, permanent.
84 C
85  INTEGER :: IERR, N, RETLBL=0, tf, tres, vf
86  SAVE :: ierr, n, retlbl, tf, tres, vf
87  REAL(rprec) :: ALPHA, BETA, ETA, RHO
88  SAVE :: alpha, beta, eta, rho
89  REAL(rprec) :: COS1, VAR, R0, RESN, TAU, UCHK, UNRM
90  SAVE :: cos1, var, r0, resn, tau, uchk, unrm
91 C
92 C Local variables, transient.
93 C
94  INTEGER INIT, REVCOM
95  REAL(rprec) DTMP
96 C
97 C Initialize some of the permanent variables.
98 C
99 C DATA RETLBL /0/
100 C
101 C Check the reverse communication flag to see where to branch.
102 C REVCOM RETLBL Comment
103 C 0 0 first call, go to label 10
104 C 1 30 returning from AXB, go to label 30
105 C 1 40 returning from AXB, go to label 40
106 C 1 60 returning from AXB, go to label 60
107 C 1 70 returning from AXB, go to label 70
108 C
109  revcom = info(2)
110  info(2) = 0
111  IF (revcom.EQ.0) THEN
112  n = 0
113  IF (retlbl.EQ.0) GO TO 10
114  ELSE IF (revcom.EQ.1) THEN
115  IF (retlbl.EQ.30) THEN
116  GO TO 30
117  ELSE IF (retlbl.EQ.40) THEN
118  GO TO 40
119  ELSE IF (retlbl.EQ.60) THEN
120  GO TO 60
121  ELSE IF (retlbl.EQ.70) THEN
122  GO TO 70
123  END IF
124  END IF
125  ierr = 1
126  GO TO 90
127 C
128 C Check whether the inputs are valid.
129 C
130  10 ierr = 0
131  IF (ndim.LT.1) ierr = 2
132  IF (nlen.LT.1) ierr = 2
133  IF (nlim.LT.1) ierr = 2
134  IF (nlen.GT.ndim) ierr = 2
135  IF (ierr.NE.0) GO TO 90
136 C
137 C Extract from INFO the output units TF and VF, the true residual
138 C flag TRES, and the left starting vector flag INIT.
139 C
140  vf = max(info(1),0)
141  init = vf / 100000
142  vf = vf - init * 100000
143  tres = vf / 10000
144  vf = vf - tres * 10000
145  tf = vf / 100
146  vf = vf - tf * 100
147 C
148 C Check the convergence tolerance.
149 C
150  IF (tol.LE.dzero) tol = sqrt(dlamch('E'))
151 C
152 C Start the trace messages and convergence history.
153 C
154  IF (vf.NE.0) THEN
155  WRITE (vf, *) ' N 2N-1 UNRM RESN'
156  WRITE (vf,'(2I8,1P,2E11.4)') 0, 0, done, done
157  ENDIF
158  IF (tf.NE.0) THEN
159  WRITE (vf, *) ' N 2N-1 UNRM RESN'
160  WRITE (tf,'(2I8,1P,2E11.4)') 0, 0, done, done
161  ENDIF
162 C
163 C Set x_0 = 0 and compute the norm of the initial residual.
164 C
165  CALL daxpby (nlen,vecs(1,5),done,vecs(1,2),dzero,vecs(1,5))
166  CALL daxpby (nlen,vecs(1,1),dzero,vecs(1,1),dzero,vecs(1,1))
167  r0 = dnrm2(nlen,vecs(1,5),1)
168  IF ((tol.GE.done).OR.(r0.EQ.dzero)) GO TO 90
169 C
170 C Check whether the auxiliary vector must be supplied.
171 C
172  IF (init.EQ.0) CALL drandn (nlen,vecs(1,3),1)
173 C
174 C Initialize the variables.
175 C
176  n = 1
177  resn = done
178  rho = done
179  var = dzero
180  eta = dzero
181  tau = r0 * r0
182  ierr = 8
183  CALL daxpby (nlen,vecs(1,8),dzero,vecs(1,8),dzero,vecs(1,8))
184  CALL daxpby (nlen,vecs(1,4),dzero,vecs(1,4),dzero,vecs(1,4))
185  CALL daxpby (nlen,vecs(1,6),dzero,vecs(1,6),dzero,vecs(1,6))
186 C
187 C This is one step of the TFQMR algorithm.
188 C Compute \beta_{n-1} and \rho_{n-1}.
189 C
190  20 dtmp = ddot(nlen,vecs(1,3),1,vecs(1,5),1)
191  beta = dtmp / rho
192  rho = dtmp
193 C
194 C Compute y_{2n-1}, v_{n-1}, and A y_{2n-1}.
195 C
196  CALL daxpby (nlen,vecs(1,4),beta,vecs(1,4),done,vecs(1,8))
197  CALL daxpby (nlen,vecs(1,6),done,vecs(1,5),beta,vecs(1,6))
198 C
199 C Have the caller carry out AXB, then return here.
200 C CALL AXB (VECS(1,6),VECS(1,9))
201 C
202  info(2) = 1
203  info(3) = 6
204  info(4) = 9
205  retlbl = 30
206  RETURN
207  30 CALL daxpby (nlen,vecs(1,4),beta,vecs(1,4),done,vecs(1,9))
208 C
209 C Compute \sigma{n-1} and check for breakdowns.
210 C
211  dtmp = ddot(nlen,vecs(1,3),1,vecs(1,4),1)
212  IF ((dtmp.EQ.dzero).OR.(rho.EQ.dzero)) THEN
213  ierr = 8
214  GO TO 90
215  END IF
216 C
217 C Compute \alpha_{n-1}, d_{2n-1} and w_{2n}.
218 C
219  alpha = rho / dtmp
220  dtmp = var * eta / alpha
221  CALL daxpby (nlen,vecs(1,7),done,vecs(1,6),dtmp,vecs(1,7))
222  CALL daxpby (nlen,vecs(1,5),done,vecs(1,5),-alpha,vecs(1,9))
223 C
224 C Compute \varepsilon_{2n-1}^2, \eta_{2n-1}^2, c_{2n-1}^2, and
225 C \tau_{2n-1}^2.
226 C
227  dtmp = dnrm2(nlen,vecs(1,5),1)
228  dtmp = dtmp * dtmp
229  var = dtmp / tau
230  cos1 = done / ( done + var )
231  tau = dtmp * cos1
232  eta = alpha * cos1
233 C
234 C Compute x_{2n-1} and the upper bound for its residual norm.
235 C
236  CALL daxpby (nlen,vecs(1,1),done,vecs(1,1),eta,vecs(1,7))
237 C
238 C Compute the residual norm upper bound.
239 C If the scaled upper bound is within one order of magnitude of the
240 C target convergence norm, compute the true residual norm.
241 C
242  unrm = sqrt((2*n) * tau) / r0
243  uchk = unrm
244  IF ((tres.EQ.0).AND.(unrm/tol.GT.dten)) GO TO 50
245 C
246 C Have the caller carry out AXB, then return here.
247 C CALL AXB (VECS(1,1),VECS(1,9))
248 C
249  info(2) = 1
250  info(3) = 1
251  info(4) = 9
252  retlbl = 40
253  RETURN
254  40 CALL daxpby (nlen,vecs(1,9),done,vecs(1,2),-done,vecs(1,9))
255  resn = dnrm2(nlen,vecs(1,9),1) / r0
256  uchk = resn
257 C
258 C Output the trace messages and convergence history.
259 C
260  50 IF (vf.NE.0) WRITE (vf,'(2I8,2E11.4)') n, 2*n-1, unrm, resn
261  IF (tf.NE.0) WRITE (tf,'(2I8,2E11.4)') n, 2*n-1, unrm, resn
262 C
263 C Check for convergence or termination. Stop if:
264 C 1. algorithm converged;
265 C 2. the residual norm upper bound is smaller than the computed
266 C residual norm by a factor of at least 100.
267 C
268  IF (resn.LE.tol) THEN
269  ierr = 0
270  GO TO 90
271  ELSE IF (unrm.LT.uchk/dhun) THEN
272  ierr = 4
273  GO TO 90
274  END IF
275 C
276 C Compute y_{2n}, A y_{2n}, d_{2n}, and w_{2n+1}.
277 C
278  CALL daxpby (nlen,vecs(1,6),done,vecs(1,6),-alpha,vecs(1,4))
279  dtmp = var * cos1
280  CALL daxpby (nlen,vecs(1,7),done,vecs(1,6),dtmp,vecs(1,7))
281 C
282 C Have the caller carry out AXB, then return here.
283 C CALL AXB (VECS(1,6),VECS(1,8))
284 C
285  info(2) = 1
286  info(3) = 6
287  info(4) = 8
288  retlbl = 60
289  RETURN
290  60 CALL daxpby (nlen,vecs(1,5),done,vecs(1,5),-alpha,vecs(1,8))
291 C
292 C Compute \varepsilon_{2n}^2, \eta_{2n}^2, c_{2n}^2, and
293 C \tau_{2n}^2.
294 C
295  dtmp = dnrm2(nlen,vecs(1,5),1)
296  dtmp = dtmp * dtmp
297  var = dtmp / tau
298  cos1 = done / ( done + var )
299  tau = dtmp * cos1
300  eta = alpha * cos1
301 C
302 C Compute x_{2n}.
303 C
304  CALL daxpby (nlen,vecs(1,1),done,vecs(1,1),eta,vecs(1,7))
305 C
306 C Compute the residual norm upper bound.
307 C If the scaled upper bound is within one order of magnitude of the
308 C target convergence norm, compute the true residual norm.
309 C
310  unrm = sqrt((2*n+1) * tau) / r0
311  uchk = unrm
312  IF ((tres.EQ.0).AND.(unrm/tol.GT.dten).AND.(n.LT.nlim)) GO TO 80
313 C
314 C Have the caller carry out AXB, then return here.
315 C CALL AXB (VECS(1,1),VECS(1,9))
316 C
317  info(2) = 1
318  info(3) = 1
319  info(4) = 9
320  retlbl = 70
321  RETURN
322  70 CALL daxpby (nlen,vecs(1,9),done,vecs(1,2),-done,vecs(1,9))
323  resn = dnrm2(nlen,vecs(1,9),1) / r0
324  uchk = unrm
325 C
326 C Output the trace messages and convergence history.
327 C
328  80 IF (vf.NE.0) WRITE (vf,'(2I8,2E11.4)') n, 2*n, unrm, resn
329  IF (tf.NE.0) WRITE (tf,'(2I8,2E11.4)') n, 2*n, unrm, resn
330 C
331 C Check for convergence or termination. Stop if:
332 C 1. algorithm converged;
333 C 2. the residual norm upper bound is smaller than the computed
334 C residual norm by a factor of at least 100;
335 C 3. algorithm exceeded the iterations limit.
336 C
337  IF (resn.LE.tol) THEN
338  ierr = 0
339  GO TO 90
340  ELSE IF (unrm.LT.uchk/dhun) THEN
341  ierr = 4
342  GO TO 90
343  ELSE IF (n.GE.nlim) THEN
344  ierr = 4
345  GO TO 90
346  END IF
347 C
348 C Update the running counter.
349 C
350  n = n + 1
351  GO TO 20
352 C
353 C Done.
354 C
355  90 nlim = n
356  retlbl = 0
357  info(1) = ierr
358 C
359  RETURN
360  END
361 C
362 C**********************************************************************
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17