1: *> \brief \b DLASQ6 computes one dqd 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 DLASQ6 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq6.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq6.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq6.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
22: * DNM1, DNM2 )
23: *
24: * .. Scalar Arguments ..
25: * INTEGER I0, N0, PP
26: * DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION Z( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DLASQ6 computes one dqd (shift equal to zero) transform in
39: *> ping-pong form, with protection against underflow and overflow.
40: *> \endverbatim
41: *
42: * Arguments:
43: * ==========
44: *
45: *> \param[in] I0
46: *> \verbatim
47: *> I0 is INTEGER
48: *> First index.
49: *> \endverbatim
50: *>
51: *> \param[in] N0
52: *> \verbatim
53: *> N0 is INTEGER
54: *> Last index.
55: *> \endverbatim
56: *>
57: *> \param[in] Z
58: *> \verbatim
59: *> Z is DOUBLE PRECISION array, dimension ( 4*N )
60: *> Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
61: *> an extra argument.
62: *> \endverbatim
63: *>
64: *> \param[in] PP
65: *> \verbatim
66: *> PP is INTEGER
67: *> PP=0 for ping, PP=1 for pong.
68: *> \endverbatim
69: *>
70: *> \param[out] DMIN
71: *> \verbatim
72: *> DMIN is DOUBLE PRECISION
73: *> Minimum value of d.
74: *> \endverbatim
75: *>
76: *> \param[out] DMIN1
77: *> \verbatim
78: *> DMIN1 is DOUBLE PRECISION
79: *> Minimum value of d, excluding D( N0 ).
80: *> \endverbatim
81: *>
82: *> \param[out] DMIN2
83: *> \verbatim
84: *> DMIN2 is DOUBLE PRECISION
85: *> Minimum value of d, excluding D( N0 ) and D( N0-1 ).
86: *> \endverbatim
87: *>
88: *> \param[out] DN
89: *> \verbatim
90: *> DN is DOUBLE PRECISION
91: *> d(N0), the last value of d.
92: *> \endverbatim
93: *>
94: *> \param[out] DNM1
95: *> \verbatim
96: *> DNM1 is DOUBLE PRECISION
97: *> d(N0-1).
98: *> \endverbatim
99: *>
100: *> \param[out] DNM2
101: *> \verbatim
102: *> DNM2 is DOUBLE PRECISION
103: *> d(N0-2).
104: *> \endverbatim
105: *
106: * Authors:
107: * ========
108: *
109: *> \author Univ. of Tennessee
110: *> \author Univ. of California Berkeley
111: *> \author Univ. of Colorado Denver
112: *> \author NAG Ltd.
113: *
114: *> \ingroup auxOTHERcomputational
115: *
116: * =====================================================================
117: SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
118: $ DNM1, DNM2 )
119: *
120: * -- LAPACK computational routine --
121: * -- LAPACK is a software package provided by Univ. of Tennessee, --
122: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
123: *
124: * .. Scalar Arguments ..
125: INTEGER I0, N0, PP
126: DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
127: * ..
128: * .. Array Arguments ..
129: DOUBLE PRECISION Z( * )
130: * ..
131: *
132: * =====================================================================
133: *
134: * .. Parameter ..
135: DOUBLE PRECISION ZERO
136: PARAMETER ( ZERO = 0.0D0 )
137: * ..
138: * .. Local Scalars ..
139: INTEGER J4, J4P2
140: DOUBLE PRECISION D, EMIN, SAFMIN, TEMP
141: * ..
142: * .. External Function ..
143: DOUBLE PRECISION DLAMCH
144: EXTERNAL DLAMCH
145: * ..
146: * .. Intrinsic Functions ..
147: INTRINSIC MIN
148: * ..
149: * .. Executable Statements ..
150: *
151: IF( ( N0-I0-1 ).LE.0 )
152: $ RETURN
153: *
154: SAFMIN = DLAMCH( 'Safe minimum' )
155: J4 = 4*I0 + PP - 3
156: EMIN = Z( J4+4 )
157: D = Z( J4 )
158: DMIN = D
159: *
160: IF( PP.EQ.0 ) THEN
161: DO 10 J4 = 4*I0, 4*( N0-3 ), 4
162: Z( J4-2 ) = D + Z( J4-1 )
163: IF( Z( J4-2 ).EQ.ZERO ) THEN
164: Z( J4 ) = ZERO
165: D = Z( J4+1 )
166: DMIN = D
167: EMIN = ZERO
168: ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
169: $ SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
170: TEMP = Z( J4+1 ) / Z( J4-2 )
171: Z( J4 ) = Z( J4-1 )*TEMP
172: D = D*TEMP
173: ELSE
174: Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
175: D = Z( J4+1 )*( D / Z( J4-2 ) )
176: END IF
177: DMIN = MIN( DMIN, D )
178: EMIN = MIN( EMIN, Z( J4 ) )
179: 10 CONTINUE
180: ELSE
181: DO 20 J4 = 4*I0, 4*( N0-3 ), 4
182: Z( J4-3 ) = D + Z( J4 )
183: IF( Z( J4-3 ).EQ.ZERO ) THEN
184: Z( J4-1 ) = ZERO
185: D = Z( J4+2 )
186: DMIN = D
187: EMIN = ZERO
188: ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
189: $ SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
190: TEMP = Z( J4+2 ) / Z( J4-3 )
191: Z( J4-1 ) = Z( J4 )*TEMP
192: D = D*TEMP
193: ELSE
194: Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
195: D = Z( J4+2 )*( D / Z( J4-3 ) )
196: END IF
197: DMIN = MIN( DMIN, D )
198: EMIN = MIN( EMIN, Z( J4-1 ) )
199: 20 CONTINUE
200: END IF
201: *
202: * Unroll last two steps.
203: *
204: DNM2 = D
205: DMIN2 = DMIN
206: J4 = 4*( N0-2 ) - PP
207: J4P2 = J4 + 2*PP - 1
208: Z( J4-2 ) = DNM2 + Z( J4P2 )
209: IF( Z( J4-2 ).EQ.ZERO ) THEN
210: Z( J4 ) = ZERO
211: DNM1 = Z( J4P2+2 )
212: DMIN = DNM1
213: EMIN = ZERO
214: ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
215: $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
216: TEMP = Z( J4P2+2 ) / Z( J4-2 )
217: Z( J4 ) = Z( J4P2 )*TEMP
218: DNM1 = DNM2*TEMP
219: ELSE
220: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
221: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
222: END IF
223: DMIN = MIN( DMIN, DNM1 )
224: *
225: DMIN1 = DMIN
226: J4 = J4 + 4
227: J4P2 = J4 + 2*PP - 1
228: Z( J4-2 ) = DNM1 + Z( J4P2 )
229: IF( Z( J4-2 ).EQ.ZERO ) THEN
230: Z( J4 ) = ZERO
231: DN = Z( J4P2+2 )
232: DMIN = DN
233: EMIN = ZERO
234: ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
235: $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
236: TEMP = Z( J4P2+2 ) / Z( J4-2 )
237: Z( J4 ) = Z( J4P2 )*TEMP
238: DN = DNM1*TEMP
239: ELSE
240: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
241: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
242: END IF
243: DMIN = MIN( DMIN, DN )
244: *
245: Z( J4+2 ) = DN
246: Z( 4*N0-PP ) = EMIN
247: RETURN
248: *
249: * End of DLASQ6
250: *
251: END
CVSweb interface <joel.bertrand@systella.fr>