20 SUBROUTINE dutfx (NDIM,NLEN,NLIM,VECS,TOL,INFO)
21 USE stel_kinds,
ONLY: rprec
71 EXTERNAL dlamch, daxpby, drandn, dnrm2, ddot
72 REAL(rprec) :: DLAMCH, DNRM2, DDOT
74 INTEGER :: INFO(4), NDIM, NLEN, NLIM
75 REAL(rprec) :: VECS(NDIM,9)
80 REAL(rprec) DHUN, DONE, DTEN, DZERO
81 parameter(dhun = 1.0d2,done = 1.0d0,dten = 1.0d1,dzero = 0.0d0)
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
111 IF (revcom.EQ.0)
THEN
113 IF (retlbl.EQ.0)
GO TO 10
114 ELSE IF (revcom.EQ.1)
THEN
115 IF (retlbl.EQ.30)
THEN
117 ELSE IF (retlbl.EQ.40)
THEN
119 ELSE IF (retlbl.EQ.60)
THEN
121 ELSE IF (retlbl.EQ.70)
THEN
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
142 vf = vf - init * 100000
144 vf = vf - tres * 10000
150 IF (tol.LE.dzero) tol = sqrt(dlamch(
'E'))
155 WRITE (vf, *)
' N 2N-1 UNRM RESN'
156 WRITE (vf,
'(2I8,1P,2E11.4)') 0, 0, done, done
159 WRITE (vf, *)
' N 2N-1 UNRM RESN'
160 WRITE (tf,
'(2I8,1P,2E11.4)') 0, 0, done, done
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
172 IF (init.EQ.0)
CALL drandn (nlen,vecs(1,3),1)
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))
190 20 dtmp = ddot(nlen,vecs(1,3),1,vecs(1,5),1)
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))
207 30
CALL daxpby (nlen,vecs(1,4),beta,vecs(1,4),done,vecs(1,9))
211 dtmp = ddot(nlen,vecs(1,3),1,vecs(1,4),1)
212 IF ((dtmp.EQ.dzero).OR.(rho.EQ.dzero))
THEN
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))
227 dtmp = dnrm2(nlen,vecs(1,5),1)
230 cos1 = done / ( done + var )
236 CALL daxpby (nlen,vecs(1,1),done,vecs(1,1),eta,vecs(1,7))
242 unrm = sqrt((2*n) * tau) / r0
244 IF ((tres.EQ.0).AND.(unrm/tol.GT.dten))
GO TO 50
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
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
268 IF (resn.LE.tol)
THEN
271 ELSE IF (unrm.LT.uchk/dhun)
THEN
278 CALL daxpby (nlen,vecs(1,6),done,vecs(1,6),-alpha,vecs(1,4))
280 CALL daxpby (nlen,vecs(1,7),done,vecs(1,6),dtmp,vecs(1,7))
290 60
CALL daxpby (nlen,vecs(1,5),done,vecs(1,5),-alpha,vecs(1,8))
295 dtmp = dnrm2(nlen,vecs(1,5),1)
298 cos1 = done / ( done + var )
304 CALL daxpby (nlen,vecs(1,1),done,vecs(1,1),eta,vecs(1,7))
310 unrm = sqrt((2*n+1) * tau) / r0
312 IF ((tres.EQ.0).AND.(unrm/tol.GT.dten).AND.(n.LT.nlim))
GO TO 80
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
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
337 IF (resn.LE.tol)
THEN
340 ELSE IF (unrm.LT.uchk/dhun)
THEN
343 ELSE IF (n.GE.nlim)
THEN