KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlasq5.f
Go to the documentation of this file.
1  SUBROUTINE dlasq5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
2  $ DNM1, DNM2, IEEE )
3 *
4 * -- LAPACK auxiliary routine (version 3.0) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * May 17, 2000
8 *
9 * .. Scalar Arguments ..
10  LOGICAL IEEE
11  INTEGER I0, N0, PP
12  DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU
13 * ..
14 * .. Array Arguments ..
15  DOUBLE PRECISION Z( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLASQ5 computes one dqds transform in ping-pong form, one
22 * version for IEEE machines another for non IEEE machines.
23 *
24 * Arguments
25 * =========
26 *
27 * I0 (input) INTEGER
28 * First index.
29 *
30 * N0 (input) INTEGER
31 * Last index.
32 *
33 * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
34 * Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
35 * an extra argument.
36 *
37 * PP (input) INTEGER
38 * PP=0 for ping, PP=1 for pong.
39 *
40 * TAU (input) DOUBLE PRECISION
41 * This is the shift.
42 *
43 * DMIN (output) DOUBLE PRECISION
44 * Minimum value of d.
45 *
46 * DMIN1 (output) DOUBLE PRECISION
47 * Minimum value of d, excluding D( N0 ).
48 *
49 * DMIN2 (output) DOUBLE PRECISION
50 * Minimum value of d, excluding D( N0 ) and D( N0-1 ).
51 *
52 * DN (output) DOUBLE PRECISION
53 * d(N0), the last value of d.
54 *
55 * DNM1 (output) DOUBLE PRECISION
56 * d(N0-1).
57 *
58 * DNM2 (output) DOUBLE PRECISION
59 * d(N0-2).
60 *
61 * IEEE (input) LOGICAL
62 * Flag for IEEE or non IEEE arithmetic.
63 *
64 * =====================================================================
65 *
66 * .. Parameter ..
67  DOUBLE PRECISION ZERO
68  parameter( zero = 0.0d0 )
69 * ..
70 * .. Local Scalars ..
71  INTEGER J4, J4P2
72  DOUBLE PRECISION D, EMIN, TEMP
73 * ..
74 * .. Intrinsic Functions ..
75  INTRINSIC min
76 * ..
77 * .. Executable Statements ..
78 *
79  IF( ( n0-i0-1 ).LE.0 )
80  $ RETURN
81 *
82  j4 = 4*i0 + pp - 3
83  emin = z( j4+4 )
84  d = z( j4 ) - tau
85  dmin = d
86  dmin1 = -z( j4 )
87 *
88  IF( ieee ) THEN
89 *
90 * Code for IEEE arithmetic.
91 *
92  IF( pp.EQ.0 ) THEN
93  DO 10 j4 = 4*i0, 4*( n0-3 ), 4
94  z( j4-2 ) = d + z( j4-1 )
95  temp = z( j4+1 ) / z( j4-2 )
96  d = d*temp - tau
97  dmin = min( dmin, d )
98  z( j4 ) = z( j4-1 )*temp
99  emin = min( z( j4 ), emin )
100  10 CONTINUE
101  ELSE
102  DO 20 j4 = 4*i0, 4*( n0-3 ), 4
103  z( j4-3 ) = d + z( j4 )
104  temp = z( j4+2 ) / z( j4-3 )
105  d = d*temp - tau
106  dmin = min( dmin, d )
107  z( j4-1 ) = z( j4 )*temp
108  emin = min( z( j4-1 ), emin )
109  20 CONTINUE
110  END IF
111 *
112 * Unroll last two steps.
113 *
114  dnm2 = d
115  dmin2 = dmin
116  j4 = 4*( n0-2 ) - pp
117  j4p2 = j4 + 2*pp - 1
118  z( j4-2 ) = dnm2 + z( j4p2 )
119  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
120  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
121  dmin = min( dmin, dnm1 )
122 *
123  dmin1 = dmin
124  j4 = j4 + 4
125  j4p2 = j4 + 2*pp - 1
126  z( j4-2 ) = dnm1 + z( j4p2 )
127  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
128  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
129  dmin = min( dmin, dn )
130 *
131  ELSE
132 *
133 * Code for non IEEE arithmetic.
134 *
135  IF( pp.EQ.0 ) THEN
136  DO 30 j4 = 4*i0, 4*( n0-3 ), 4
137  z( j4-2 ) = d + z( j4-1 )
138  IF( d.LT.zero ) THEN
139  RETURN
140  ELSE
141  z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
142  d = z( j4+1 )*( d / z( j4-2 ) ) - tau
143  END IF
144  dmin = min( dmin, d )
145  emin = min( emin, z( j4 ) )
146  30 CONTINUE
147  ELSE
148  DO 40 j4 = 4*i0, 4*( n0-3 ), 4
149  z( j4-3 ) = d + z( j4 )
150  IF( d.LT.zero ) THEN
151  RETURN
152  ELSE
153  z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
154  d = z( j4+2 )*( d / z( j4-3 ) ) - tau
155  END IF
156  dmin = min( dmin, d )
157  emin = min( emin, z( j4-1 ) )
158  40 CONTINUE
159  END IF
160 *
161 * Unroll last two steps.
162 *
163  dnm2 = d
164  dmin2 = dmin
165  j4 = 4*( n0-2 ) - pp
166  j4p2 = j4 + 2*pp - 1
167  z( j4-2 ) = dnm2 + z( j4p2 )
168  IF( dnm2.LT.zero ) THEN
169  RETURN
170  ELSE
171  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
172  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
173  END IF
174  dmin = min( dmin, dnm1 )
175 *
176  dmin1 = dmin
177  j4 = j4 + 4
178  j4p2 = j4 + 2*pp - 1
179  z( j4-2 ) = dnm1 + z( j4p2 )
180  IF( dnm1.LT.zero ) THEN
181  RETURN
182  ELSE
183  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
184  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
185  END IF
186  dmin = min( dmin, dn )
187 *
188  END IF
189 *
190  z( j4+2 ) = dn
191  z( 4*n0-pp ) = emin
192  RETURN
193 *
194 * End of DLASQ5
195 *
196  END
subroutine dlasq5(I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE)
Definition: dlasq5.f:3