Annotation of rpl/lapack/lapack/dlasq5.f, revision 1.10
1.8 bertrand 1: *> \brief \b DLASQ5
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.10 ! bertrand 139: *> \date April 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.10 ! bertrand 147: * -- LAPACK computational routine (version 3.4.1) --
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.10 ! bertrand 150: * April 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
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