LAPACK: SRC/claic1.f Source File

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