diff options
Diffstat (limited to '2.3-1/src/fortran/lapack/dlasq2.f')
| -rw-r--r-- | 2.3-1/src/fortran/lapack/dlasq2.f | 448 |
1 files changed, 448 insertions, 0 deletions
diff --git a/2.3-1/src/fortran/lapack/dlasq2.f b/2.3-1/src/fortran/lapack/dlasq2.f new file mode 100644 index 00000000..b6b79aeb --- /dev/null +++ b/2.3-1/src/fortran/lapack/dlasq2.f @@ -0,0 +1,448 @@ + SUBROUTINE DLASQ2( N, Z, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* Modified to call DLAZQ3 in place of DLASQ3, 13 Feb 03, SJH. +* +* .. Scalar Arguments .. + INTEGER INFO, N +* .. +* .. Array Arguments .. + 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 DLAZQ3. +* +* Arguments +* ========= +* +* N (input) INTEGER +* The number of rows and columns in the matrix. N >= 0. +* +* Z (workspace) 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 .. + 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, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K, + $ N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE + DOUBLE PRECISION D, DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, E, + $ EMAX, EMIN, EPS, OLDEMN, QMAX, QMIN, S, SAFMIN, + $ SIGMA, T, TAU, TEMP, TOL, TOL2, TRACE, ZMAX +* .. +* .. External Subroutines .. + EXTERNAL DLAZQ3, 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 DLAZQ3 +* + TTYPE = 0 + DMIN1 = ZERO + DMIN2 = ZERO + DN = ZERO + DN1 = ZERO + DN2 = ZERO + TAU = ZERO +* + ITER = 2 + NFAIL = 0 + NDIV = 2*( N0-I0 ) +* + DO 140 IWHILA = 1, N + 1 + IF( N0.LT.1 ) + $ GO TO 150 +* +* 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 +* +* Store EMIN for passing to DLAZQ3. +* + Z( 4*N0-1 ) = EMIN +* +* 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 = 0 +* + NBIG = 30*( N0-I0+1 ) + DO 120 IWHILB = 1, NBIG + IF( I0.GT.N0 ) + $ GO TO 130 +* +* While submatrix unfinished take a good dqds step. +* + CALL DLAZQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, + $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, + $ DN2, 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 110 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 + 110 CONTINUE + Z( 4*N0-1 ) = EMIN + Z( 4*N0 ) = OLDEMN + I0 = SPLT + 1 + END IF + END IF +* + 120 CONTINUE +* + INFO = 2 + RETURN +* +* end IWHILB +* + 130 CONTINUE +* + 140 CONTINUE +* + INFO = 3 + RETURN +* +* end IWHILA +* + 150 CONTINUE +* +* Move q's to the front. +* + DO 160 K = 2, N + Z( K ) = Z( 4*K-3 ) + 160 CONTINUE +* +* Sort and compute sum of eigenvalues. +* + CALL DLASRT( 'D', N, Z, IINFO ) +* + E = ZERO + DO 170 K = N, 1, -1 + E = E + Z( K ) + 170 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 |