#include "ESMF_LapackBlas.inc" !> \brief \b DLASQ2 ! ! =========== DOCUMENTATION =========== ! ! Online html documentation available at ! http://www.netlib.org/lapack/explore-html/ ! !> \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 November 2011 ! !> \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.4.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! November 2011 ! ! .. Scalar Arguments .. INTEGER INFO, N ! .. ! .. Array Arguments .. DOUBLE PRECISION Z( * ) ! .. ! ! ===================================================================== ! ! .. Parameters .. DOUBLE PRECISION CBIAS PARAMETER ( CBIAS = 1.50D0 ) DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, HUNDRD PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, & & TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 ) ! .. ! .. Local Scalars .. LOGICAL IEEE INTEGER I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, & & K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT, & & TTYPE DOUBLE PRECISION D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN, & & DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX, & & QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL, & & TOL2, TRACE, ZMAX, TEMPE, TEMPQ ! .. ! .. External Subroutines .. EXTERNAL DLASQ3, DLASRT, XERBLA ! .. ! .. External Functions .. INTEGER ILAENV DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH, ILAENV ! .. ! .. Intrinsic Functions .. INTRINSIC ABS, DBLE, MAX, MIN, SQRT ! .. ! .. Executable Statements .. ! ! Test the input arguments. ! (in case DLASQ2 is not called by DLASQ1) ! INFO = 0 EPS = DLAMCH( 'Precision' ) SAFMIN = DLAMCH( 'Safe minimum' ) TOL = EPS*HUNDRD TOL2 = TOL**2 ! IF( N.LT.0 ) THEN INFO = -1 CALL XERBLA( 'DLASQ2', 1 ) RETURN ELSE IF( N.EQ.0 ) THEN RETURN ELSE IF( N.EQ.1 ) THEN ! ! 1-by-1 case. ! IF( Z( 1 ).LT.ZERO ) THEN INFO = -201 CALL XERBLA( 'DLASQ2', 2 ) END IF RETURN ELSE IF( N.EQ.2 ) THEN ! ! 2-by-2 case. ! IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN INFO = -2 CALL XERBLA( 'DLASQ2', 2 ) RETURN ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN D = Z( 3 ) Z( 3 ) = Z( 1 ) Z( 1 ) = D END IF Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 ) IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) S = Z( 3 )*( Z( 2 ) / T ) IF( S.LE.T ) THEN S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) ) ELSE S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) ) END IF T = Z( 1 ) + ( S+Z( 2 ) ) Z( 3 ) = Z( 3 )*( Z( 1 ) / T ) Z( 1 ) = T END IF Z( 2 ) = Z( 3 ) Z( 6 ) = Z( 2 ) + Z( 1 ) RETURN END IF ! ! Check for negative data and compute sums of q's and e's. ! Z( 2*N ) = ZERO EMIN = Z( 2 ) QMAX = ZERO ZMAX = ZERO D = ZERO E = ZERO ! DO 10 K = 1, 2*( N-1 ), 2 IF( Z( K ).LT.ZERO ) THEN INFO = -( 200+K ) CALL XERBLA( 'DLASQ2', 2 ) RETURN ELSE IF( Z( K+1 ).LT.ZERO ) THEN INFO = -( 200+K+1 ) CALL XERBLA( 'DLASQ2', 2 ) RETURN END IF D = D + Z( K ) E = E + Z( K+1 ) QMAX = MAX( QMAX, Z( K ) ) EMIN = MIN( EMIN, Z( K+1 ) ) ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) ) 10 CONTINUE IF( Z( 2*N-1 ).LT.ZERO ) THEN INFO = -( 200+2*N-1 ) CALL XERBLA( 'DLASQ2', 2 ) RETURN END IF D = D + Z( 2*N-1 ) QMAX = MAX( QMAX, Z( 2*N-1 ) ) ZMAX = MAX( QMAX, ZMAX ) ! ! Check for diagonality. ! IF( E.EQ.ZERO ) THEN DO 20 K = 2, N Z( K ) = Z( 2*K-1 ) 20 CONTINUE CALL DLASRT( 'D', N, Z, IINFO ) Z( 2*N-1 ) = D RETURN END IF ! TRACE = D + E ! ! Check for zero data. ! IF( TRACE.EQ.ZERO ) THEN Z( 2*N-1 ) = ZERO RETURN END IF ! ! Check whether the machine is IEEE conformable. ! IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND. & & ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 ! ! Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). ! DO 30 K = 2*N, 2, -2 Z( 2*K ) = ZERO Z( 2*K-1 ) = Z( K ) Z( 2*K-2 ) = ZERO Z( 2*K-3 ) = Z( K-1 ) 30 CONTINUE ! I0 = 1 N0 = N ! ! Reverse the qd-array, if warranted. ! IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN IPN4 = 4*( I0+N0 ) DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4 TEMP = Z( I4-3 ) Z( I4-3 ) = Z( IPN4-I4-3 ) Z( IPN4-I4-3 ) = TEMP TEMP = Z( I4-1 ) Z( I4-1 ) = Z( IPN4-I4-5 ) Z( IPN4-I4-5 ) = TEMP 40 CONTINUE END IF ! ! Initial split checking via dqd and Li's test. ! PP = 0 ! DO 80 K = 1, 2 ! D = Z( 4*N0+PP-3 ) DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4 IF( Z( I4-1 ).LE.TOL2*D ) THEN Z( I4-1 ) = -ZERO D = Z( I4-3 ) ELSE D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) ) END IF 50 CONTINUE ! ! dqd maps Z to ZZ plus Li's test. ! EMIN = Z( 4*I0+PP+1 ) D = Z( 4*I0+PP-3 ) DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4 Z( I4-2*PP-2 ) = D + Z( I4-1 ) IF( Z( I4-1 ).LE.TOL2*D ) THEN Z( I4-1 ) = -ZERO Z( I4-2*PP-2 ) = D Z( I4-2*PP ) = ZERO D = Z( I4+1 ) ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND. & & SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN TEMP = Z( I4+1 ) / Z( I4-2*PP-2 ) Z( I4-2*PP ) = Z( I4-1 )*TEMP D = D*TEMP ELSE Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) ) D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) ) END IF EMIN = MIN( EMIN, Z( I4-2*PP ) ) 60 CONTINUE Z( 4*N0-PP-2 ) = D ! ! Now find qmax. ! QMAX = Z( 4*I0-PP-2 ) DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4 QMAX = MAX( QMAX, Z( I4 ) ) 70 CONTINUE ! ! Prepare for the next iteration on K. ! PP = 1 - PP 80 CONTINUE ! ! Initialise variables to pass to DLASQ3. ! TTYPE = 0 DMIN1 = ZERO DMIN2 = ZERO DN = ZERO DN1 = ZERO DN2 = ZERO G = ZERO TAU = ZERO ! ITER = 2 NFAIL = 0 NDIV = 2*( N0-I0 ) ! DO 160 IWHILA = 1, N + 1 IF( N0.LT.1 ) & & GO TO 170 ! ! While array unfinished do ! ! E(N0) holds the value of SIGMA when submatrix in I0:N0 ! splits from the rest of the array, but is negated. ! DESIG = ZERO IF( N0.EQ.N ) THEN SIGMA = ZERO ELSE SIGMA = -Z( 4*N0-1 ) END IF IF( SIGMA.LT.ZERO ) THEN INFO = 1 RETURN END IF ! ! Find last unreduced submatrix's top index I0, find QMAX and ! EMIN. Find Gershgorin-type bound if Q's much greater than E's. ! EMAX = ZERO IF( N0.GT.I0 ) THEN EMIN = ABS( Z( 4*N0-5 ) ) ELSE EMIN = ZERO END IF QMIN = Z( 4*N0-3 ) QMAX = QMIN DO 90 I4 = 4*N0, 8, -4 IF( Z( I4-5 ).LE.ZERO ) & & GO TO 100 IF( QMIN.GE.FOUR*EMAX ) THEN QMIN = MIN( QMIN, Z( I4-3 ) ) EMAX = MAX( EMAX, Z( I4-5 ) ) END IF QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) ) EMIN = MIN( EMIN, Z( I4-5 ) ) 90 CONTINUE I4 = 4 ! 100 CONTINUE I0 = I4 / 4 PP = 0 ! IF( N0-I0.GT.1 ) THEN DEE = Z( 4*I0-3 ) DEEMIN = DEE KMIN = I0 DO 110 I4 = 4*I0+1, 4*N0-3, 4 DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) ) IF( DEE.LE.DEEMIN ) THEN DEEMIN = DEE KMIN = ( I4+3 )/4 END IF 110 CONTINUE IF( (KMIN-I0)*2.LT.N0-KMIN .AND. & & DEEMIN.LE.HALF*Z(4*N0-3) ) THEN IPN4 = 4*( I0+N0 ) PP = 2 DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4 TEMP = Z( I4-3 ) Z( I4-3 ) = Z( IPN4-I4-3 ) Z( IPN4-I4-3 ) = TEMP TEMP = Z( I4-2 ) Z( I4-2 ) = Z( IPN4-I4-2 ) Z( IPN4-I4-2 ) = TEMP TEMP = Z( I4-1 ) Z( I4-1 ) = Z( IPN4-I4-5 ) Z( IPN4-I4-5 ) = TEMP TEMP = Z( I4 ) Z( I4 ) = Z( IPN4-I4-4 ) Z( IPN4-I4-4 ) = TEMP 120 CONTINUE END IF END IF ! ! Put -(initial shift) into DMIN. ! DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) ) ! ! Now I0:N0 is unreduced. ! PP = 0 for ping, PP = 1 for pong. ! PP = 2 indicates that flipping was applied to the Z array and ! and that the tests for deflation upon entry in DLASQ3 ! should not be performed. ! NBIG = 100*( N0-I0+1 ) DO 140 IWHILB = 1, NBIG IF( I0.GT.N0 ) & & GO TO 150 ! ! While submatrix unfinished take a good dqds step. ! CALL DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, & & ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, & & DN2, G, TAU ) ! PP = 1 - PP ! ! When EMIN is very small check for splits. ! IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN IF( Z( 4*N0 ).LE.TOL2*QMAX .OR. & & Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN SPLT = I0 - 1 QMAX = Z( 4*I0-3 ) EMIN = Z( 4*I0-1 ) OLDEMN = Z( 4*I0 ) DO 130 I4 = 4*I0, 4*( N0-3 ), 4 IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR. & & Z( I4-1 ).LE.TOL2*SIGMA ) THEN Z( I4-1 ) = -SIGMA SPLT = I4 / 4 QMAX = ZERO EMIN = Z( I4+3 ) OLDEMN = Z( I4+4 ) ELSE QMAX = MAX( QMAX, Z( I4+1 ) ) EMIN = MIN( EMIN, Z( I4-1 ) ) OLDEMN = MIN( OLDEMN, Z( I4 ) ) END IF 130 CONTINUE Z( 4*N0-1 ) = EMIN Z( 4*N0 ) = OLDEMN I0 = SPLT + 1 END IF END IF ! 140 CONTINUE ! 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 ! ! end IWHILB ! 150 CONTINUE ! 160 CONTINUE ! INFO = 3 RETURN ! ! end IWHILA ! 170 CONTINUE ! ! Move q's to the front. ! DO 180 K = 2, N Z( K ) = Z( 4*K-3 ) 180 CONTINUE ! ! Sort and compute sum of eigenvalues. ! CALL DLASRT( 'D', N, Z, IINFO ) ! E = ZERO DO 190 K = N, 1, -1 E = E + Z( K ) 190 CONTINUE ! ! Store trace, sum(eigenvalues) and information on performance. ! Z( 2*N+1 ) = TRACE Z( 2*N+2 ) = E Z( 2*N+3 ) = DBLE( ITER ) Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 ) Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER ) RETURN ! ! End of DLASQ2 ! END