LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlasq4.f
Go to the documentation of this file.
00001 *> \brief \b DLASQ4
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLASQ4 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq4.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq4.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq4.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
00022 *                          DN1, DN2, TAU, TTYPE, G )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            I0, N0, N0IN, PP, TTYPE
00026 *       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   Z( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> DLASQ4 computes an approximation TAU to the smallest eigenvalue
00039 *> using values of d from the previous transform.
00040 *> \endverbatim
00041 *
00042 *  Arguments:
00043 *  ==========
00044 *
00045 *> \param[in] I0
00046 *> \verbatim
00047 *>          I0 is INTEGER
00048 *>        First index.
00049 *> \endverbatim
00050 *>
00051 *> \param[in] N0
00052 *> \verbatim
00053 *>          N0 is INTEGER
00054 *>        Last index.
00055 *> \endverbatim
00056 *>
00057 *> \param[in] Z
00058 *> \verbatim
00059 *>          Z is DOUBLE PRECISION array, dimension ( 4*N )
00060 *>        Z holds the qd array.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] PP
00064 *> \verbatim
00065 *>          PP is INTEGER
00066 *>        PP=0 for ping, PP=1 for pong.
00067 *> \endverbatim
00068 *>
00069 *> \param[in] N0IN
00070 *> \verbatim
00071 *>          N0IN is INTEGER
00072 *>        The value of N0 at start of EIGTEST.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] DMIN
00076 *> \verbatim
00077 *>          DMIN is DOUBLE PRECISION
00078 *>        Minimum value of d.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] DMIN1
00082 *> \verbatim
00083 *>          DMIN1 is DOUBLE PRECISION
00084 *>        Minimum value of d, excluding D( N0 ).
00085 *> \endverbatim
00086 *>
00087 *> \param[in] DMIN2
00088 *> \verbatim
00089 *>          DMIN2 is DOUBLE PRECISION
00090 *>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
00091 *> \endverbatim
00092 *>
00093 *> \param[in] DN
00094 *> \verbatim
00095 *>          DN is DOUBLE PRECISION
00096 *>        d(N)
00097 *> \endverbatim
00098 *>
00099 *> \param[in] DN1
00100 *> \verbatim
00101 *>          DN1 is DOUBLE PRECISION
00102 *>        d(N-1)
00103 *> \endverbatim
00104 *>
00105 *> \param[in] DN2
00106 *> \verbatim
00107 *>          DN2 is DOUBLE PRECISION
00108 *>        d(N-2)
00109 *> \endverbatim
00110 *>
00111 *> \param[out] TAU
00112 *> \verbatim
00113 *>          TAU is DOUBLE PRECISION
00114 *>        This is the shift.
00115 *> \endverbatim
00116 *>
00117 *> \param[out] TTYPE
00118 *> \verbatim
00119 *>          TTYPE is INTEGER
00120 *>        Shift type.
00121 *> \endverbatim
00122 *>
00123 *> \param[in,out] G
00124 *> \verbatim
00125 *>          G is REAL
00126 *>        G is passed as an argument in order to save its value between
00127 *>        calls to DLASQ4.
00128 *> \endverbatim
00129 *
00130 *  Authors:
00131 *  ========
00132 *
00133 *> \author Univ. of Tennessee 
00134 *> \author Univ. of California Berkeley 
00135 *> \author Univ. of Colorado Denver 
00136 *> \author NAG Ltd. 
00137 *
00138 *> \date November 2011
00139 *
00140 *> \ingroup auxOTHERcomputational
00141 *
00142 *> \par Further Details:
00143 *  =====================
00144 *>
00145 *> \verbatim
00146 *>
00147 *>  CNST1 = 9/16
00148 *> \endverbatim
00149 *>
00150 *  =====================================================================
00151       SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
00152      $                   DN1, DN2, TAU, TTYPE, G )
00153 *
00154 *  -- LAPACK computational routine (version 3.4.0) --
00155 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00156 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00157 *     November 2011
00158 *
00159 *     .. Scalar Arguments ..
00160       INTEGER            I0, N0, N0IN, PP, TTYPE
00161       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
00162 *     ..
00163 *     .. Array Arguments ..
00164       DOUBLE PRECISION   Z( * )
00165 *     ..
00166 *
00167 *  =====================================================================
00168 *
00169 *     .. Parameters ..
00170       DOUBLE PRECISION   CNST1, CNST2, CNST3
00171       PARAMETER          ( CNST1 = 0.5630D0, CNST2 = 1.010D0,
00172      $                   CNST3 = 1.050D0 )
00173       DOUBLE PRECISION   QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
00174       PARAMETER          ( QURTR = 0.250D0, THIRD = 0.3330D0,
00175      $                   HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0,
00176      $                   TWO = 2.0D0, HUNDRD = 100.0D0 )
00177 *     ..
00178 *     .. Local Scalars ..
00179       INTEGER            I4, NN, NP
00180       DOUBLE PRECISION   A2, B1, B2, GAM, GAP1, GAP2, S
00181 *     ..
00182 *     .. Intrinsic Functions ..
00183       INTRINSIC          MAX, MIN, SQRT
00184 *     ..
00185 *     .. Executable Statements ..
00186 *
00187 *     A negative DMIN forces the shift to take that absolute value
00188 *     TTYPE records the type of shift.
00189 *
00190       IF( DMIN.LE.ZERO ) THEN
00191          TAU = -DMIN
00192          TTYPE = -1
00193          RETURN
00194       END IF
00195 *       
00196       NN = 4*N0 + PP
00197       IF( N0IN.EQ.N0 ) THEN
00198 *
00199 *        No eigenvalues deflated.
00200 *
00201          IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN
00202 *
00203             B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) )
00204             B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) )
00205             A2 = Z( NN-7 ) + Z( NN-5 )
00206 *
00207 *           Cases 2 and 3.
00208 *
00209             IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN
00210                GAP2 = DMIN2 - A2 - DMIN2*QURTR
00211                IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN
00212                   GAP1 = A2 - DN - ( B2 / GAP2 )*B2
00213                ELSE
00214                   GAP1 = A2 - DN - ( B1+B2 )
00215                END IF
00216                IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN
00217                   S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN )
00218                   TTYPE = -2
00219                ELSE
00220                   S = ZERO
00221                   IF( DN.GT.B1 )
00222      $               S = DN - B1
00223                   IF( A2.GT.( B1+B2 ) )
00224      $               S = MIN( S, A2-( B1+B2 ) )
00225                   S = MAX( S, THIRD*DMIN )
00226                   TTYPE = -3
00227                END IF
00228             ELSE
00229 *
00230 *              Case 4.
00231 *
00232                TTYPE = -4
00233                S = QURTR*DMIN
00234                IF( DMIN.EQ.DN ) THEN
00235                   GAM = DN
00236                   A2 = ZERO
00237                   IF( Z( NN-5 ) .GT. Z( NN-7 ) )
00238      $               RETURN
00239                   B2 = Z( NN-5 ) / Z( NN-7 )
00240                   NP = NN - 9
00241                ELSE
00242                   NP = NN - 2*PP
00243                   B2 = Z( NP-2 )
00244                   GAM = DN1
00245                   IF( Z( NP-4 ) .GT. Z( NP-2 ) )
00246      $               RETURN
00247                   A2 = Z( NP-4 ) / Z( NP-2 )
00248                   IF( Z( NN-9 ) .GT. Z( NN-11 ) )
00249      $               RETURN
00250                   B2 = Z( NN-9 ) / Z( NN-11 )
00251                   NP = NN - 13
00252                END IF
00253 *
00254 *              Approximate contribution to norm squared from I < NN-1.
00255 *
00256                A2 = A2 + B2
00257                DO 10 I4 = NP, 4*I0 - 1 + PP, -4
00258                   IF( B2.EQ.ZERO )
00259      $               GO TO 20
00260                   B1 = B2
00261                   IF( Z( I4 ) .GT. Z( I4-2 ) )
00262      $               RETURN
00263                   B2 = B2*( Z( I4 ) / Z( I4-2 ) )
00264                   A2 = A2 + B2
00265                   IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) 
00266      $               GO TO 20
00267    10          CONTINUE
00268    20          CONTINUE
00269                A2 = CNST3*A2
00270 *
00271 *              Rayleigh quotient residual bound.
00272 *
00273                IF( A2.LT.CNST1 )
00274      $            S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
00275             END IF
00276          ELSE IF( DMIN.EQ.DN2 ) THEN
00277 *
00278 *           Case 5.
00279 *
00280             TTYPE = -5
00281             S = QURTR*DMIN
00282 *
00283 *           Compute contribution to norm squared from I > NN-2.
00284 *
00285             NP = NN - 2*PP
00286             B1 = Z( NP-2 )
00287             B2 = Z( NP-6 )
00288             GAM = DN2
00289             IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
00290      $         RETURN
00291             A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
00292 *
00293 *           Approximate contribution to norm squared from I < NN-2.
00294 *
00295             IF( N0-I0.GT.2 ) THEN
00296                B2 = Z( NN-13 ) / Z( NN-15 )
00297                A2 = A2 + B2
00298                DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
00299                   IF( B2.EQ.ZERO )
00300      $               GO TO 40
00301                   B1 = B2
00302                   IF( Z( I4 ) .GT. Z( I4-2 ) )
00303      $               RETURN
00304                   B2 = B2*( Z( I4 ) / Z( I4-2 ) )
00305                   A2 = A2 + B2
00306                   IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) 
00307      $               GO TO 40
00308    30          CONTINUE
00309    40          CONTINUE
00310                A2 = CNST3*A2
00311             END IF
00312 *
00313             IF( A2.LT.CNST1 )
00314      $         S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
00315          ELSE
00316 *
00317 *           Case 6, no information to guide us.
00318 *
00319             IF( TTYPE.EQ.-6 ) THEN
00320                G = G + THIRD*( ONE-G )
00321             ELSE IF( TTYPE.EQ.-18 ) THEN
00322                G = QURTR*THIRD
00323             ELSE
00324                G = QURTR
00325             END IF
00326             S = G*DMIN
00327             TTYPE = -6
00328          END IF
00329 *
00330       ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
00331 *
00332 *        One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
00333 *
00334          IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN 
00335 *
00336 *           Cases 7 and 8.
00337 *
00338             TTYPE = -7
00339             S = THIRD*DMIN1
00340             IF( Z( NN-5 ).GT.Z( NN-7 ) )
00341      $         RETURN
00342             B1 = Z( NN-5 ) / Z( NN-7 )
00343             B2 = B1
00344             IF( B2.EQ.ZERO )
00345      $         GO TO 60
00346             DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
00347                A2 = B1
00348                IF( Z( I4 ).GT.Z( I4-2 ) )
00349      $            RETURN
00350                B1 = B1*( Z( I4 ) / Z( I4-2 ) )
00351                B2 = B2 + B1
00352                IF( HUNDRD*MAX( B1, A2 ).LT.B2 ) 
00353      $            GO TO 60
00354    50       CONTINUE
00355    60       CONTINUE
00356             B2 = SQRT( CNST3*B2 )
00357             A2 = DMIN1 / ( ONE+B2**2 )
00358             GAP2 = HALF*DMIN2 - A2
00359             IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
00360                S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
00361             ELSE 
00362                S = MAX( S, A2*( ONE-CNST2*B2 ) )
00363                TTYPE = -8
00364             END IF
00365          ELSE
00366 *
00367 *           Case 9.
00368 *
00369             S = QURTR*DMIN1
00370             IF( DMIN1.EQ.DN1 )
00371      $         S = HALF*DMIN1
00372             TTYPE = -9
00373          END IF
00374 *
00375       ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
00376 *
00377 *        Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
00378 *
00379 *        Cases 10 and 11.
00380 *
00381          IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN 
00382             TTYPE = -10
00383             S = THIRD*DMIN2
00384             IF( Z( NN-5 ).GT.Z( NN-7 ) )
00385      $         RETURN
00386             B1 = Z( NN-5 ) / Z( NN-7 )
00387             B2 = B1
00388             IF( B2.EQ.ZERO )
00389      $         GO TO 80
00390             DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
00391                IF( Z( I4 ).GT.Z( I4-2 ) )
00392      $            RETURN
00393                B1 = B1*( Z( I4 ) / Z( I4-2 ) )
00394                B2 = B2 + B1
00395                IF( HUNDRD*B1.LT.B2 )
00396      $            GO TO 80
00397    70       CONTINUE
00398    80       CONTINUE
00399             B2 = SQRT( CNST3*B2 )
00400             A2 = DMIN2 / ( ONE+B2**2 )
00401             GAP2 = Z( NN-7 ) + Z( NN-9 ) -
00402      $             SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
00403             IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
00404                S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
00405             ELSE 
00406                S = MAX( S, A2*( ONE-CNST2*B2 ) )
00407             END IF
00408          ELSE
00409             S = QURTR*DMIN2
00410             TTYPE = -11
00411          END IF
00412       ELSE IF( N0IN.GT.( N0+2 ) ) THEN
00413 *
00414 *        Case 12, more than two eigenvalues deflated. No information.
00415 *
00416          S = ZERO 
00417          TTYPE = -12
00418       END IF
00419 *
00420       TAU = S
00421       RETURN
00422 *
00423 *     End of DLASQ4
00424 *
00425       END
 All Files Functions