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: *> \date December 2016
115: *
116: *> \ingroup auxOTHERcomputational
117: *
118: * =====================================================================
119: SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
120: $ DNM1, DNM2 )
121: *
122: * -- LAPACK computational routine (version 3.7.0) --
123: * -- LAPACK is a software package provided by Univ. of Tennessee, --
124: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
125: * December 2016
126: *
127: * .. Scalar Arguments ..
128: INTEGER I0, N0, PP
129: DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
130: * ..
131: * .. Array Arguments ..
132: DOUBLE PRECISION Z( * )
133: * ..
134: *
135: * =====================================================================
136: *
137: * .. Parameter ..
138: DOUBLE PRECISION ZERO
139: PARAMETER ( ZERO = 0.0D0 )
140: * ..
141: * .. Local Scalars ..
142: INTEGER J4, J4P2
143: DOUBLE PRECISION D, EMIN, SAFMIN, TEMP
144: * ..
145: * .. External Function ..
146: DOUBLE PRECISION DLAMCH
147: EXTERNAL DLAMCH
148: * ..
149: * .. Intrinsic Functions ..
150: INTRINSIC MIN
151: * ..
152: * .. Executable Statements ..
153: *
154: IF( ( N0-I0-1 ).LE.0 )
155: $ RETURN
156: *
157: SAFMIN = DLAMCH( 'Safe minimum' )
158: J4 = 4*I0 + PP - 3
159: EMIN = Z( J4+4 )
160: D = Z( J4 )
161: DMIN = D
162: *
163: IF( PP.EQ.0 ) THEN
164: DO 10 J4 = 4*I0, 4*( N0-3 ), 4
165: Z( J4-2 ) = D + Z( J4-1 )
166: IF( Z( J4-2 ).EQ.ZERO ) THEN
167: Z( J4 ) = ZERO
168: D = Z( J4+1 )
169: DMIN = D
170: EMIN = ZERO
171: ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
172: $ SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
173: TEMP = Z( J4+1 ) / Z( J4-2 )
174: Z( J4 ) = Z( J4-1 )*TEMP
175: D = D*TEMP
176: ELSE
177: Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
178: D = Z( J4+1 )*( D / Z( J4-2 ) )
179: END IF
180: DMIN = MIN( DMIN, D )
181: EMIN = MIN( EMIN, Z( J4 ) )
182: 10 CONTINUE
183: ELSE
184: DO 20 J4 = 4*I0, 4*( N0-3 ), 4
185: Z( J4-3 ) = D + Z( J4 )
186: IF( Z( J4-3 ).EQ.ZERO ) THEN
187: Z( J4-1 ) = ZERO
188: D = Z( J4+2 )
189: DMIN = D
190: EMIN = ZERO
191: ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
192: $ SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
193: TEMP = Z( J4+2 ) / Z( J4-3 )
194: Z( J4-1 ) = Z( J4 )*TEMP
195: D = D*TEMP
196: ELSE
197: Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
198: D = Z( J4+2 )*( D / Z( J4-3 ) )
199: END IF
200: DMIN = MIN( DMIN, D )
201: EMIN = MIN( EMIN, Z( J4-1 ) )
202: 20 CONTINUE
203: END IF
204: *
205: * Unroll last two steps.
206: *
207: DNM2 = D
208: DMIN2 = DMIN
209: J4 = 4*( N0-2 ) - PP
210: J4P2 = J4 + 2*PP - 1
211: Z( J4-2 ) = DNM2 + Z( J4P2 )
212: IF( Z( J4-2 ).EQ.ZERO ) THEN
213: Z( J4 ) = ZERO
214: DNM1 = Z( J4P2+2 )
215: DMIN = DNM1
216: EMIN = ZERO
217: ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
218: $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
219: TEMP = Z( J4P2+2 ) / Z( J4-2 )
220: Z( J4 ) = Z( J4P2 )*TEMP
221: DNM1 = DNM2*TEMP
222: ELSE
223: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
224: DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
225: END IF
226: DMIN = MIN( DMIN, DNM1 )
227: *
228: DMIN1 = DMIN
229: J4 = J4 + 4
230: J4P2 = J4 + 2*PP - 1
231: Z( J4-2 ) = DNM1 + Z( J4P2 )
232: IF( Z( J4-2 ).EQ.ZERO ) THEN
233: Z( J4 ) = ZERO
234: DN = Z( J4P2+2 )
235: DMIN = DN
236: EMIN = ZERO
237: ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
238: $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
239: TEMP = Z( J4P2+2 ) / Z( J4-2 )
240: Z( J4 ) = Z( J4P2 )*TEMP
241: DN = DNM1*TEMP
242: ELSE
243: Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
244: DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
245: END IF
246: DMIN = MIN( DMIN, DN )
247: *
248: Z( J4+2 ) = DN
249: Z( 4*N0-PP ) = EMIN
250: RETURN
251: *
252: * End of DLASQ6
253: *
254: END
CVSweb interface <joel.bertrand@systella.fr>