LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlaic1.f
Go to the documentation of this file.
00001 *> \brief \b DLAIC1
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLAIC1 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaic1.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaic1.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaic1.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            J, JOB
00025 *       DOUBLE PRECISION   C, GAMMA, S, SEST, SESTPR
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       DOUBLE PRECISION   W( J ), X( J )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> DLAIC1 applies one step of incremental condition estimation in
00038 *> its simplest version:
00039 *>
00040 *> Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
00041 *> lower triangular matrix L, such that
00042 *>          twonorm(L*x) = sest
00043 *> Then DLAIC1 computes sestpr, s, c such that
00044 *> the vector
00045 *>                 [ s*x ]
00046 *>          xhat = [  c  ]
00047 *> is an approximate singular vector of
00048 *>                 [ L       0  ]
00049 *>          Lhat = [ w**T gamma ]
00050 *> in the sense that
00051 *>          twonorm(Lhat*xhat) = sestpr.
00052 *>
00053 *> Depending on JOB, an estimate for the largest or smallest singular
00054 *> value is computed.
00055 *>
00056 *> Note that [s c]**T and sestpr**2 is an eigenpair of the system
00057 *>
00058 *>     diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]
00059 *>                                           [ gamma ]
00060 *>
00061 *> where  alpha =  x**T*w.
00062 *> \endverbatim
00063 *
00064 *  Arguments:
00065 *  ==========
00066 *
00067 *> \param[in] JOB
00068 *> \verbatim
00069 *>          JOB is INTEGER
00070 *>          = 1: an estimate for the largest singular value is computed.
00071 *>          = 2: an estimate for the smallest singular value is computed.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] J
00075 *> \verbatim
00076 *>          J is INTEGER
00077 *>          Length of X and W
00078 *> \endverbatim
00079 *>
00080 *> \param[in] X
00081 *> \verbatim
00082 *>          X is DOUBLE PRECISION array, dimension (J)
00083 *>          The j-vector x.
00084 *> \endverbatim
00085 *>
00086 *> \param[in] SEST
00087 *> \verbatim
00088 *>          SEST is DOUBLE PRECISION
00089 *>          Estimated singular value of j by j matrix L
00090 *> \endverbatim
00091 *>
00092 *> \param[in] W
00093 *> \verbatim
00094 *>          W is DOUBLE PRECISION array, dimension (J)
00095 *>          The j-vector w.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] GAMMA
00099 *> \verbatim
00100 *>          GAMMA is DOUBLE PRECISION
00101 *>          The diagonal element gamma.
00102 *> \endverbatim
00103 *>
00104 *> \param[out] SESTPR
00105 *> \verbatim
00106 *>          SESTPR is DOUBLE PRECISION
00107 *>          Estimated singular value of (j+1) by (j+1) matrix Lhat.
00108 *> \endverbatim
00109 *>
00110 *> \param[out] S
00111 *> \verbatim
00112 *>          S is DOUBLE PRECISION
00113 *>          Sine needed in forming xhat.
00114 *> \endverbatim
00115 *>
00116 *> \param[out] C
00117 *> \verbatim
00118 *>          C is DOUBLE PRECISION
00119 *>          Cosine needed in forming xhat.
00120 *> \endverbatim
00121 *
00122 *  Authors:
00123 *  ========
00124 *
00125 *> \author Univ. of Tennessee 
00126 *> \author Univ. of California Berkeley 
00127 *> \author Univ. of Colorado Denver 
00128 *> \author NAG Ltd. 
00129 *
00130 *> \date November 2011
00131 *
00132 *> \ingroup doubleOTHERauxiliary
00133 *
00134 *  =====================================================================
00135       SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
00136 *
00137 *  -- LAPACK auxiliary routine (version 3.4.0) --
00138 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00139 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00140 *     November 2011
00141 *
00142 *     .. Scalar Arguments ..
00143       INTEGER            J, JOB
00144       DOUBLE PRECISION   C, GAMMA, S, SEST, SESTPR
00145 *     ..
00146 *     .. Array Arguments ..
00147       DOUBLE PRECISION   W( J ), X( J )
00148 *     ..
00149 *
00150 *  =====================================================================
00151 *
00152 *     .. Parameters ..
00153       DOUBLE PRECISION   ZERO, ONE, TWO
00154       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
00155       DOUBLE PRECISION   HALF, FOUR
00156       PARAMETER          ( HALF = 0.5D0, FOUR = 4.0D0 )
00157 *     ..
00158 *     .. Local Scalars ..
00159       DOUBLE PRECISION   ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS,
00160      $                   NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
00161 *     ..
00162 *     .. Intrinsic Functions ..
00163       INTRINSIC          ABS, MAX, SIGN, SQRT
00164 *     ..
00165 *     .. External Functions ..
00166       DOUBLE PRECISION   DDOT, DLAMCH
00167       EXTERNAL           DDOT, DLAMCH
00168 *     ..
00169 *     .. Executable Statements ..
00170 *
00171       EPS = DLAMCH( 'Epsilon' )
00172       ALPHA = DDOT( J, X, 1, W, 1 )
00173 *
00174       ABSALP = ABS( ALPHA )
00175       ABSGAM = ABS( GAMMA )
00176       ABSEST = ABS( SEST )
00177 *
00178       IF( JOB.EQ.1 ) THEN
00179 *
00180 *        Estimating largest singular value
00181 *
00182 *        special cases
00183 *
00184          IF( SEST.EQ.ZERO ) THEN
00185             S1 = MAX( ABSGAM, ABSALP )
00186             IF( S1.EQ.ZERO ) THEN
00187                S = ZERO
00188                C = ONE
00189                SESTPR = ZERO
00190             ELSE
00191                S = ALPHA / S1
00192                C = GAMMA / S1
00193                TMP = SQRT( S*S+C*C )
00194                S = S / TMP
00195                C = C / TMP
00196                SESTPR = S1*TMP
00197             END IF
00198             RETURN
00199          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
00200             S = ONE
00201             C = ZERO
00202             TMP = MAX( ABSEST, ABSALP )
00203             S1 = ABSEST / TMP
00204             S2 = ABSALP / TMP
00205             SESTPR = TMP*SQRT( S1*S1+S2*S2 )
00206             RETURN
00207          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
00208             S1 = ABSGAM
00209             S2 = ABSEST
00210             IF( S1.LE.S2 ) THEN
00211                S = ONE
00212                C = ZERO
00213                SESTPR = S2
00214             ELSE
00215                S = ZERO
00216                C = ONE
00217                SESTPR = S1
00218             END IF
00219             RETURN
00220          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
00221             S1 = ABSGAM
00222             S2 = ABSALP
00223             IF( S1.LE.S2 ) THEN
00224                TMP = S1 / S2
00225                S = SQRT( ONE+TMP*TMP )
00226                SESTPR = S2*S
00227                C = ( GAMMA / S2 ) / S
00228                S = SIGN( ONE, ALPHA ) / S
00229             ELSE
00230                TMP = S2 / S1
00231                C = SQRT( ONE+TMP*TMP )
00232                SESTPR = S1*C
00233                S = ( ALPHA / S1 ) / C
00234                C = SIGN( ONE, GAMMA ) / C
00235             END IF
00236             RETURN
00237          ELSE
00238 *
00239 *           normal case
00240 *
00241             ZETA1 = ALPHA / ABSEST
00242             ZETA2 = GAMMA / ABSEST
00243 *
00244             B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
00245             C = ZETA1*ZETA1
00246             IF( B.GT.ZERO ) THEN
00247                T = C / ( B+SQRT( B*B+C ) )
00248             ELSE
00249                T = SQRT( B*B+C ) - B
00250             END IF
00251 *
00252             SINE = -ZETA1 / T
00253             COSINE = -ZETA2 / ( ONE+T )
00254             TMP = SQRT( SINE*SINE+COSINE*COSINE )
00255             S = SINE / TMP
00256             C = COSINE / TMP
00257             SESTPR = SQRT( T+ONE )*ABSEST
00258             RETURN
00259          END IF
00260 *
00261       ELSE IF( JOB.EQ.2 ) THEN
00262 *
00263 *        Estimating smallest singular value
00264 *
00265 *        special cases
00266 *
00267          IF( SEST.EQ.ZERO ) THEN
00268             SESTPR = ZERO
00269             IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
00270                SINE = ONE
00271                COSINE = ZERO
00272             ELSE
00273                SINE = -GAMMA
00274                COSINE = ALPHA
00275             END IF
00276             S1 = MAX( ABS( SINE ), ABS( COSINE ) )
00277             S = SINE / S1
00278             C = COSINE / S1
00279             TMP = SQRT( S*S+C*C )
00280             S = S / TMP
00281             C = C / TMP
00282             RETURN
00283          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
00284             S = ZERO
00285             C = ONE
00286             SESTPR = ABSGAM
00287             RETURN
00288          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
00289             S1 = ABSGAM
00290             S2 = ABSEST
00291             IF( S1.LE.S2 ) THEN
00292                S = ZERO
00293                C = ONE
00294                SESTPR = S1
00295             ELSE
00296                S = ONE
00297                C = ZERO
00298                SESTPR = S2
00299             END IF
00300             RETURN
00301          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
00302             S1 = ABSGAM
00303             S2 = ABSALP
00304             IF( S1.LE.S2 ) THEN
00305                TMP = S1 / S2
00306                C = SQRT( ONE+TMP*TMP )
00307                SESTPR = ABSEST*( TMP / C )
00308                S = -( GAMMA / S2 ) / C
00309                C = SIGN( ONE, ALPHA ) / C
00310             ELSE
00311                TMP = S2 / S1
00312                S = SQRT( ONE+TMP*TMP )
00313                SESTPR = ABSEST / S
00314                C = ( ALPHA / S1 ) / S
00315                S = -SIGN( ONE, GAMMA ) / S
00316             END IF
00317             RETURN
00318          ELSE
00319 *
00320 *           normal case
00321 *
00322             ZETA1 = ALPHA / ABSEST
00323             ZETA2 = GAMMA / ABSEST
00324 *
00325             NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ),
00326      $              ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 )
00327 *
00328 *           See if root is closer to zero or to ONE
00329 *
00330             TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
00331             IF( TEST.GE.ZERO ) THEN
00332 *
00333 *              root is close to zero, compute directly
00334 *
00335                B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
00336                C = ZETA2*ZETA2
00337                T = C / ( B+SQRT( ABS( B*B-C ) ) )
00338                SINE = ZETA1 / ( ONE-T )
00339                COSINE = -ZETA2 / T
00340                SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
00341             ELSE
00342 *
00343 *              root is closer to ONE, shift by that amount
00344 *
00345                B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
00346                C = ZETA1*ZETA1
00347                IF( B.GE.ZERO ) THEN
00348                   T = -C / ( B+SQRT( B*B+C ) )
00349                ELSE
00350                   T = B - SQRT( B*B+C )
00351                END IF
00352                SINE = -ZETA1 / T
00353                COSINE = -ZETA2 / ( ONE+T )
00354                SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
00355             END IF
00356             TMP = SQRT( SINE*SINE+COSINE*COSINE )
00357             S = SINE / TMP
00358             C = COSINE / TMP
00359             RETURN
00360 *
00361          END IF
00362       END IF
00363       RETURN
00364 *
00365 *     End of DLAIC1
00366 *
00367       END
 All Files Functions