LAPACK: SRC/dlaic1.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            J, JOB
00010       DOUBLE PRECISION   C, GAMMA, S, SEST, SESTPR
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   W( J ), X( J )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DLAIC1 applies one step of incremental condition estimation in
00020 *  its simplest version:
00021 *
00022 *  Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
00023 *  lower triangular matrix L, such that
00024 *           twonorm(L*x) = sest
00025 *  Then DLAIC1 computes sestpr, s, c such that
00026 *  the vector
00027 *                  [ s*x ]
00028 *           xhat = [  c  ]
00029 *  is an approximate singular vector of
00030 *                  [ L     0  ]
00031 *           Lhat = [ w' gamma ]
00032 *  in the sense that
00033 *           twonorm(Lhat*xhat) = sestpr.
00034 *
00035 *  Depending on JOB, an estimate for the largest or smallest singular
00036 *  value is computed.
00037 *
00038 *  Note that [s c]' and sestpr**2 is an eigenpair of the system
00039 *
00040 *      diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]
00041 *                                            [ gamma ]
00042 *
00043 *  where  alpha =  x'*w.
00044 *
00045 *  Arguments
00046 *  =========
00047 *
00048 *  JOB     (input) INTEGER
00049 *          = 1: an estimate for the largest singular value is computed.
00050 *          = 2: an estimate for the smallest singular value is computed.
00051 *
00052 *  J       (input) INTEGER
00053 *          Length of X and W
00054 *
00055 *  X       (input) DOUBLE PRECISION array, dimension (J)
00056 *          The j-vector x.
00057 *
00058 *  SEST    (input) DOUBLE PRECISION
00059 *          Estimated singular value of j by j matrix L
00060 *
00061 *  W       (input) DOUBLE PRECISION array, dimension (J)
00062 *          The j-vector w.
00063 *
00064 *  GAMMA   (input) DOUBLE PRECISION
00065 *          The diagonal element gamma.
00066 *
00067 *  SESTPR  (output) DOUBLE PRECISION
00068 *          Estimated singular value of (j+1) by (j+1) matrix Lhat.
00069 *
00070 *  S       (output) DOUBLE PRECISION
00071 *          Sine needed in forming xhat.
00072 *
00073 *  C       (output) DOUBLE PRECISION
00074 *          Cosine needed in forming xhat.
00075 *
00076 *  =====================================================================
00077 *
00078 *     .. Parameters ..
00079       DOUBLE PRECISION   ZERO, ONE, TWO
00080       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
00081       DOUBLE PRECISION   HALF, FOUR
00082       PARAMETER          ( HALF = 0.5D0, FOUR = 4.0D0 )
00083 *     ..
00084 *     .. Local Scalars ..
00085       DOUBLE PRECISION   ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS,
00086      $                   NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
00087 *     ..
00088 *     .. Intrinsic Functions ..
00089       INTRINSIC          ABS, MAX, SIGN, SQRT
00090 *     ..
00091 *     .. External Functions ..
00092       DOUBLE PRECISION   DDOT, DLAMCH
00093       EXTERNAL           DDOT, DLAMCH
00094 *     ..
00095 *     .. Executable Statements ..
00096 *
00097       EPS = DLAMCH( 'Epsilon' )
00098       ALPHA = DDOT( J, X, 1, W, 1 )
00099 *
00100       ABSALP = ABS( ALPHA )
00101       ABSGAM = ABS( GAMMA )
00102       ABSEST = ABS( SEST )
00103 *
00104       IF( JOB.EQ.1 ) THEN
00105 *
00106 *        Estimating largest singular value
00107 *
00108 *        special cases
00109 *
00110          IF( SEST.EQ.ZERO ) THEN
00111             S1 = MAX( ABSGAM, ABSALP )
00112             IF( S1.EQ.ZERO ) THEN
00113                S = ZERO
00114                C = ONE
00115                SESTPR = ZERO
00116             ELSE
00117                S = ALPHA / S1
00118                C = GAMMA / S1
00119                TMP = SQRT( S*S+C*C )
00120                S = S / TMP
00121                C = C / TMP
00122                SESTPR = S1*TMP
00123             END IF
00124             RETURN
00125          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
00126             S = ONE
00127             C = ZERO
00128             TMP = MAX( ABSEST, ABSALP )
00129             S1 = ABSEST / TMP
00130             S2 = ABSALP / TMP
00131             SESTPR = TMP*SQRT( S1*S1+S2*S2 )
00132             RETURN
00133          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
00134             S1 = ABSGAM
00135             S2 = ABSEST
00136             IF( S1.LE.S2 ) THEN
00137                S = ONE
00138                C = ZERO
00139                SESTPR = S2
00140             ELSE
00141                S = ZERO
00142                C = ONE
00143                SESTPR = S1
00144             END IF
00145             RETURN
00146          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
00147             S1 = ABSGAM
00148             S2 = ABSALP
00149             IF( S1.LE.S2 ) THEN
00150                TMP = S1 / S2
00151                S = SQRT( ONE+TMP*TMP )
00152                SESTPR = S2*S
00153                C = ( GAMMA / S2 ) / S
00154                S = SIGN( ONE, ALPHA ) / S
00155             ELSE
00156                TMP = S2 / S1
00157                C = SQRT( ONE+TMP*TMP )
00158                SESTPR = S1*C
00159                S = ( ALPHA / S1 ) / C
00160                C = SIGN( ONE, GAMMA ) / C
00161             END IF
00162             RETURN
00163          ELSE
00164 *
00165 *           normal case
00166 *
00167             ZETA1 = ALPHA / ABSEST
00168             ZETA2 = GAMMA / ABSEST
00169 *
00170             B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
00171             C = ZETA1*ZETA1
00172             IF( B.GT.ZERO ) THEN
00173                T = C / ( B+SQRT( B*B+C ) )
00174             ELSE
00175                T = SQRT( B*B+C ) - B
00176             END IF
00177 *
00178             SINE = -ZETA1 / T
00179             COSINE = -ZETA2 / ( ONE+T )
00180             TMP = SQRT( SINE*SINE+COSINE*COSINE )
00181             S = SINE / TMP
00182             C = COSINE / TMP
00183             SESTPR = SQRT( T+ONE )*ABSEST
00184             RETURN
00185          END IF
00186 *
00187       ELSE IF( JOB.EQ.2 ) THEN
00188 *
00189 *        Estimating smallest singular value
00190 *
00191 *        special cases
00192 *
00193          IF( SEST.EQ.ZERO ) THEN
00194             SESTPR = ZERO
00195             IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
00196                SINE = ONE
00197                COSINE = ZERO
00198             ELSE
00199                SINE = -GAMMA
00200                COSINE = ALPHA
00201             END IF
00202             S1 = MAX( ABS( SINE ), ABS( COSINE ) )
00203             S = SINE / S1
00204             C = COSINE / S1
00205             TMP = SQRT( S*S+C*C )
00206             S = S / TMP
00207             C = C / TMP
00208             RETURN
00209          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
00210             S = ZERO
00211             C = ONE
00212             SESTPR = ABSGAM
00213             RETURN
00214          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
00215             S1 = ABSGAM
00216             S2 = ABSEST
00217             IF( S1.LE.S2 ) THEN
00218                S = ZERO
00219                C = ONE
00220                SESTPR = S1
00221             ELSE
00222                S = ONE
00223                C = ZERO
00224                SESTPR = S2
00225             END IF
00226             RETURN
00227          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
00228             S1 = ABSGAM
00229             S2 = ABSALP
00230             IF( S1.LE.S2 ) THEN
00231                TMP = S1 / S2
00232                C = SQRT( ONE+TMP*TMP )
00233                SESTPR = ABSEST*( TMP / C )
00234                S = -( GAMMA / S2 ) / C
00235                C = SIGN( ONE, ALPHA ) / C
00236             ELSE
00237                TMP = S2 / S1
00238                S = SQRT( ONE+TMP*TMP )
00239                SESTPR = ABSEST / S
00240                C = ( ALPHA / S1 ) / S
00241                S = -SIGN( ONE, GAMMA ) / S
00242             END IF
00243             RETURN
00244          ELSE
00245 *
00246 *           normal case
00247 *
00248             ZETA1 = ALPHA / ABSEST
00249             ZETA2 = GAMMA / ABSEST
00250 *
00251             NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ),
00252      $              ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 )
00253 *
00254 *           See if root is closer to zero or to ONE
00255 *
00256             TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
00257             IF( TEST.GE.ZERO ) THEN
00258 *
00259 *              root is close to zero, compute directly
00260 *
00261                B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
00262                C = ZETA2*ZETA2
00263                T = C / ( B+SQRT( ABS( B*B-C ) ) )
00264                SINE = ZETA1 / ( ONE-T )
00265                COSINE = -ZETA2 / T
00266                SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
00267             ELSE
00268 *
00269 *              root is closer to ONE, shift by that amount
00270 *
00271                B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
00272                C = ZETA1*ZETA1
00273                IF( B.GE.ZERO ) THEN
00274                   T = -C / ( B+SQRT( B*B+C ) )
00275                ELSE
00276                   T = B - SQRT( B*B+C )
00277                END IF
00278                SINE = -ZETA1 / T
00279                COSINE = -ZETA2 / ( ONE+T )
00280                SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
00281             END IF
00282             TMP = SQRT( SINE*SINE+COSINE*COSINE )
00283             S = SINE / TMP
00284             C = COSINE / TMP
00285             RETURN
00286 *
00287          END IF
00288       END IF
00289       RETURN
00290 *
00291 *     End of DLAIC1
00292 *
00293       END