SRC/zlaic1.f

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