LAPACK: SRC/slasq2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SLASQ2( N, Z, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2)                                    --
00004 *
00005 *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
00006 *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
00007 *  -- Berkeley                                                        --
00008 *  -- November 2008                                                   --
00009 *
00010 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00011 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00012 *
00013 *     .. Scalar Arguments ..
00014       INTEGER            INFO, N
00015 *     ..
00016 *     .. Array Arguments ..
00017       REAL               Z( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  SLASQ2 computes all the eigenvalues of the symmetric positive 
00024 *  definite tridiagonal matrix associated with the qd array Z to high
00025 *  relative accuracy are computed to high relative accuracy, in the
00026 *  absence of denormalization, underflow and overflow.
00027 *
00028 *  To see the relation of Z to the tridiagonal matrix, let L be a
00029 *  unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
00030 *  let U be an upper bidiagonal matrix with 1's above and diagonal
00031 *  Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
00032 *  symmetric tridiagonal to which it is similar.
00033 *
00034 *  Note : SLASQ2 defines a logical variable, IEEE, which is true
00035 *  on machines which follow ieee-754 floating-point standard in their
00036 *  handling of infinities and NaNs, and false otherwise. This variable
00037 *  is passed to SLASQ3.
00038 *
00039 *  Arguments
00040 *  =========
00041 *
00042 *  N     (input) INTEGER
00043 *        The number of rows and columns in the matrix. N >= 0.
00044 *
00045 *  Z     (input/output) REAL array, dimension ( 4*N )
00046 *        On entry Z holds the qd array. On exit, entries 1 to N hold
00047 *        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
00048 *        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
00049 *        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
00050 *        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
00051 *        shifts that failed.
00052 *
00053 *  INFO  (output) INTEGER
00054 *        = 0: successful exit
00055 *        < 0: if the i-th argument is a scalar and had an illegal
00056 *             value, then INFO = -i, if the i-th argument is an
00057 *             array and the j-entry had an illegal value, then
00058 *             INFO = -(i*100+j)
00059 *        > 0: the algorithm failed
00060 *              = 1, a split was marked by a positive value in E
00061 *              = 2, current block of Z not diagonalized after 30*N
00062 *                   iterations (in inner while loop)
00063 *              = 3, termination criterion of outer while loop not met 
00064 *                   (program created more than N unreduced blocks)
00065 *
00066 *  Further Details
00067 *  ===============
00068 *  Local Variables: I0:N0 defines a current unreduced segment of Z.
00069 *  The shifts are accumulated in SIGMA. Iteration count is in ITER.
00070 *  Ping-pong is controlled by PP (alternates between 0 and 1).
00071 *
00072 *  =====================================================================
00073 *
00074 *     .. Parameters ..
00075       REAL               CBIAS
00076       PARAMETER          ( CBIAS = 1.50E0 )
00077       REAL               ZERO, HALF, ONE, TWO, FOUR, HUNDRD
00078       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0,
00079      $                     TWO = 2.0E0, FOUR = 4.0E0, HUNDRD = 100.0E0 )
00080 *     ..
00081 *     .. Local Scalars ..
00082       LOGICAL            IEEE
00083       INTEGER            I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K,
00084      $                   KMIN, N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE
00085       REAL               D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN,
00086      $                   DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX,
00087      $                   QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL,
00088      $                   TOL2, TRACE, ZMAX
00089 *     ..
00090 *     .. External Subroutines ..
00091       EXTERNAL           SLASQ3, SLASRT, XERBLA
00092 *     ..
00093 *     .. External Functions ..
00094       INTEGER            ILAENV
00095       REAL               SLAMCH
00096       EXTERNAL           ILAENV, SLAMCH
00097 *     ..
00098 *     .. Intrinsic Functions ..
00099       INTRINSIC          ABS, MAX, MIN, REAL, SQRT
00100 *     ..
00101 *     .. Executable Statements ..
00102 *      
00103 *     Test the input arguments.
00104 *     (in case SLASQ2 is not called by SLASQ1)
00105 *
00106       INFO = 0
00107       EPS = SLAMCH( 'Precision' )
00108       SAFMIN = SLAMCH( 'Safe minimum' )
00109       TOL = EPS*HUNDRD
00110       TOL2 = TOL**2
00111 *
00112       IF( N.LT.0 ) THEN
00113          INFO = -1
00114          CALL XERBLA( 'SLASQ2', 1 )
00115          RETURN
00116       ELSE IF( N.EQ.0 ) THEN
00117          RETURN
00118       ELSE IF( N.EQ.1 ) THEN
00119 *
00120 *        1-by-1 case.
00121 *
00122          IF( Z( 1 ).LT.ZERO ) THEN
00123             INFO = -201
00124             CALL XERBLA( 'SLASQ2', 2 )
00125          END IF
00126          RETURN
00127       ELSE IF( N.EQ.2 ) THEN
00128 *
00129 *        2-by-2 case.
00130 *
00131          IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN
00132             INFO = -2
00133             CALL XERBLA( 'SLASQ2', 2 )
00134             RETURN
00135          ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
00136             D = Z( 3 )
00137             Z( 3 ) = Z( 1 )
00138             Z( 1 ) = D
00139          END IF
00140          Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
00141          IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
00142             T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) 
00143             S = Z( 3 )*( Z( 2 ) / T )
00144             IF( S.LE.T ) THEN
00145                S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
00146             ELSE
00147                S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
00148             END IF
00149             T = Z( 1 ) + ( S+Z( 2 ) )
00150             Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
00151             Z( 1 ) = T
00152          END IF
00153          Z( 2 ) = Z( 3 )
00154          Z( 6 ) = Z( 2 ) + Z( 1 )
00155          RETURN
00156       END IF
00157 *
00158 *     Check for negative data and compute sums of q's and e's.
00159 *
00160       Z( 2*N ) = ZERO
00161       EMIN = Z( 2 )
00162       QMAX = ZERO
00163       ZMAX = ZERO
00164       D = ZERO
00165       E = ZERO
00166 *
00167       DO 10 K = 1, 2*( N-1 ), 2
00168          IF( Z( K ).LT.ZERO ) THEN
00169             INFO = -( 200+K )
00170             CALL XERBLA( 'SLASQ2', 2 )
00171             RETURN
00172          ELSE IF( Z( K+1 ).LT.ZERO ) THEN
00173             INFO = -( 200+K+1 )
00174             CALL XERBLA( 'SLASQ2', 2 )
00175             RETURN
00176          END IF
00177          D = D + Z( K )
00178          E = E + Z( K+1 )
00179          QMAX = MAX( QMAX, Z( K ) )
00180          EMIN = MIN( EMIN, Z( K+1 ) )
00181          ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
00182    10 CONTINUE
00183       IF( Z( 2*N-1 ).LT.ZERO ) THEN
00184          INFO = -( 200+2*N-1 )
00185          CALL XERBLA( 'SLASQ2', 2 )
00186          RETURN
00187       END IF
00188       D = D + Z( 2*N-1 )
00189       QMAX = MAX( QMAX, Z( 2*N-1 ) )
00190       ZMAX = MAX( QMAX, ZMAX )
00191 *
00192 *     Check for diagonality.
00193 *
00194       IF( E.EQ.ZERO ) THEN
00195          DO 20 K = 2, N
00196             Z( K ) = Z( 2*K-1 )
00197    20    CONTINUE
00198          CALL SLASRT( 'D', N, Z, IINFO )
00199          Z( 2*N-1 ) = D
00200          RETURN
00201       END IF
00202 *
00203       TRACE = D + E
00204 *
00205 *     Check for zero data.
00206 *
00207       IF( TRACE.EQ.ZERO ) THEN
00208          Z( 2*N-1 ) = ZERO
00209          RETURN
00210       END IF
00211 *         
00212 *     Check whether the machine is IEEE conformable.
00213 *         
00214 *     IEEE = ILAENV( 10, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
00215 *    $       ILAENV( 11, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1      
00216 *
00217 *     [11/15/2008] The case IEEE=.TRUE. has a problem in single precision with
00218 *     some the test matrices of type 16. The double precision code is fine.
00219 *
00220       IEEE = .FALSE.
00221 *         
00222 *     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
00223 *
00224       DO 30 K = 2*N, 2, -2
00225          Z( 2*K ) = ZERO 
00226          Z( 2*K-1 ) = Z( K ) 
00227          Z( 2*K-2 ) = ZERO 
00228          Z( 2*K-3 ) = Z( K-1 ) 
00229    30 CONTINUE
00230 *
00231       I0 = 1
00232       N0 = N
00233 *
00234 *     Reverse the qd-array, if warranted.
00235 *
00236       IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN
00237          IPN4 = 4*( I0+N0 )
00238          DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4
00239             TEMP = Z( I4-3 )
00240             Z( I4-3 ) = Z( IPN4-I4-3 )
00241             Z( IPN4-I4-3 ) = TEMP
00242             TEMP = Z( I4-1 )
00243             Z( I4-1 ) = Z( IPN4-I4-5 )
00244             Z( IPN4-I4-5 ) = TEMP
00245    40    CONTINUE
00246       END IF
00247 *
00248 *     Initial split checking via dqd and Li's test.
00249 *
00250       PP = 0
00251 *
00252       DO 80 K = 1, 2
00253 *
00254          D = Z( 4*N0+PP-3 )
00255          DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4
00256             IF( Z( I4-1 ).LE.TOL2*D ) THEN
00257                Z( I4-1 ) = -ZERO
00258                D = Z( I4-3 )
00259             ELSE
00260                D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) )
00261             END IF
00262    50    CONTINUE
00263 *
00264 *        dqd maps Z to ZZ plus Li's test.
00265 *
00266          EMIN = Z( 4*I0+PP+1 )
00267          D = Z( 4*I0+PP-3 )
00268          DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4
00269             Z( I4-2*PP-2 ) = D + Z( I4-1 )
00270             IF( Z( I4-1 ).LE.TOL2*D ) THEN
00271                Z( I4-1 ) = -ZERO
00272                Z( I4-2*PP-2 ) = D
00273                Z( I4-2*PP ) = ZERO
00274                D = Z( I4+1 )
00275             ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND.
00276      $               SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN
00277                TEMP = Z( I4+1 ) / Z( I4-2*PP-2 )
00278                Z( I4-2*PP ) = Z( I4-1 )*TEMP
00279                D = D*TEMP
00280             ELSE
00281                Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) )
00282                D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )
00283             END IF
00284             EMIN = MIN( EMIN, Z( I4-2*PP ) )
00285    60    CONTINUE 
00286          Z( 4*N0-PP-2 ) = D
00287 *
00288 *        Now find qmax.
00289 *
00290          QMAX = Z( 4*I0-PP-2 )
00291          DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4
00292             QMAX = MAX( QMAX, Z( I4 ) )
00293    70    CONTINUE
00294 *
00295 *        Prepare for the next iteration on K.
00296 *
00297          PP = 1 - PP
00298    80 CONTINUE
00299 *
00300 *     Initialise variables to pass to SLASQ3.
00301 *
00302       TTYPE = 0
00303       DMIN1 = ZERO
00304       DMIN2 = ZERO
00305       DN    = ZERO
00306       DN1   = ZERO
00307       DN2   = ZERO
00308       G     = ZERO
00309       TAU   = ZERO
00310 *
00311       ITER = 2
00312       NFAIL = 0
00313       NDIV = 2*( N0-I0 )
00314 *
00315       DO 160 IWHILA = 1, N + 1
00316          IF( N0.LT.1 ) 
00317      $      GO TO 170
00318 *
00319 *        While array unfinished do 
00320 *
00321 *        E(N0) holds the value of SIGMA when submatrix in I0:N0
00322 *        splits from the rest of the array, but is negated.
00323 *      
00324          DESIG = ZERO
00325          IF( N0.EQ.N ) THEN
00326             SIGMA = ZERO
00327          ELSE
00328             SIGMA = -Z( 4*N0-1 )
00329          END IF
00330          IF( SIGMA.LT.ZERO ) THEN
00331             INFO = 1
00332             RETURN
00333          END IF
00334 *
00335 *        Find last unreduced submatrix's top index I0, find QMAX and
00336 *        EMIN. Find Gershgorin-type bound if Q's much greater than E's.
00337 *
00338          EMAX = ZERO 
00339          IF( N0.GT.I0 ) THEN
00340             EMIN = ABS( Z( 4*N0-5 ) )
00341          ELSE
00342             EMIN = ZERO
00343          END IF
00344          QMIN = Z( 4*N0-3 )
00345          QMAX = QMIN
00346          DO 90 I4 = 4*N0, 8, -4
00347             IF( Z( I4-5 ).LE.ZERO )
00348      $         GO TO 100
00349             IF( QMIN.GE.FOUR*EMAX ) THEN
00350                QMIN = MIN( QMIN, Z( I4-3 ) )
00351                EMAX = MAX( EMAX, Z( I4-5 ) )
00352             END IF
00353             QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )
00354             EMIN = MIN( EMIN, Z( I4-5 ) )
00355    90    CONTINUE
00356          I4 = 4 
00357 *
00358   100    CONTINUE
00359          I0 = I4 / 4
00360          PP = 0
00361 *
00362          IF( N0-I0.GT.1 ) THEN
00363             DEE = Z( 4*I0-3 )
00364             DEEMIN = DEE
00365             KMIN = I0
00366             DO 110 I4 = 4*I0+1, 4*N0-3, 4
00367                DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) )
00368                IF( DEE.LE.DEEMIN ) THEN
00369                   DEEMIN = DEE
00370                   KMIN = ( I4+3 )/4
00371                END IF
00372   110       CONTINUE
00373             IF( (KMIN-I0)*2.LT.N0-KMIN .AND. 
00374      $         DEEMIN.LE.HALF*Z(4*N0-3) ) THEN
00375                IPN4 = 4*( I0+N0 )
00376                PP = 2
00377                DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4
00378                   TEMP = Z( I4-3 )
00379                   Z( I4-3 ) = Z( IPN4-I4-3 )
00380                   Z( IPN4-I4-3 ) = TEMP
00381                   TEMP = Z( I4-2 )
00382                   Z( I4-2 ) = Z( IPN4-I4-2 )
00383                   Z( IPN4-I4-2 ) = TEMP
00384                   TEMP = Z( I4-1 )
00385                   Z( I4-1 ) = Z( IPN4-I4-5 )
00386                   Z( IPN4-I4-5 ) = TEMP
00387                   TEMP = Z( I4 )
00388                   Z( I4 ) = Z( IPN4-I4-4 )
00389                   Z( IPN4-I4-4 ) = TEMP
00390   120          CONTINUE
00391             END IF
00392          END IF
00393 *
00394 *        Put -(initial shift) into DMIN.
00395 *
00396          DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )
00397 *
00398 *        Now I0:N0 is unreduced. 
00399 *        PP = 0 for ping, PP = 1 for pong.
00400 *        PP = 2 indicates that flipping was applied to the Z array and
00401 *               and that the tests for deflation upon entry in SLASQ3 
00402 *               should not be performed.
00403 *
00404          NBIG = 30*( N0-I0+1 )
00405          DO 140 IWHILB = 1, NBIG
00406             IF( I0.GT.N0 ) 
00407      $         GO TO 150
00408 *
00409 *           While submatrix unfinished take a good dqds step.
00410 *
00411             CALL SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
00412      $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
00413      $                   DN2, G, TAU )
00414 *
00415             PP = 1 - PP
00416 *
00417 *           When EMIN is very small check for splits.
00418 *
00419             IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN
00420                IF( Z( 4*N0 ).LE.TOL2*QMAX .OR.
00421      $             Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN
00422                   SPLT = I0 - 1
00423                   QMAX = Z( 4*I0-3 )
00424                   EMIN = Z( 4*I0-1 )
00425                   OLDEMN = Z( 4*I0 )
00426                   DO 130 I4 = 4*I0, 4*( N0-3 ), 4
00427                      IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR.
00428      $                   Z( I4-1 ).LE.TOL2*SIGMA ) THEN
00429                         Z( I4-1 ) = -SIGMA
00430                         SPLT = I4 / 4
00431                         QMAX = ZERO
00432                         EMIN = Z( I4+3 )
00433                         OLDEMN = Z( I4+4 )
00434                      ELSE
00435                         QMAX = MAX( QMAX, Z( I4+1 ) )
00436                         EMIN = MIN( EMIN, Z( I4-1 ) )
00437                         OLDEMN = MIN( OLDEMN, Z( I4 ) )
00438                      END IF
00439   130             CONTINUE
00440                   Z( 4*N0-1 ) = EMIN
00441                   Z( 4*N0 ) = OLDEMN
00442                   I0 = SPLT + 1
00443                END IF
00444             END IF
00445 *
00446   140    CONTINUE
00447 *
00448          INFO = 2
00449          RETURN
00450 *
00451 *        end IWHILB
00452 *
00453   150    CONTINUE
00454 *
00455   160 CONTINUE
00456 *
00457       INFO = 3
00458       RETURN
00459 *
00460 *     end IWHILA   
00461 *
00462   170 CONTINUE
00463 *      
00464 *     Move q's to the front.
00465 *      
00466       DO 180 K = 2, N
00467          Z( K ) = Z( 4*K-3 )
00468   180 CONTINUE
00469 *      
00470 *     Sort and compute sum of eigenvalues.
00471 *
00472       CALL SLASRT( 'D', N, Z, IINFO )
00473 *
00474       E = ZERO
00475       DO 190 K = N, 1, -1
00476          E = E + Z( K )
00477   190 CONTINUE
00478 *
00479 *     Store trace, sum(eigenvalues) and information on performance.
00480 *
00481       Z( 2*N+1 ) = TRACE 
00482       Z( 2*N+2 ) = E
00483       Z( 2*N+3 ) = REAL( ITER )
00484       Z( 2*N+4 ) = REAL( NDIV ) / REAL( N**2 )
00485       Z( 2*N+5 ) = HUNDRD*NFAIL / REAL( ITER )
00486       RETURN
00487 *
00488 *     End of SLASQ2
00489 *
00490       END