version 1.7, 2010/12/21 13:53:33
|
version 1.16, 2017/06/17 11:06:27
|
Line 1
|
Line 1
|
SUBROUTINE DLASQ2( N, Z, INFO ) |
*> \brief \b DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr. |
* |
* |
* -- LAPACK routine (version 3.2) -- |
* =========== DOCUMENTATION =========== |
* |
* |
* -- Contributed by Osni Marques of the Lawrence Berkeley National -- |
* Online html documentation available at |
* -- Laboratory and Beresford Parlett of the Univ. of California at -- |
* http://www.netlib.org/lapack/explore-html/ |
* -- Berkeley -- |
|
* -- November 2008 -- |
|
* |
* |
|
*> \htmlonly |
|
*> Download DLASQ2 + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DLASQ2( N, Z, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* INTEGER INFO, N |
|
* .. |
|
* .. Array Arguments .. |
|
* DOUBLE PRECISION Z( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DLASQ2 computes all the eigenvalues of the symmetric positive |
|
*> definite tridiagonal matrix associated with the qd array Z to high |
|
*> relative accuracy are computed to high relative accuracy, in the |
|
*> absence of denormalization, underflow and overflow. |
|
*> |
|
*> To see the relation of Z to the tridiagonal matrix, let L be a |
|
*> unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and |
|
*> let U be an upper bidiagonal matrix with 1's above and diagonal |
|
*> Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the |
|
*> symmetric tridiagonal to which it is similar. |
|
*> |
|
*> Note : DLASQ2 defines a logical variable, IEEE, which is true |
|
*> on machines which follow ieee-754 floating-point standard in their |
|
*> handling of infinities and NaNs, and false otherwise. This variable |
|
*> is passed to DLASQ3. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] N |
|
*> \verbatim |
|
*> N is INTEGER |
|
*> The number of rows and columns in the matrix. N >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] Z |
|
*> \verbatim |
|
*> Z is DOUBLE PRECISION array, dimension ( 4*N ) |
|
*> On entry Z holds the qd array. On exit, entries 1 to N hold |
|
*> the eigenvalues in decreasing order, Z( 2*N+1 ) holds the |
|
*> trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If |
|
*> N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) |
|
*> holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of |
|
*> shifts that failed. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit |
|
*> < 0: if the i-th argument is a scalar and had an illegal |
|
*> value, then INFO = -i, if the i-th argument is an |
|
*> array and the j-entry had an illegal value, then |
|
*> INFO = -(i*100+j) |
|
*> > 0: the algorithm failed |
|
*> = 1, a split was marked by a positive value in E |
|
*> = 2, current block of Z not diagonalized after 100*N |
|
*> iterations (in inner while loop). On exit Z holds |
|
*> a qd array with the same eigenvalues as the given Z. |
|
*> = 3, termination criterion of outer while loop not met |
|
*> (program created more than N unreduced blocks) |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date December 2016 |
|
* |
|
*> \ingroup auxOTHERcomputational |
|
* |
|
*> \par Further Details: |
|
* ===================== |
|
*> |
|
*> \verbatim |
|
*> |
|
*> Local Variables: I0:N0 defines a current unreduced segment of Z. |
|
*> The shifts are accumulated in SIGMA. Iteration count is in ITER. |
|
*> Ping-pong is controlled by PP (alternates between 0 and 1). |
|
*> \endverbatim |
|
*> |
|
* ===================================================================== |
|
SUBROUTINE DLASQ2( N, Z, INFO ) |
|
* |
|
* -- LAPACK computational routine (version 3.7.0) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
|
* December 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER INFO, N |
INTEGER INFO, N |
Line 17
|
Line 124
|
DOUBLE PRECISION Z( * ) |
DOUBLE PRECISION Z( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* DLASQ2 computes all the eigenvalues of the symmetric positive |
|
* definite tridiagonal matrix associated with the qd array Z to high |
|
* relative accuracy are computed to high relative accuracy, in the |
|
* absence of denormalization, underflow and overflow. |
|
* |
|
* To see the relation of Z to the tridiagonal matrix, let L be a |
|
* unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and |
|
* let U be an upper bidiagonal matrix with 1's above and diagonal |
|
* Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the |
|
* symmetric tridiagonal to which it is similar. |
|
* |
|
* Note : DLASQ2 defines a logical variable, IEEE, which is true |
|
* on machines which follow ieee-754 floating-point standard in their |
|
* handling of infinities and NaNs, and false otherwise. This variable |
|
* is passed to DLASQ3. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* N (input) INTEGER |
|
* The number of rows and columns in the matrix. N >= 0. |
|
* |
|
* Z (input/output) DOUBLE PRECISION array, dimension ( 4*N ) |
|
* On entry Z holds the qd array. On exit, entries 1 to N hold |
|
* the eigenvalues in decreasing order, Z( 2*N+1 ) holds the |
|
* trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If |
|
* N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) |
|
* holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of |
|
* shifts that failed. |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit |
|
* < 0: if the i-th argument is a scalar and had an illegal |
|
* value, then INFO = -i, if the i-th argument is an |
|
* array and the j-entry had an illegal value, then |
|
* INFO = -(i*100+j) |
|
* > 0: the algorithm failed |
|
* = 1, a split was marked by a positive value in E |
|
* = 2, current block of Z not diagonalized after 30*N |
|
* iterations (in inner while loop) |
|
* = 3, termination criterion of outer while loop not met |
|
* (program created more than N unreduced blocks) |
|
* |
|
* Further Details |
|
* =============== |
|
* Local Variables: I0:N0 defines a current unreduced segment of Z. |
|
* The shifts are accumulated in SIGMA. Iteration count is in ITER. |
|
* Ping-pong is controlled by PP (alternates between 0 and 1). |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 80
|
Line 135
|
* .. |
* .. |
* .. Local Scalars .. |
* .. Local Scalars .. |
LOGICAL IEEE |
LOGICAL IEEE |
INTEGER I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K, |
INTEGER I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, |
$ KMIN, N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE |
$ K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT, |
|
$ TTYPE |
DOUBLE PRECISION D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN, |
DOUBLE PRECISION D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN, |
$ DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX, |
$ DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX, |
$ QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL, |
$ QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL, |
$ TOL2, TRACE, ZMAX |
$ TOL2, TRACE, ZMAX, TEMPE, TEMPQ |
* .. |
* .. |
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL DLASQ3, DLASRT, XERBLA |
EXTERNAL DLASQ3, DLASRT, XERBLA |
Line 99
|
Line 155
|
INTRINSIC ABS, DBLE, MAX, MIN, SQRT |
INTRINSIC ABS, DBLE, MAX, MIN, SQRT |
* .. |
* .. |
* .. Executable Statements .. |
* .. Executable Statements .. |
* |
* |
* Test the input arguments. |
* Test the input arguments. |
* (in case DLASQ2 is not called by DLASQ1) |
* (in case DLASQ2 is not called by DLASQ1) |
* |
* |
Line 139
|
Line 195
|
END IF |
END IF |
Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 ) |
Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 ) |
IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN |
IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN |
T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) |
T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) |
S = Z( 3 )*( Z( 2 ) / T ) |
S = Z( 3 )*( Z( 2 ) / T ) |
IF( S.LE.T ) THEN |
IF( S.LE.T ) THEN |
S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) ) |
S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) ) |
Line 208
|
Line 264
|
Z( 2*N-1 ) = ZERO |
Z( 2*N-1 ) = ZERO |
RETURN |
RETURN |
END IF |
END IF |
* |
* |
* Check whether the machine is IEEE conformable. |
* Check whether the machine is IEEE conformable. |
* |
* |
IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND. |
IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND. |
$ ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 |
$ ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 |
* |
* |
* Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). |
* Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). |
* |
* |
DO 30 K = 2*N, 2, -2 |
DO 30 K = 2*N, 2, -2 |
Z( 2*K ) = ZERO |
Z( 2*K ) = ZERO |
Z( 2*K-1 ) = Z( K ) |
Z( 2*K-1 ) = Z( K ) |
Z( 2*K-2 ) = ZERO |
Z( 2*K-2 ) = ZERO |
Z( 2*K-3 ) = Z( K-1 ) |
Z( 2*K-3 ) = Z( K-1 ) |
30 CONTINUE |
30 CONTINUE |
* |
* |
I0 = 1 |
I0 = 1 |
Line 277
|
Line 333
|
D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) ) |
D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) ) |
END IF |
END IF |
EMIN = MIN( EMIN, Z( I4-2*PP ) ) |
EMIN = MIN( EMIN, Z( I4-2*PP ) ) |
60 CONTINUE |
60 CONTINUE |
Z( 4*N0-PP-2 ) = D |
Z( 4*N0-PP-2 ) = D |
* |
* |
* Now find qmax. |
* Now find qmax. |
Line 308
|
Line 364
|
NDIV = 2*( N0-I0 ) |
NDIV = 2*( N0-I0 ) |
* |
* |
DO 160 IWHILA = 1, N + 1 |
DO 160 IWHILA = 1, N + 1 |
IF( N0.LT.1 ) |
IF( N0.LT.1 ) |
$ GO TO 170 |
$ GO TO 170 |
* |
* |
* While array unfinished do |
* While array unfinished do |
* |
* |
* E(N0) holds the value of SIGMA when submatrix in I0:N0 |
* E(N0) holds the value of SIGMA when submatrix in I0:N0 |
* splits from the rest of the array, but is negated. |
* splits from the rest of the array, but is negated. |
* |
* |
DESIG = ZERO |
DESIG = ZERO |
IF( N0.EQ.N ) THEN |
IF( N0.EQ.N ) THEN |
SIGMA = ZERO |
SIGMA = ZERO |
Line 330
|
Line 386
|
* Find last unreduced submatrix's top index I0, find QMAX and |
* Find last unreduced submatrix's top index I0, find QMAX and |
* EMIN. Find Gershgorin-type bound if Q's much greater than E's. |
* EMIN. Find Gershgorin-type bound if Q's much greater than E's. |
* |
* |
EMAX = ZERO |
EMAX = ZERO |
IF( N0.GT.I0 ) THEN |
IF( N0.GT.I0 ) THEN |
EMIN = ABS( Z( 4*N0-5 ) ) |
EMIN = ABS( Z( 4*N0-5 ) ) |
ELSE |
ELSE |
Line 348
|
Line 404
|
QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) ) |
QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) ) |
EMIN = MIN( EMIN, Z( I4-5 ) ) |
EMIN = MIN( EMIN, Z( I4-5 ) ) |
90 CONTINUE |
90 CONTINUE |
I4 = 4 |
I4 = 4 |
* |
* |
100 CONTINUE |
100 CONTINUE |
I0 = I4 / 4 |
I0 = I4 / 4 |
Line 365
|
Line 421
|
KMIN = ( I4+3 )/4 |
KMIN = ( I4+3 )/4 |
END IF |
END IF |
110 CONTINUE |
110 CONTINUE |
IF( (KMIN-I0)*2.LT.N0-KMIN .AND. |
IF( (KMIN-I0)*2.LT.N0-KMIN .AND. |
$ DEEMIN.LE.HALF*Z(4*N0-3) ) THEN |
$ DEEMIN.LE.HALF*Z(4*N0-3) ) THEN |
IPN4 = 4*( I0+N0 ) |
IPN4 = 4*( I0+N0 ) |
PP = 2 |
PP = 2 |
Line 390
|
Line 446
|
* |
* |
DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) ) |
DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) ) |
* |
* |
* Now I0:N0 is unreduced. |
* Now I0:N0 is unreduced. |
* PP = 0 for ping, PP = 1 for pong. |
* PP = 0 for ping, PP = 1 for pong. |
* PP = 2 indicates that flipping was applied to the Z array and |
* PP = 2 indicates that flipping was applied to the Z array and |
* and that the tests for deflation upon entry in DLASQ3 |
* and that the tests for deflation upon entry in DLASQ3 |
* should not be performed. |
* should not be performed. |
* |
* |
NBIG = 30*( N0-I0+1 ) |
NBIG = 100*( N0-I0+1 ) |
DO 140 IWHILB = 1, NBIG |
DO 140 IWHILB = 1, NBIG |
IF( I0.GT.N0 ) |
IF( I0.GT.N0 ) |
$ GO TO 150 |
$ GO TO 150 |
* |
* |
* While submatrix unfinished take a good dqds step. |
* While submatrix unfinished take a good dqds step. |
Line 441
|
Line 497
|
140 CONTINUE |
140 CONTINUE |
* |
* |
INFO = 2 |
INFO = 2 |
|
* |
|
* Maximum number of iterations exceeded, restore the shift |
|
* SIGMA and place the new d's and e's in a qd array. |
|
* This might need to be done for several blocks |
|
* |
|
I1 = I0 |
|
N1 = N0 |
|
145 CONTINUE |
|
TEMPQ = Z( 4*I0-3 ) |
|
Z( 4*I0-3 ) = Z( 4*I0-3 ) + SIGMA |
|
DO K = I0+1, N0 |
|
TEMPE = Z( 4*K-5 ) |
|
Z( 4*K-5 ) = Z( 4*K-5 ) * (TEMPQ / Z( 4*K-7 )) |
|
TEMPQ = Z( 4*K-3 ) |
|
Z( 4*K-3 ) = Z( 4*K-3 ) + SIGMA + TEMPE - Z( 4*K-5 ) |
|
END DO |
|
* |
|
* Prepare to do this on the previous block if there is one |
|
* |
|
IF( I1.GT.1 ) THEN |
|
N1 = I1-1 |
|
DO WHILE( ( I1.GE.2 ) .AND. ( Z(4*I1-5).GE.ZERO ) ) |
|
I1 = I1 - 1 |
|
END DO |
|
SIGMA = -Z(4*N1-1) |
|
GO TO 145 |
|
END IF |
|
|
|
DO K = 1, N |
|
Z( 2*K-1 ) = Z( 4*K-3 ) |
|
* |
|
* Only the block 1..N0 is unfinished. The rest of the e's |
|
* must be essentially zero, although sometimes other data |
|
* has been stored in them. |
|
* |
|
IF( K.LT.N0 ) THEN |
|
Z( 2*K ) = Z( 4*K-1 ) |
|
ELSE |
|
Z( 2*K ) = 0 |
|
END IF |
|
END DO |
RETURN |
RETURN |
* |
* |
* end IWHILB |
* end IWHILB |
Line 452
|
Line 549
|
INFO = 3 |
INFO = 3 |
RETURN |
RETURN |
* |
* |
* end IWHILA |
* end IWHILA |
* |
* |
170 CONTINUE |
170 CONTINUE |
* |
* |
* Move q's to the front. |
* Move q's to the front. |
* |
* |
DO 180 K = 2, N |
DO 180 K = 2, N |
Z( K ) = Z( 4*K-3 ) |
Z( K ) = Z( 4*K-3 ) |
180 CONTINUE |
180 CONTINUE |
* |
* |
* Sort and compute sum of eigenvalues. |
* Sort and compute sum of eigenvalues. |
* |
* |
CALL DLASRT( 'D', N, Z, IINFO ) |
CALL DLASRT( 'D', N, Z, IINFO ) |
Line 473
|
Line 570
|
* |
* |
* Store trace, sum(eigenvalues) and information on performance. |
* Store trace, sum(eigenvalues) and information on performance. |
* |
* |
Z( 2*N+1 ) = TRACE |
Z( 2*N+1 ) = TRACE |
Z( 2*N+2 ) = E |
Z( 2*N+2 ) = E |
Z( 2*N+3 ) = DBLE( ITER ) |
Z( 2*N+3 ) = DBLE( ITER ) |
Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 ) |
Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 ) |