SRC/dlaic1.f

   1 *> \brief \b DLAIC1 applies one step of incremental condition estimation.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAIC1 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaic1.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaic1.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaic1.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            J, JOB
  25 *       DOUBLE PRECISION   C, GAMMA, S, SEST, SESTPR
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   W( J ), X( J )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DLAIC1 applies one step of incremental condition estimation in
  38 *> its simplest version:
  39 *>
  40 *> Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
  41 *> lower triangular matrix L, such that
  42 *>          twonorm(L*x) = sest
  43 *> Then DLAIC1 computes sestpr, s, c such that
  44 *> the vector
  45 *>                 [ s*x ]
  46 *>          xhat = [  c  ]
  47 *> is an approximate singular vector of
  48 *>                 [ L       0  ]
  49 *>          Lhat = [ w**T gamma ]
  50 *> in the sense that
  51 *>          twonorm(Lhat*xhat) = sestpr.
  52 *>
  53 *> Depending on JOB, an estimate for the largest or smallest singular
  54 *> value is computed.
  55 *>
  56 *> Note that [s c]**T and sestpr**2 is an eigenpair of the system
  57 *>
  58 *>     diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]
  59 *>                                           [ gamma ]
  60 *>
  61 *> where  alpha =  x**T*w.
  62 *> \endverbatim
  63 *
  64 *  Arguments:
  65 *  ==========
  66 *
  67 *> \param[in] JOB
  68 *> \verbatim
  69 *>          JOB is INTEGER
  70 *>          = 1: an estimate for the largest singular value is computed.
  71 *>          = 2: an estimate for the smallest singular value is computed.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] J
  75 *> \verbatim
  76 *>          J is INTEGER
  77 *>          Length of X and W
  78 *> \endverbatim
  79 *>
  80 *> \param[in] X
  81 *> \verbatim
  82 *>          X is DOUBLE PRECISION array, dimension (J)
  83 *>          The j-vector x.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] SEST
  87 *> \verbatim
  88 *>          SEST is DOUBLE PRECISION
  89 *>          Estimated singular value of j by j matrix L
  90 *> \endverbatim
  91 *>
  92 *> \param[in] W
  93 *> \verbatim
  94 *>          W is DOUBLE PRECISION array, dimension (J)
  95 *>          The j-vector w.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] GAMMA
  99 *> \verbatim
 100 *>          GAMMA is DOUBLE PRECISION
 101 *>          The diagonal element gamma.
 102 *> \endverbatim
 103 *>
 104 *> \param[out] SESTPR
 105 *> \verbatim
 106 *>          SESTPR is DOUBLE PRECISION
 107 *>          Estimated singular value of (j+1) by (j+1) matrix Lhat.
 108 *> \endverbatim
 109 *>
 110 *> \param[out] S
 111 *> \verbatim
 112 *>          S is DOUBLE PRECISION
 113 *>          Sine needed in forming xhat.
 114 *> \endverbatim
 115 *>
 116 *> \param[out] C
 117 *> \verbatim
 118 *>          C is DOUBLE PRECISION
 119 *>          Cosine needed in forming xhat.
 120 *> \endverbatim
 121 *
 122 *  Authors:
 123 *  ========
 124 *
 125 *> \author Univ. of Tennessee
 126 *> \author Univ. of California Berkeley
 127 *> \author Univ. of Colorado Denver
 128 *> \author NAG Ltd.
 129 *
 130 *> \date September 2012
 131 *
 132 *> \ingroup doubleOTHERauxiliary
 133 *
 134 *  =====================================================================
 135       SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
 136 *
 137 *  -- LAPACK auxiliary routine (version 3.4.2) --
 138 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 139 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 140 *     September 2012
 141 *
 142 *     .. Scalar Arguments ..
 143       INTEGER            J, JOB
 144       DOUBLE PRECISION   C, GAMMA, S, SEST, SESTPR
 145 *     ..
 146 *     .. Array Arguments ..
 147       DOUBLE PRECISION   W( J ), X( J )
 148 *     ..
 149 *
 150 *  =====================================================================
 151 *
 152 *     .. Parameters ..
 153       DOUBLE PRECISION   ZERO, ONE, TWO
 154       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
 155       DOUBLE PRECISION   HALF, FOUR
 156       PARAMETER          ( HALF = 0.5D0, FOUR = 4.0D0 )
 157 *     ..
 158 *     .. Local Scalars ..
 159       DOUBLE PRECISION   ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS,
 160      $                   NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
 161 *     ..
 162 *     .. Intrinsic Functions ..
 163       INTRINSIC          ABS, MAX, SIGN, SQRT
 164 *     ..
 165 *     .. External Functions ..
 166       DOUBLE PRECISION   DDOT, DLAMCH
 167       EXTERNAL           DDOT, DLAMCH
 168 *     ..
 169 *     .. Executable Statements ..
 170 *
 171       EPS = DLAMCH( 'Epsilon' )
 172       ALPHA = DDOT( J, X, 1, W, 1 )
 173 *
 174       ABSALP = ABS( ALPHA )
 175       ABSGAM = ABS( GAMMA )
 176       ABSEST = ABS( SEST )
 177 *
 178       IF( JOB.EQ.1 ) THEN
 179 *
 180 *        Estimating largest singular value
 181 *
 182 *        special cases
 183 *
 184          IF( SEST.EQ.ZERO ) THEN
 185             S1 = MAX( ABSGAM, ABSALP )
 186             IF( S1.EQ.ZERO ) THEN
 187                S = ZERO
 188                C = ONE
 189                SESTPR = ZERO
 190             ELSE
 191                S = ALPHA / S1
 192                C = GAMMA / S1
 193                TMP = SQRT( S*S+C*C )
 194                S = S / TMP
 195                C = C / TMP
 196                SESTPR = S1*TMP
 197             END IF
 198             RETURN
 199          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
 200             S = ONE
 201             C = ZERO
 202             TMP = MAX( ABSEST, ABSALP )
 203             S1 = ABSEST / TMP
 204             S2 = ABSALP / TMP
 205             SESTPR = TMP*SQRT( S1*S1+S2*S2 )
 206             RETURN
 207          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
 208             S1 = ABSGAM
 209             S2 = ABSEST
 210             IF( S1.LE.S2 ) THEN
 211                S = ONE
 212                C = ZERO
 213                SESTPR = S2
 214             ELSE
 215                S = ZERO
 216                C = ONE
 217                SESTPR = S1
 218             END IF
 219             RETURN
 220          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
 221             S1 = ABSGAM
 222             S2 = ABSALP
 223             IF( S1.LE.S2 ) THEN
 224                TMP = S1 / S2
 225                S = SQRT( ONE+TMP*TMP )
 226                SESTPR = S2*S
 227                C = ( GAMMA / S2 ) / S
 228                S = SIGN( ONE, ALPHA ) / S
 229             ELSE
 230                TMP = S2 / S1
 231                C = SQRT( ONE+TMP*TMP )
 232                SESTPR = S1*C
 233                S = ( ALPHA / S1 ) / C
 234                C = SIGN( ONE, GAMMA ) / C
 235             END IF
 236             RETURN
 237          ELSE
 238 *
 239 *           normal case
 240 *
 241             ZETA1 = ALPHA / ABSEST
 242             ZETA2 = GAMMA / ABSEST
 243 *
 244             B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
 245             C = ZETA1*ZETA1
 246             IF( B.GT.ZERO ) THEN
 247                T = C / ( B+SQRT( B*B+C ) )
 248             ELSE
 249                T = SQRT( B*B+C ) - B
 250             END IF
 251 *
 252             SINE = -ZETA1 / T
 253             COSINE = -ZETA2 / ( ONE+T )
 254             TMP = SQRT( SINE*SINE+COSINE*COSINE )
 255             S = SINE / TMP
 256             C = COSINE / TMP
 257             SESTPR = SQRT( T+ONE )*ABSEST
 258             RETURN
 259          END IF
 260 *
 261       ELSE IF( JOB.EQ.2 ) THEN
 262 *
 263 *        Estimating smallest singular value
 264 *
 265 *        special cases
 266 *
 267          IF( SEST.EQ.ZERO ) THEN
 268             SESTPR = ZERO
 269             IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
 270                SINE = ONE
 271                COSINE = ZERO
 272             ELSE
 273                SINE = -GAMMA
 274                COSINE = ALPHA
 275             END IF
 276             S1 = MAX( ABS( SINE ), ABS( COSINE ) )
 277             S = SINE / S1
 278             C = COSINE / S1
 279             TMP = SQRT( S*S+C*C )
 280             S = S / TMP
 281             C = C / TMP
 282             RETURN
 283          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
 284             S = ZERO
 285             C = ONE
 286             SESTPR = ABSGAM
 287             RETURN
 288          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
 289             S1 = ABSGAM
 290             S2 = ABSEST
 291             IF( S1.LE.S2 ) THEN
 292                S = ZERO
 293                C = ONE
 294                SESTPR = S1
 295             ELSE
 296                S = ONE
 297                C = ZERO
 298                SESTPR = S2
 299             END IF
 300             RETURN
 301          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
 302             S1 = ABSGAM
 303             S2 = ABSALP
 304             IF( S1.LE.S2 ) THEN
 305                TMP = S1 / S2
 306                C = SQRT( ONE+TMP*TMP )
 307                SESTPR = ABSEST*( TMP / C )
 308                S = -( GAMMA / S2 ) / C
 309                C = SIGN( ONE, ALPHA ) / C
 310             ELSE
 311                TMP = S2 / S1
 312                S = SQRT( ONE+TMP*TMP )
 313                SESTPR = ABSEST / S
 314                C = ( ALPHA / S1 ) / S
 315                S = -SIGN( ONE, GAMMA ) / S
 316             END IF
 317             RETURN
 318          ELSE
 319 *
 320 *           normal case
 321 *
 322             ZETA1 = ALPHA / ABSEST
 323             ZETA2 = GAMMA / ABSEST
 324 *
 325             NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ),
 326      $              ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 )
 327 *
 328 *           See if root is closer to zero or to ONE
 329 *
 330             TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
 331             IF( TEST.GE.ZERO ) THEN
 332 *
 333 *              root is close to zero, compute directly
 334 *
 335                B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
 336                C = ZETA2*ZETA2
 337                T = C / ( B+SQRT( ABS( B*B-C ) ) )
 338                SINE = ZETA1 / ( ONE-T )
 339                COSINE = -ZETA2 / T
 340                SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
 341             ELSE
 342 *
 343 *              root is closer to ONE, shift by that amount
 344 *
 345                B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
 346                C = ZETA1*ZETA1
 347                IF( B.GE.ZERO ) THEN
 348                   T = -C / ( B+SQRT( B*B+C ) )
 349                ELSE
 350                   T = B - SQRT( B*B+C )
 351                END IF
 352                SINE = -ZETA1 / T
 353                COSINE = -ZETA2 / ( ONE+T )
 354                SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
 355             END IF
 356             TMP = SQRT( SINE*SINE+COSINE*COSINE )
 357             S = SINE / TMP
 358             C = COSINE / TMP
 359             RETURN
 360 *
 361          END IF
 362       END IF
 363       RETURN
 364 *
 365 *     End of DLAIC1
 366 *
 367       END