1: *> \brief \b DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLASQ5 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq5.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq5.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq5.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
22: * DNM1, DNM2, IEEE, EPS )
23: *
24: * .. Scalar Arguments ..
25: * LOGICAL IEEE
26: * INTEGER I0, N0, PP
27: * DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS
28: * ..
29: * .. Array Arguments ..
30: * DOUBLE PRECISION Z( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DLASQ5 computes one dqds transform in ping-pong form, one
40: *> version for IEEE machines another for non IEEE machines.
41: *> \endverbatim
42: *
43: * Arguments:
44: * ==========
45: *
46: *> \param[in] I0
47: *> \verbatim
48: *> I0 is INTEGER
49: *> First index.
50: *> \endverbatim
51: *>
52: *> \param[in] N0
53: *> \verbatim
54: *> N0 is INTEGER
55: *> Last index.
56: *> \endverbatim
57: *>
58: *> \param[in] Z
59: *> \verbatim
60: *> Z is DOUBLE PRECISION array, dimension ( 4*N )
61: *> Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
62: *> an extra argument.
63: *> \endverbatim
64: *>
65: *> \param[in] PP
66: *> \verbatim
67: *> PP is INTEGER
68: *> PP=0 for ping, PP=1 for pong.
69: *> \endverbatim
70: *>
71: *> \param[in] TAU
72: *> \verbatim
73: *> TAU is DOUBLE PRECISION
74: *> This is the shift.
75: *> \endverbatim
76: *>
77: *> \param[in] SIGMA
78: *> \verbatim
79: *> SIGMA is DOUBLE PRECISION
80: *> This is the accumulated shift up to this step.
81: *> \endverbatim
82: *>
83: *> \param[out] DMIN
84: *> \verbatim
85: *> DMIN is DOUBLE PRECISION
86: *> Minimum value of d.
87: *> \endverbatim
88: *>
89: *> \param[out] DMIN1
90: *> \verbatim
91: *> DMIN1 is DOUBLE PRECISION
92: *> Minimum value of d, excluding D( N0 ).
93: *> \endverbatim
94: *>
95: *> \param[out] DMIN2
96: *> \verbatim
97: *> DMIN2 is DOUBLE PRECISION
98: *> Minimum value of d, excluding D( N0 ) and D( N0-1 ).
99: *> \endverbatim
100: *>
101: *> \param[out] DN
102: *> \verbatim
103: *> DN is DOUBLE PRECISION
104: *> d(N0), the last value of d.
105: *> \endverbatim
106: *>
107: *> \param[out] DNM1
108: *> \verbatim
109: *> DNM1 is DOUBLE PRECISION
110: *> d(N0-1).
111: *> \endverbatim
112: *>
113: *> \param[out] DNM2
114: *> \verbatim
115: *> DNM2 is DOUBLE PRECISION
116: *> d(N0-2).
117: *> \endverbatim
118: *>
119: *> \param[in] IEEE
120: *> \verbatim
121: *> IEEE is LOGICAL
122: *> Flag for IEEE or non IEEE arithmetic.
123: *> \endverbatim
124: *>
125: *> \param[in] EPS
126: *> \verbatim
127: *> EPS is DOUBLE PRECISION
128: *> This is the value of epsilon used.
129: *> \endverbatim
130: *>
131: * Authors:
132: * ========
133: *
134: *> \author Univ. of Tennessee
135: *> \author Univ. of California Berkeley
136: *> \author Univ. of Colorado Denver
137: *> \author NAG Ltd.
138: *
139: *> \ingroup auxOTHERcomputational
140: *
141: * =====================================================================
142: SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
143: $ DN, DNM1, DNM2, IEEE, EPS )
144: *
145: * -- LAPACK computational routine --
146: * -- LAPACK is a software package provided by Univ. of Tennessee, --
147: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148: *
149: * .. Scalar Arguments ..
150: LOGICAL IEEE
151: INTEGER I0, N0, PP
152: DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
153: $ SIGMA, EPS
154: * ..
155: * .. Array Arguments ..
156: DOUBLE PRECISION Z( * )
157: * ..
158: *
159: * =====================================================================
160: *
161: * .. Parameter ..
162: DOUBLE PRECISION ZERO, HALF
163: PARAMETER ( ZERO = 0.0D0, HALF = 0.5 )
164: * ..
165: * .. Local Scalars ..
166: INTEGER J4, J4P2
167: DOUBLE PRECISION D, EMIN, TEMP, DTHRESH
168: * ..
169: * .. Intrinsic Functions ..
170: INTRINSIC MIN
171: * ..
172: * .. Executable Statements ..
173: *
174: IF( ( N0-I0-1 ).LE.0 )
175: $ RETURN
176: *
177: DTHRESH = EPS*(SIGMA+TAU)
178: IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
179: IF( TAU.NE.ZERO ) THEN
180: J4 = 4*I0 + PP - 3
181: EMIN = Z( J4+4 )
182: D = Z( J4 ) - TAU
183: DMIN = D
184: DMIN1 = -Z( J4 )
185: *
186: IF( IEEE ) THEN
187: *
188: * Code for IEEE arithmetic.
189: *
190: IF( PP.EQ.0 ) THEN
191: DO 10 J4 = 4*I0, 4*( N0-3 ), 4
192: Z( J4-2 ) = D + Z( J4-1 )
193: TEMP = Z( J4+1 ) / Z( J4-2 )
194: D = D*TEMP - TAU
195: DMIN = MIN( DMIN, D )
196: Z( J4 ) = Z( J4-1 )*TEMP
197: EMIN = MIN( Z( J4 ), EMIN )
198: 10 CONTINUE
199: ELSE
200: DO 20 J4 = 4*I0, 4*( N0-3 ), 4
201: Z( J4-3 ) = D + Z( J4 )
202: TEMP = Z( J4+2 ) / Z( J4-3 )
203: D = D*TEMP - TAU
204: DMIN = MIN( DMIN, D )
205: Z( J4-1 ) = Z( J4 )*TEMP
206: EMIN = MIN( Z( J4-1 ), EMIN )
207: 20 CONTINUE
208: END IF
209: *
210: * Unroll last two steps.
211: *
212: DNM2 = D
213: DMIN2 = DMIN
214: J4 = 4*( N0-2 ) - PP
215: J4P2 = J4 + 2*PP - 1
216: Z( J4-2 ) = DNM2 + Z( J4P2 )
217: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
218: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
219: DMIN = MIN( DMIN, DNM1 )
220: *
221: DMIN1 = DMIN
222: J4 = J4 + 4
223: J4P2 = J4 + 2*PP - 1
224: Z( J4-2 ) = DNM1 + Z( J4P2 )
225: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
226: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
227: DMIN = MIN( DMIN, DN )
228: *
229: ELSE
230: *
231: * Code for non IEEE arithmetic.
232: *
233: IF( PP.EQ.0 ) THEN
234: DO 30 J4 = 4*I0, 4*( N0-3 ), 4
235: Z( J4-2 ) = D + Z( J4-1 )
236: IF( D.LT.ZERO ) THEN
237: RETURN
238: ELSE
239: Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
240: D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
241: END IF
242: DMIN = MIN( DMIN, D )
243: EMIN = MIN( EMIN, Z( J4 ) )
244: 30 CONTINUE
245: ELSE
246: DO 40 J4 = 4*I0, 4*( N0-3 ), 4
247: Z( J4-3 ) = D + Z( J4 )
248: IF( D.LT.ZERO ) THEN
249: RETURN
250: ELSE
251: Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
252: D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
253: END IF
254: DMIN = MIN( DMIN, D )
255: EMIN = MIN( EMIN, Z( J4-1 ) )
256: 40 CONTINUE
257: END IF
258: *
259: * Unroll last two steps.
260: *
261: DNM2 = D
262: DMIN2 = DMIN
263: J4 = 4*( N0-2 ) - PP
264: J4P2 = J4 + 2*PP - 1
265: Z( J4-2 ) = DNM2 + Z( J4P2 )
266: IF( DNM2.LT.ZERO ) THEN
267: RETURN
268: ELSE
269: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
270: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
271: END IF
272: DMIN = MIN( DMIN, DNM1 )
273: *
274: DMIN1 = DMIN
275: J4 = J4 + 4
276: J4P2 = J4 + 2*PP - 1
277: Z( J4-2 ) = DNM1 + Z( J4P2 )
278: IF( DNM1.LT.ZERO ) THEN
279: RETURN
280: ELSE
281: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
282: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
283: END IF
284: DMIN = MIN( DMIN, DN )
285: *
286: END IF
287: ELSE
288: * This is the version that sets d's to zero if they are small enough
289: J4 = 4*I0 + PP - 3
290: EMIN = Z( J4+4 )
291: D = Z( J4 ) - TAU
292: DMIN = D
293: DMIN1 = -Z( J4 )
294: IF( IEEE ) THEN
295: *
296: * Code for IEEE arithmetic.
297: *
298: IF( PP.EQ.0 ) THEN
299: DO 50 J4 = 4*I0, 4*( N0-3 ), 4
300: Z( J4-2 ) = D + Z( J4-1 )
301: TEMP = Z( J4+1 ) / Z( J4-2 )
302: D = D*TEMP - TAU
303: IF( D.LT.DTHRESH ) D = ZERO
304: DMIN = MIN( DMIN, D )
305: Z( J4 ) = Z( J4-1 )*TEMP
306: EMIN = MIN( Z( J4 ), EMIN )
307: 50 CONTINUE
308: ELSE
309: DO 60 J4 = 4*I0, 4*( N0-3 ), 4
310: Z( J4-3 ) = D + Z( J4 )
311: TEMP = Z( J4+2 ) / Z( J4-3 )
312: D = D*TEMP - TAU
313: IF( D.LT.DTHRESH ) D = ZERO
314: DMIN = MIN( DMIN, D )
315: Z( J4-1 ) = Z( J4 )*TEMP
316: EMIN = MIN( Z( J4-1 ), EMIN )
317: 60 CONTINUE
318: END IF
319: *
320: * Unroll last two steps.
321: *
322: DNM2 = D
323: DMIN2 = DMIN
324: J4 = 4*( N0-2 ) - PP
325: J4P2 = J4 + 2*PP - 1
326: Z( J4-2 ) = DNM2 + Z( J4P2 )
327: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
328: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
329: DMIN = MIN( DMIN, DNM1 )
330: *
331: DMIN1 = DMIN
332: J4 = J4 + 4
333: J4P2 = J4 + 2*PP - 1
334: Z( J4-2 ) = DNM1 + Z( J4P2 )
335: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
336: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
337: DMIN = MIN( DMIN, DN )
338: *
339: ELSE
340: *
341: * Code for non IEEE arithmetic.
342: *
343: IF( PP.EQ.0 ) THEN
344: DO 70 J4 = 4*I0, 4*( N0-3 ), 4
345: Z( J4-2 ) = D + Z( J4-1 )
346: IF( D.LT.ZERO ) THEN
347: RETURN
348: ELSE
349: Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
350: D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
351: END IF
352: IF( D.LT.DTHRESH) D = ZERO
353: DMIN = MIN( DMIN, D )
354: EMIN = MIN( EMIN, Z( J4 ) )
355: 70 CONTINUE
356: ELSE
357: DO 80 J4 = 4*I0, 4*( N0-3 ), 4
358: Z( J4-3 ) = D + Z( J4 )
359: IF( D.LT.ZERO ) THEN
360: RETURN
361: ELSE
362: Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
363: D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
364: END IF
365: IF( D.LT.DTHRESH) D = ZERO
366: DMIN = MIN( DMIN, D )
367: EMIN = MIN( EMIN, Z( J4-1 ) )
368: 80 CONTINUE
369: END IF
370: *
371: * Unroll last two steps.
372: *
373: DNM2 = D
374: DMIN2 = DMIN
375: J4 = 4*( N0-2 ) - PP
376: J4P2 = J4 + 2*PP - 1
377: Z( J4-2 ) = DNM2 + Z( J4P2 )
378: IF( DNM2.LT.ZERO ) THEN
379: RETURN
380: ELSE
381: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
382: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
383: END IF
384: DMIN = MIN( DMIN, DNM1 )
385: *
386: DMIN1 = DMIN
387: J4 = J4 + 4
388: J4P2 = J4 + 2*PP - 1
389: Z( J4-2 ) = DNM1 + Z( J4P2 )
390: IF( DNM1.LT.ZERO ) THEN
391: RETURN
392: ELSE
393: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
394: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
395: END IF
396: DMIN = MIN( DMIN, DN )
397: *
398: END IF
399: END IF
400: *
401: Z( J4+2 ) = DN
402: Z( 4*N0-PP ) = EMIN
403: RETURN
404: *
405: * End of DLASQ5
406: *
407: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq5.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack