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