version 1.14, 2016/08/27 15:34:31
|
version 1.17, 2018/05/29 07:18:01
|
Line 2
|
Line 2
|
* |
* |
* =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION =========== |
* |
* |
* Online html documentation available at |
* Online html documentation available at |
* http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
*> \htmlonly |
*> \htmlonly |
*> Download DLASQ2 + dependencies |
*> Download DLASQ2 + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
* |
* |
* SUBROUTINE DLASQ2( N, Z, INFO ) |
* SUBROUTINE DLASQ2( N, Z, INFO ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* INTEGER INFO, N |
* INTEGER INFO, N |
* .. |
* .. |
* .. Array Arguments .. |
* .. Array Arguments .. |
* DOUBLE PRECISION Z( * ) |
* DOUBLE PRECISION Z( * ) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
*> |
*> |
*> \verbatim |
*> \verbatim |
*> |
*> |
*> DLASQ2 computes all the eigenvalues of the symmetric positive |
*> DLASQ2 computes all the eigenvalues of the symmetric positive |
*> definite tridiagonal matrix associated with the qd array Z to high |
*> definite tridiagonal matrix associated with the qd array Z to high |
*> relative accuracy are computed to high relative accuracy, in the |
*> relative accuracy are computed to high relative accuracy, in the |
*> absence of denormalization, underflow and overflow. |
*> absence of denormalization, underflow and overflow. |
Line 83
|
Line 83
|
*> = 2, current block of Z not diagonalized after 100*N |
*> = 2, current block of Z not diagonalized after 100*N |
*> iterations (in inner while loop). On exit Z holds |
*> iterations (in inner while loop). On exit Z holds |
*> a qd array with the same eigenvalues as the given Z. |
*> a qd array with the same eigenvalues as the given Z. |
*> = 3, termination criterion of outer while loop not met |
*> = 3, termination criterion of outer while loop not met |
*> (program created more than N unreduced blocks) |
*> (program created more than N unreduced blocks) |
*> \endverbatim |
*> \endverbatim |
* |
* |
* Authors: |
* Authors: |
* ======== |
* ======== |
* |
* |
*> \author Univ. of Tennessee |
*> \author Univ. of Tennessee |
*> \author Univ. of California Berkeley |
*> \author Univ. of California Berkeley |
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date September 2012 |
*> \date December 2016 |
* |
* |
*> \ingroup auxOTHERcomputational |
*> \ingroup auxOTHERcomputational |
* |
* |
Line 112
|
Line 112
|
* ===================================================================== |
* ===================================================================== |
SUBROUTINE DLASQ2( N, Z, INFO ) |
SUBROUTINE DLASQ2( N, Z, INFO ) |
* |
* |
* -- LAPACK computational routine (version 3.4.2) -- |
* -- 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..-- |
* September 2012 |
* December 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER INFO, N |
INTEGER INFO, N |
Line 136
|
Line 136
|
* .. Local Scalars .. |
* .. Local Scalars .. |
LOGICAL IEEE |
LOGICAL IEEE |
INTEGER I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, |
INTEGER I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, |
$ K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT, |
$ K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT, |
$ TTYPE |
$ 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, |
Line 155
|
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 195
|
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 264
|
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 333
|
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 364
|
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 386
|
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 404
|
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 421
|
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 446
|
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 = 100*( 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 497
|
Line 497
|
140 CONTINUE |
140 CONTINUE |
* |
* |
INFO = 2 |
INFO = 2 |
* |
* |
* Maximum number of iterations exceeded, restore the shift |
* Maximum number of iterations exceeded, restore the shift |
* SIGMA and place the new d's and e's in a qd array. |
* SIGMA and place the new d's and e's in a qd array. |
* This might need to be done for several blocks |
* This might need to be done for several blocks |
* |
* |
Line 549
|
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 570
|
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 ) |