![[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|
Return to dlasq5.f CVS log | Up to [local] / rpl / lapack / lapack |
1.12 bertrand 1: *> \brief \b DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
1.8 bertrand 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: *
1.10 bertrand 21: * SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
22: * DNM1, DNM2, IEEE, EPS )
1.8 bertrand 23: *
24: * .. Scalar Arguments ..
25: * LOGICAL IEEE
26: * INTEGER I0, N0, PP
1.10 bertrand 27: * DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS
1.8 bertrand 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: *>
1.10 bertrand 77: *> \param[in] SIGMA
78: *> \verbatim
79: *> SIGMA is DOUBLE PRECISION
80: *> This is the accumulated shift up to this step.
81: *> \endverbatim
82: *>
1.8 bertrand 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: *
1.10 bertrand 125: *> \param[in] EPS
126: *> \verbatim
127: *> EPS is DOUBLE PRECISION
128: *> This is the value of epsilon used.
129: *> \endverbatim
130: *>
1.8 bertrand 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: *
1.12 bertrand 139: *> \date September 2012
1.8 bertrand 140: *
141: *> \ingroup auxOTHERcomputational
142: *
143: * =====================================================================
1.10 bertrand 144: SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
145: $ DN, DNM1, DNM2, IEEE, EPS )
1.1 bertrand 146: *
1.12 bertrand 147: * -- LAPACK computational routine (version 3.4.2) --
1.1 bertrand 148: * -- LAPACK is a software package provided by Univ. of Tennessee, --
149: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.12 bertrand 150: * September 2012
1.1 bertrand 151: *
152: * .. Scalar Arguments ..
153: LOGICAL IEEE
154: INTEGER I0, N0, PP
1.10 bertrand 155: DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
156: $ SIGMA, EPS
1.1 bertrand 157: * ..
158: * .. Array Arguments ..
159: DOUBLE PRECISION Z( * )
160: * ..
161: *
162: * =====================================================================
163: *
164: * .. Parameter ..
1.10 bertrand 165: DOUBLE PRECISION ZERO, HALF
166: PARAMETER ( ZERO = 0.0D0, HALF = 0.5 )
1.1 bertrand 167: * ..
168: * .. Local Scalars ..
169: INTEGER J4, J4P2
1.10 bertrand 170: DOUBLE PRECISION D, EMIN, TEMP, DTHRESH
1.1 bertrand 171: * ..
172: * .. Intrinsic Functions ..
173: INTRINSIC MIN
174: * ..
175: * .. Executable Statements ..
176: *
177: IF( ( N0-I0-1 ).LE.0 )
178: $ RETURN
179: *
1.10 bertrand 180: DTHRESH = EPS*(SIGMA+TAU)
181: IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
182: IF( TAU.NE.ZERO ) THEN
1.1 bertrand 183: J4 = 4*I0 + PP - 3
184: EMIN = Z( J4+4 )
185: D = Z( J4 ) - TAU
186: DMIN = D
187: DMIN1 = -Z( J4 )
188: *
189: IF( IEEE ) THEN
190: *
191: * Code for IEEE arithmetic.
192: *
193: IF( PP.EQ.0 ) THEN
194: DO 10 J4 = 4*I0, 4*( N0-3 ), 4
195: Z( J4-2 ) = D + Z( J4-1 )
196: TEMP = Z( J4+1 ) / Z( J4-2 )
197: D = D*TEMP - TAU
198: DMIN = MIN( DMIN, D )
199: Z( J4 ) = Z( J4-1 )*TEMP
200: EMIN = MIN( Z( J4 ), EMIN )
201: 10 CONTINUE
202: ELSE
203: DO 20 J4 = 4*I0, 4*( N0-3 ), 4
204: Z( J4-3 ) = D + Z( J4 )
205: TEMP = Z( J4+2 ) / Z( J4-3 )
206: D = D*TEMP - TAU
207: DMIN = MIN( DMIN, D )
208: Z( J4-1 ) = Z( J4 )*TEMP
209: EMIN = MIN( Z( J4-1 ), EMIN )
210: 20 CONTINUE
211: END IF
212: *
213: * Unroll last two steps.
214: *
215: DNM2 = D
216: DMIN2 = DMIN
217: J4 = 4*( N0-2 ) - PP
218: J4P2 = J4 + 2*PP - 1
219: Z( J4-2 ) = DNM2 + Z( J4P2 )
220: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
221: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
222: DMIN = MIN( DMIN, DNM1 )
223: *
224: DMIN1 = DMIN
225: J4 = J4 + 4
226: J4P2 = J4 + 2*PP - 1
227: Z( J4-2 ) = DNM1 + Z( J4P2 )
228: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
229: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
230: DMIN = MIN( DMIN, DN )
231: *
232: ELSE
233: *
234: * Code for non IEEE arithmetic.
235: *
236: IF( PP.EQ.0 ) THEN
237: DO 30 J4 = 4*I0, 4*( N0-3 ), 4
238: Z( J4-2 ) = D + Z( J4-1 )
239: IF( D.LT.ZERO ) THEN
240: RETURN
241: ELSE
242: Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
243: D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
244: END IF
245: DMIN = MIN( DMIN, D )
246: EMIN = MIN( EMIN, Z( J4 ) )
247: 30 CONTINUE
248: ELSE
249: DO 40 J4 = 4*I0, 4*( N0-3 ), 4
250: Z( J4-3 ) = D + Z( J4 )
251: IF( D.LT.ZERO ) THEN
252: RETURN
253: ELSE
254: Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
255: D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
256: END IF
257: DMIN = MIN( DMIN, D )
258: EMIN = MIN( EMIN, Z( J4-1 ) )
259: 40 CONTINUE
260: END IF
261: *
262: * Unroll last two steps.
263: *
264: DNM2 = D
265: DMIN2 = DMIN
266: J4 = 4*( N0-2 ) - PP
267: J4P2 = J4 + 2*PP - 1
268: Z( J4-2 ) = DNM2 + Z( J4P2 )
269: IF( DNM2.LT.ZERO ) THEN
270: RETURN
271: ELSE
272: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
273: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
274: END IF
275: DMIN = MIN( DMIN, DNM1 )
276: *
277: DMIN1 = DMIN
278: J4 = J4 + 4
279: J4P2 = J4 + 2*PP - 1
280: Z( J4-2 ) = DNM1 + Z( J4P2 )
281: IF( DNM1.LT.ZERO ) THEN
282: RETURN
283: ELSE
284: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
285: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
286: END IF
287: DMIN = MIN( DMIN, DN )
288: *
289: END IF
1.10 bertrand 290: ELSE
291: * This is the version that sets d's to zero if they are small enough
292: J4 = 4*I0 + PP - 3
293: EMIN = Z( J4+4 )
294: D = Z( J4 ) - TAU
295: DMIN = D
296: DMIN1 = -Z( J4 )
297: IF( IEEE ) THEN
298: *
299: * Code for IEEE arithmetic.
300: *
301: IF( PP.EQ.0 ) THEN
302: DO 50 J4 = 4*I0, 4*( N0-3 ), 4
303: Z( J4-2 ) = D + Z( J4-1 )
304: TEMP = Z( J4+1 ) / Z( J4-2 )
305: D = D*TEMP - TAU
306: IF( D.LT.DTHRESH ) D = ZERO
307: DMIN = MIN( DMIN, D )
308: Z( J4 ) = Z( J4-1 )*TEMP
309: EMIN = MIN( Z( J4 ), EMIN )
310: 50 CONTINUE
311: ELSE
312: DO 60 J4 = 4*I0, 4*( N0-3 ), 4
313: Z( J4-3 ) = D + Z( J4 )
314: TEMP = Z( J4+2 ) / Z( J4-3 )
315: D = D*TEMP - TAU
316: IF( D.LT.DTHRESH ) D = ZERO
317: DMIN = MIN( DMIN, D )
318: Z( J4-1 ) = Z( J4 )*TEMP
319: EMIN = MIN( Z( J4-1 ), EMIN )
320: 60 CONTINUE
321: END IF
322: *
323: * Unroll last two steps.
324: *
325: DNM2 = D
326: DMIN2 = DMIN
327: J4 = 4*( N0-2 ) - PP
328: J4P2 = J4 + 2*PP - 1
329: Z( J4-2 ) = DNM2 + Z( J4P2 )
330: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
331: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
332: DMIN = MIN( DMIN, DNM1 )
333: *
334: DMIN1 = DMIN
335: J4 = J4 + 4
336: J4P2 = J4 + 2*PP - 1
337: Z( J4-2 ) = DNM1 + Z( J4P2 )
338: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
339: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
340: DMIN = MIN( DMIN, DN )
341: *
342: ELSE
343: *
344: * Code for non IEEE arithmetic.
345: *
346: IF( PP.EQ.0 ) THEN
347: DO 70 J4 = 4*I0, 4*( N0-3 ), 4
348: Z( J4-2 ) = D + Z( J4-1 )
349: IF( D.LT.ZERO ) THEN
350: RETURN
351: ELSE
352: Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
353: D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
354: END IF
355: IF( D.LT.DTHRESH) D = ZERO
356: DMIN = MIN( DMIN, D )
357: EMIN = MIN( EMIN, Z( J4 ) )
358: 70 CONTINUE
359: ELSE
360: DO 80 J4 = 4*I0, 4*( N0-3 ), 4
361: Z( J4-3 ) = D + Z( J4 )
362: IF( D.LT.ZERO ) THEN
363: RETURN
364: ELSE
365: Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
366: D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
367: END IF
368: IF( D.LT.DTHRESH) D = ZERO
369: DMIN = MIN( DMIN, D )
370: EMIN = MIN( EMIN, Z( J4-1 ) )
371: 80 CONTINUE
372: END IF
373: *
374: * Unroll last two steps.
375: *
376: DNM2 = D
377: DMIN2 = DMIN
378: J4 = 4*( N0-2 ) - PP
379: J4P2 = J4 + 2*PP - 1
380: Z( J4-2 ) = DNM2 + Z( J4P2 )
381: IF( DNM2.LT.ZERO ) THEN
382: RETURN
383: ELSE
384: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
385: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
386: END IF
387: DMIN = MIN( DMIN, DNM1 )
388: *
389: DMIN1 = DMIN
390: J4 = J4 + 4
391: J4P2 = J4 + 2*PP - 1
392: Z( J4-2 ) = DNM1 + Z( J4P2 )
393: IF( DNM1.LT.ZERO ) THEN
394: RETURN
395: ELSE
396: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
397: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
398: END IF
399: DMIN = MIN( DMIN, DN )
400: *
401: END IF
402: END IF
403: *
1.1 bertrand 404: Z( J4+2 ) = DN
405: Z( 4*N0-PP ) = EMIN
406: RETURN
407: *
408: * End of DLASQ5
409: *
410: END