LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlasq5.f
Go to the documentation of this file.
00001 *> \brief \b DLASQ5
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLASQ5 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq5.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq5.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq5.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
00022 *                          DNM1, DNM2, IEEE, EPS )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       LOGICAL            IEEE
00026 *       INTEGER            I0, N0, PP
00027 *       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       DOUBLE PRECISION   Z( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DLASQ5 computes one dqds transform in ping-pong form, one
00040 *> version for IEEE machines another for non IEEE machines.
00041 *> \endverbatim
00042 *
00043 *  Arguments:
00044 *  ==========
00045 *
00046 *> \param[in] I0
00047 *> \verbatim
00048 *>          I0 is INTEGER
00049 *>        First index.
00050 *> \endverbatim
00051 *>
00052 *> \param[in] N0
00053 *> \verbatim
00054 *>          N0 is INTEGER
00055 *>        Last index.
00056 *> \endverbatim
00057 *>
00058 *> \param[in] Z
00059 *> \verbatim
00060 *>          Z is DOUBLE PRECISION array, dimension ( 4*N )
00061 *>        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
00062 *>        an extra argument.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] PP
00066 *> \verbatim
00067 *>          PP is INTEGER
00068 *>        PP=0 for ping, PP=1 for pong.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] TAU
00072 *> \verbatim
00073 *>          TAU is DOUBLE PRECISION
00074 *>        This is the shift.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] SIGMA
00078 *> \verbatim
00079 *>          SIGMA is DOUBLE PRECISION
00080 *>        This is the accumulated shift up to this step.
00081 *> \endverbatim
00082 *>
00083 *> \param[out] DMIN
00084 *> \verbatim
00085 *>          DMIN is DOUBLE PRECISION
00086 *>        Minimum value of d.
00087 *> \endverbatim
00088 *>
00089 *> \param[out] DMIN1
00090 *> \verbatim
00091 *>          DMIN1 is DOUBLE PRECISION
00092 *>        Minimum value of d, excluding D( N0 ).
00093 *> \endverbatim
00094 *>
00095 *> \param[out] DMIN2
00096 *> \verbatim
00097 *>          DMIN2 is DOUBLE PRECISION
00098 *>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
00099 *> \endverbatim
00100 *>
00101 *> \param[out] DN
00102 *> \verbatim
00103 *>          DN is DOUBLE PRECISION
00104 *>        d(N0), the last value of d.
00105 *> \endverbatim
00106 *>
00107 *> \param[out] DNM1
00108 *> \verbatim
00109 *>          DNM1 is DOUBLE PRECISION
00110 *>        d(N0-1).
00111 *> \endverbatim
00112 *>
00113 *> \param[out] DNM2
00114 *> \verbatim
00115 *>          DNM2 is DOUBLE PRECISION
00116 *>        d(N0-2).
00117 *> \endverbatim
00118 *>
00119 *> \param[in] IEEE
00120 *> \verbatim
00121 *>          IEEE is LOGICAL
00122 *>        Flag for IEEE or non IEEE arithmetic.
00123 *> \endverbatim
00124 *
00125 *> \param[in] EPS
00126 *> \verbatim
00127 *>          EPS is DOUBLE PRECISION
00128 *>        This is the value of epsilon used.
00129 *> \endverbatim
00130 *>
00131 *  Authors:
00132 *  ========
00133 *
00134 *> \author Univ. of Tennessee 
00135 *> \author Univ. of California Berkeley 
00136 *> \author Univ. of Colorado Denver 
00137 *> \author NAG Ltd. 
00138 *
00139 *> \date April 2012
00140 *
00141 *> \ingroup auxOTHERcomputational
00142 *
00143 *  =====================================================================
00144       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
00145      $                   DN, DNM1, DNM2, IEEE, EPS )
00146 *
00147 *  -- LAPACK computational routine (version 3.4.1) --
00148 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00149 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00150 *     April 2012
00151 *
00152 *     .. Scalar Arguments ..
00153       LOGICAL            IEEE
00154       INTEGER            I0, N0, PP
00155       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
00156      $                   SIGMA, EPS
00157 *     ..
00158 *     .. Array Arguments ..
00159       DOUBLE PRECISION   Z( * )
00160 *     ..
00161 *
00162 *  =====================================================================
00163 *
00164 *     .. Parameter ..
00165       DOUBLE PRECISION   ZERO, HALF
00166       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5 )
00167 *     ..
00168 *     .. Local Scalars ..
00169       INTEGER            J4, J4P2
00170       DOUBLE PRECISION   D, EMIN, TEMP, DTHRESH
00171 *     ..
00172 *     .. Intrinsic Functions ..
00173       INTRINSIC          MIN
00174 *     ..
00175 *     .. Executable Statements ..
00176 *
00177       IF( ( N0-I0-1 ).LE.0 )
00178      $   RETURN
00179 *
00180       DTHRESH = EPS*(SIGMA+TAU)
00181       IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
00182       IF( TAU.NE.ZERO ) THEN
00183       J4 = 4*I0 + PP - 3
00184       EMIN = Z( J4+4 ) 
00185       D = Z( J4 ) - TAU
00186       DMIN = D
00187       DMIN1 = -Z( J4 )
00188 *
00189       IF( IEEE ) THEN
00190 *
00191 *        Code for IEEE arithmetic.
00192 *
00193          IF( PP.EQ.0 ) THEN
00194             DO 10 J4 = 4*I0, 4*( N0-3 ), 4
00195                Z( J4-2 ) = D + Z( J4-1 ) 
00196                TEMP = Z( J4+1 ) / Z( J4-2 )
00197                D = D*TEMP - TAU
00198                DMIN = MIN( DMIN, D )
00199                Z( J4 ) = Z( J4-1 )*TEMP
00200                EMIN = MIN( Z( J4 ), EMIN )
00201    10       CONTINUE
00202          ELSE
00203             DO 20 J4 = 4*I0, 4*( N0-3 ), 4
00204                Z( J4-3 ) = D + Z( J4 ) 
00205                TEMP = Z( J4+2 ) / Z( J4-3 )
00206                D = D*TEMP - TAU
00207                DMIN = MIN( DMIN, D )
00208                Z( J4-1 ) = Z( J4 )*TEMP
00209                EMIN = MIN( Z( J4-1 ), EMIN )
00210    20       CONTINUE
00211          END IF
00212 *
00213 *        Unroll last two steps. 
00214 *
00215          DNM2 = D
00216          DMIN2 = DMIN
00217          J4 = 4*( N0-2 ) - PP
00218          J4P2 = J4 + 2*PP - 1
00219          Z( J4-2 ) = DNM2 + Z( J4P2 )
00220          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00221          DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
00222          DMIN = MIN( DMIN, DNM1 )
00223 *
00224          DMIN1 = DMIN
00225          J4 = J4 + 4
00226          J4P2 = J4 + 2*PP - 1
00227          Z( J4-2 ) = DNM1 + Z( J4P2 )
00228          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00229          DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
00230          DMIN = MIN( DMIN, DN )
00231 *
00232       ELSE
00233 *
00234 *        Code for non IEEE arithmetic.
00235 *
00236          IF( PP.EQ.0 ) THEN
00237             DO 30 J4 = 4*I0, 4*( N0-3 ), 4
00238                Z( J4-2 ) = D + Z( J4-1 ) 
00239                IF( D.LT.ZERO ) THEN
00240                   RETURN
00241                ELSE 
00242                   Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
00243                   D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
00244                END IF
00245                DMIN = MIN( DMIN, D )
00246                EMIN = MIN( EMIN, Z( J4 ) )
00247    30       CONTINUE
00248          ELSE
00249             DO 40 J4 = 4*I0, 4*( N0-3 ), 4
00250                Z( J4-3 ) = D + Z( J4 ) 
00251                IF( D.LT.ZERO ) THEN
00252                   RETURN
00253                ELSE 
00254                   Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
00255                   D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
00256                END IF
00257                DMIN = MIN( DMIN, D )
00258                EMIN = MIN( EMIN, Z( J4-1 ) )
00259    40       CONTINUE
00260          END IF
00261 *
00262 *        Unroll last two steps. 
00263 *
00264          DNM2 = D
00265          DMIN2 = DMIN
00266          J4 = 4*( N0-2 ) - PP
00267          J4P2 = J4 + 2*PP - 1
00268          Z( J4-2 ) = DNM2 + Z( J4P2 )
00269          IF( DNM2.LT.ZERO ) THEN
00270             RETURN
00271          ELSE
00272             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00273             DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
00274          END IF
00275          DMIN = MIN( DMIN, DNM1 )
00276 *
00277          DMIN1 = DMIN
00278          J4 = J4 + 4
00279          J4P2 = J4 + 2*PP - 1
00280          Z( J4-2 ) = DNM1 + Z( J4P2 )
00281          IF( DNM1.LT.ZERO ) THEN
00282             RETURN
00283          ELSE
00284             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00285             DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
00286          END IF
00287          DMIN = MIN( DMIN, DN )
00288 *
00289       END IF
00290       ELSE
00291 *     This is the version that sets d's to zero if they are small enough
00292          J4 = 4*I0 + PP - 3
00293          EMIN = Z( J4+4 ) 
00294          D = Z( J4 ) - TAU
00295          DMIN = D
00296          DMIN1 = -Z( J4 )
00297          IF( IEEE ) THEN
00298 *     
00299 *     Code for IEEE arithmetic.
00300 *     
00301             IF( PP.EQ.0 ) THEN
00302                DO 50 J4 = 4*I0, 4*( N0-3 ), 4
00303                   Z( J4-2 ) = D + Z( J4-1 ) 
00304                   TEMP = Z( J4+1 ) / Z( J4-2 )
00305                   D = D*TEMP - TAU
00306                   IF( D.LT.DTHRESH ) D = ZERO
00307                   DMIN = MIN( DMIN, D )
00308                   Z( J4 ) = Z( J4-1 )*TEMP
00309                   EMIN = MIN( Z( J4 ), EMIN )
00310  50            CONTINUE
00311             ELSE
00312                DO 60 J4 = 4*I0, 4*( N0-3 ), 4
00313                   Z( J4-3 ) = D + Z( J4 ) 
00314                   TEMP = Z( J4+2 ) / Z( J4-3 )
00315                   D = D*TEMP - TAU
00316                   IF( D.LT.DTHRESH ) D = ZERO
00317                   DMIN = MIN( DMIN, D )
00318                   Z( J4-1 ) = Z( J4 )*TEMP
00319                   EMIN = MIN( Z( J4-1 ), EMIN )
00320  60            CONTINUE
00321             END IF
00322 *     
00323 *     Unroll last two steps. 
00324 *     
00325             DNM2 = D
00326             DMIN2 = DMIN
00327             J4 = 4*( N0-2 ) - PP
00328             J4P2 = J4 + 2*PP - 1
00329             Z( J4-2 ) = DNM2 + Z( J4P2 )
00330             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00331             DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
00332             DMIN = MIN( DMIN, DNM1 )
00333 *     
00334             DMIN1 = DMIN
00335             J4 = J4 + 4
00336             J4P2 = J4 + 2*PP - 1
00337             Z( J4-2 ) = DNM1 + Z( J4P2 )
00338             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00339             DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
00340             DMIN = MIN( DMIN, DN )
00341 *     
00342          ELSE
00343 *     
00344 *     Code for non IEEE arithmetic.
00345 *     
00346             IF( PP.EQ.0 ) THEN
00347                DO 70 J4 = 4*I0, 4*( N0-3 ), 4
00348                   Z( J4-2 ) = D + Z( J4-1 ) 
00349                   IF( D.LT.ZERO ) THEN
00350                      RETURN
00351                   ELSE 
00352                      Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
00353                      D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
00354                   END IF
00355                   IF( D.LT.DTHRESH) D = ZERO
00356                   DMIN = MIN( DMIN, D )
00357                   EMIN = MIN( EMIN, Z( J4 ) )
00358  70            CONTINUE
00359             ELSE
00360                DO 80 J4 = 4*I0, 4*( N0-3 ), 4
00361                   Z( J4-3 ) = D + Z( J4 ) 
00362                   IF( D.LT.ZERO ) THEN
00363                      RETURN
00364                   ELSE 
00365                      Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
00366                      D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
00367                   END IF
00368                   IF( D.LT.DTHRESH) D = ZERO
00369                   DMIN = MIN( DMIN, D )
00370                   EMIN = MIN( EMIN, Z( J4-1 ) )
00371  80            CONTINUE
00372             END IF
00373 *     
00374 *     Unroll last two steps. 
00375 *     
00376             DNM2 = D
00377             DMIN2 = DMIN
00378             J4 = 4*( N0-2 ) - PP
00379             J4P2 = J4 + 2*PP - 1
00380             Z( J4-2 ) = DNM2 + Z( J4P2 )
00381             IF( DNM2.LT.ZERO ) THEN
00382                RETURN
00383             ELSE
00384                Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00385                DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
00386             END IF
00387             DMIN = MIN( DMIN, DNM1 )
00388 *     
00389             DMIN1 = DMIN
00390             J4 = J4 + 4
00391             J4P2 = J4 + 2*PP - 1
00392             Z( J4-2 ) = DNM1 + Z( J4P2 )
00393             IF( DNM1.LT.ZERO ) THEN
00394                RETURN
00395             ELSE
00396                Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00397                DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
00398             END IF
00399             DMIN = MIN( DMIN, DN )
00400 *     
00401          END IF
00402       END IF
00403 *     
00404       Z( J4+2 ) = DN
00405       Z( 4*N0-PP ) = EMIN
00406       RETURN
00407 *
00408 *     End of DLASQ5
00409 *
00410       END
 All Files Functions