Annotation of rpl/lapack/lapack/dlasq3.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLASQ3
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLASQ3 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq3.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq3.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq3.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
! 22: * ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
! 23: * DN2, G, TAU )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * LOGICAL IEEE
! 27: * INTEGER I0, ITER, N0, NDIV, NFAIL, PP
! 28: * DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
! 29: * $ QMAX, SIGMA, TAU
! 30: * ..
! 31: * .. Array Arguments ..
! 32: * DOUBLE PRECISION Z( * )
! 33: * ..
! 34: *
! 35: *
! 36: *> \par Purpose:
! 37: * =============
! 38: *>
! 39: *> \verbatim
! 40: *>
! 41: *> DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
! 42: *> In case of failure it changes shifts, and tries again until output
! 43: *> is positive.
! 44: *> \endverbatim
! 45: *
! 46: * Arguments:
! 47: * ==========
! 48: *
! 49: *> \param[in] I0
! 50: *> \verbatim
! 51: *> I0 is INTEGER
! 52: *> First index.
! 53: *> \endverbatim
! 54: *>
! 55: *> \param[in,out] N0
! 56: *> \verbatim
! 57: *> N0 is INTEGER
! 58: *> Last index.
! 59: *> \endverbatim
! 60: *>
! 61: *> \param[in] Z
! 62: *> \verbatim
! 63: *> Z is DOUBLE PRECISION array, dimension ( 4*N )
! 64: *> Z holds the qd array.
! 65: *> \endverbatim
! 66: *>
! 67: *> \param[in,out] PP
! 68: *> \verbatim
! 69: *> PP is INTEGER
! 70: *> PP=0 for ping, PP=1 for pong.
! 71: *> PP=2 indicates that flipping was applied to the Z array
! 72: *> and that the initial tests for deflation should not be
! 73: *> performed.
! 74: *> \endverbatim
! 75: *>
! 76: *> \param[out] DMIN
! 77: *> \verbatim
! 78: *> DMIN is DOUBLE PRECISION
! 79: *> Minimum value of d.
! 80: *> \endverbatim
! 81: *>
! 82: *> \param[out] SIGMA
! 83: *> \verbatim
! 84: *> SIGMA is DOUBLE PRECISION
! 85: *> Sum of shifts used in current segment.
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[in,out] DESIG
! 89: *> \verbatim
! 90: *> DESIG is DOUBLE PRECISION
! 91: *> Lower order part of SIGMA
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[in] QMAX
! 95: *> \verbatim
! 96: *> QMAX is DOUBLE PRECISION
! 97: *> Maximum value of q.
! 98: *> \endverbatim
! 99: *>
! 100: *> \param[out] NFAIL
! 101: *> \verbatim
! 102: *> NFAIL is INTEGER
! 103: *> Number of times shift was too big.
! 104: *> \endverbatim
! 105: *>
! 106: *> \param[out] ITER
! 107: *> \verbatim
! 108: *> ITER is INTEGER
! 109: *> Number of iterations.
! 110: *> \endverbatim
! 111: *>
! 112: *> \param[out] NDIV
! 113: *> \verbatim
! 114: *> NDIV is INTEGER
! 115: *> Number of divisions.
! 116: *> \endverbatim
! 117: *>
! 118: *> \param[in] IEEE
! 119: *> \verbatim
! 120: *> IEEE is LOGICAL
! 121: *> Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
! 122: *> \endverbatim
! 123: *>
! 124: *> \param[in,out] TTYPE
! 125: *> \verbatim
! 126: *> TTYPE is INTEGER
! 127: *> Shift type.
! 128: *> \endverbatim
! 129: *>
! 130: *> \param[in,out] DMIN1
! 131: *> \verbatim
! 132: *> DMIN1 is DOUBLE PRECISION
! 133: *> \endverbatim
! 134: *>
! 135: *> \param[in,out] DMIN2
! 136: *> \verbatim
! 137: *> DMIN2 is DOUBLE PRECISION
! 138: *> \endverbatim
! 139: *>
! 140: *> \param[in,out] DN
! 141: *> \verbatim
! 142: *> DN is DOUBLE PRECISION
! 143: *> \endverbatim
! 144: *>
! 145: *> \param[in,out] DN1
! 146: *> \verbatim
! 147: *> DN1 is DOUBLE PRECISION
! 148: *> \endverbatim
! 149: *>
! 150: *> \param[in,out] DN2
! 151: *> \verbatim
! 152: *> DN2 is DOUBLE PRECISION
! 153: *> \endverbatim
! 154: *>
! 155: *> \param[in,out] G
! 156: *> \verbatim
! 157: *> G is DOUBLE PRECISION
! 158: *> \endverbatim
! 159: *>
! 160: *> \param[in,out] TAU
! 161: *> \verbatim
! 162: *> TAU is DOUBLE PRECISION
! 163: *>
! 164: *> These are passed as arguments in order to save their values
! 165: *> between calls to DLASQ3.
! 166: *> \endverbatim
! 167: *
! 168: * Authors:
! 169: * ========
! 170: *
! 171: *> \author Univ. of Tennessee
! 172: *> \author Univ. of California Berkeley
! 173: *> \author Univ. of Colorado Denver
! 174: *> \author NAG Ltd.
! 175: *
! 176: *> \date November 2011
! 177: *
! 178: *> \ingroup auxOTHERcomputational
! 179: *
! 180: * =====================================================================
1.1 bertrand 181: SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
182: $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
183: $ DN2, G, TAU )
184: *
1.9 ! bertrand 185: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 186: * -- LAPACK is a software package provided by Univ. of Tennessee, --
187: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 188: * November 2011
1.1 bertrand 189: *
190: * .. Scalar Arguments ..
191: LOGICAL IEEE
192: INTEGER I0, ITER, N0, NDIV, NFAIL, PP
193: DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
194: $ QMAX, SIGMA, TAU
195: * ..
196: * .. Array Arguments ..
197: DOUBLE PRECISION Z( * )
198: * ..
199: *
200: * =====================================================================
201: *
202: * .. Parameters ..
203: DOUBLE PRECISION CBIAS
204: PARAMETER ( CBIAS = 1.50D0 )
205: DOUBLE PRECISION ZERO, QURTR, HALF, ONE, TWO, HUNDRD
206: PARAMETER ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
207: $ ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
208: * ..
209: * .. Local Scalars ..
210: INTEGER IPN4, J4, N0IN, NN, TTYPE
211: DOUBLE PRECISION EPS, S, T, TEMP, TOL, TOL2
212: * ..
213: * .. External Subroutines ..
214: EXTERNAL DLASQ4, DLASQ5, DLASQ6
215: * ..
216: * .. External Function ..
217: DOUBLE PRECISION DLAMCH
218: LOGICAL DISNAN
219: EXTERNAL DISNAN, DLAMCH
220: * ..
221: * .. Intrinsic Functions ..
222: INTRINSIC ABS, MAX, MIN, SQRT
223: * ..
224: * .. Executable Statements ..
225: *
226: N0IN = N0
227: EPS = DLAMCH( 'Precision' )
228: TOL = EPS*HUNDRD
229: TOL2 = TOL**2
230: *
231: * Check for deflation.
232: *
233: 10 CONTINUE
234: *
235: IF( N0.LT.I0 )
236: $ RETURN
237: IF( N0.EQ.I0 )
238: $ GO TO 20
239: NN = 4*N0 + PP
240: IF( N0.EQ.( I0+1 ) )
241: $ GO TO 40
242: *
243: * Check whether E(N0-1) is negligible, 1 eigenvalue.
244: *
245: IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
246: $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
247: $ GO TO 30
248: *
249: 20 CONTINUE
250: *
251: Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
252: N0 = N0 - 1
253: GO TO 10
254: *
255: * Check whether E(N0-2) is negligible, 2 eigenvalues.
256: *
257: 30 CONTINUE
258: *
259: IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
260: $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
261: $ GO TO 50
262: *
263: 40 CONTINUE
264: *
265: IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
266: S = Z( NN-3 )
267: Z( NN-3 ) = Z( NN-7 )
268: Z( NN-7 ) = S
269: END IF
270: IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
271: T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
272: S = Z( NN-3 )*( Z( NN-5 ) / T )
273: IF( S.LE.T ) THEN
274: S = Z( NN-3 )*( Z( NN-5 ) /
275: $ ( T*( ONE+SQRT( ONE+S / T ) ) ) )
276: ELSE
277: S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
278: END IF
279: T = Z( NN-7 ) + ( S+Z( NN-5 ) )
280: Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
281: Z( NN-7 ) = T
282: END IF
283: Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
284: Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
285: N0 = N0 - 2
286: GO TO 10
287: *
288: 50 CONTINUE
289: IF( PP.EQ.2 )
290: $ PP = 0
291: *
292: * Reverse the qd-array, if warranted.
293: *
294: IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
295: IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
296: IPN4 = 4*( I0+N0 )
297: DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
298: TEMP = Z( J4-3 )
299: Z( J4-3 ) = Z( IPN4-J4-3 )
300: Z( IPN4-J4-3 ) = TEMP
301: TEMP = Z( J4-2 )
302: Z( J4-2 ) = Z( IPN4-J4-2 )
303: Z( IPN4-J4-2 ) = TEMP
304: TEMP = Z( J4-1 )
305: Z( J4-1 ) = Z( IPN4-J4-5 )
306: Z( IPN4-J4-5 ) = TEMP
307: TEMP = Z( J4 )
308: Z( J4 ) = Z( IPN4-J4-4 )
309: Z( IPN4-J4-4 ) = TEMP
310: 60 CONTINUE
311: IF( N0-I0.LE.4 ) THEN
312: Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
313: Z( 4*N0-PP ) = Z( 4*I0-PP )
314: END IF
315: DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
316: Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
317: $ Z( 4*I0+PP+3 ) )
318: Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
319: $ Z( 4*I0-PP+4 ) )
320: QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
321: DMIN = -ZERO
322: END IF
323: END IF
324: *
325: * Choose a shift.
326: *
327: CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
328: $ DN2, TAU, TTYPE, G )
329: *
330: * Call dqds until DMIN > 0.
331: *
332: 70 CONTINUE
333: *
334: CALL DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
335: $ DN1, DN2, IEEE )
336: *
337: NDIV = NDIV + ( N0-I0+2 )
338: ITER = ITER + 1
339: *
340: * Check status.
341: *
342: IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
343: *
344: * Success.
345: *
346: GO TO 90
347: *
348: ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
349: $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
350: $ ABS( DN ).LT.TOL*SIGMA ) THEN
351: *
352: * Convergence hidden by negative DN.
353: *
354: Z( 4*( N0-1 )-PP+2 ) = ZERO
355: DMIN = ZERO
356: GO TO 90
357: ELSE IF( DMIN.LT.ZERO ) THEN
358: *
359: * TAU too big. Select new TAU and try again.
360: *
361: NFAIL = NFAIL + 1
362: IF( TTYPE.LT.-22 ) THEN
363: *
364: * Failed twice. Play it safe.
365: *
366: TAU = ZERO
367: ELSE IF( DMIN1.GT.ZERO ) THEN
368: *
369: * Late failure. Gives excellent shift.
370: *
371: TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
372: TTYPE = TTYPE - 11
373: ELSE
374: *
375: * Early failure. Divide by 4.
376: *
377: TAU = QURTR*TAU
378: TTYPE = TTYPE - 12
379: END IF
380: GO TO 70
381: ELSE IF( DISNAN( DMIN ) ) THEN
382: *
383: * NaN.
384: *
385: IF( TAU.EQ.ZERO ) THEN
386: GO TO 80
387: ELSE
388: TAU = ZERO
389: GO TO 70
390: END IF
391: ELSE
392: *
393: * Possible underflow. Play it safe.
394: *
395: GO TO 80
396: END IF
397: *
398: * Risk of underflow.
399: *
400: 80 CONTINUE
401: CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
402: NDIV = NDIV + ( N0-I0+2 )
403: ITER = ITER + 1
404: TAU = ZERO
405: *
406: 90 CONTINUE
407: IF( TAU.LT.SIGMA ) THEN
408: DESIG = DESIG + TAU
409: T = SIGMA + DESIG
410: DESIG = DESIG - ( T-SIGMA )
411: ELSE
412: T = SIGMA + TAU
413: DESIG = SIGMA - ( T-TAU ) + DESIG
414: END IF
415: SIGMA = T
416: *
417: RETURN
418: *
419: * End of DLASQ3
420: *
421: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack