Return to zlaic1.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    1:       SUBROUTINE ZLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
    2: *
    3: *  -- LAPACK auxiliary routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            J, JOB
   10:       DOUBLE PRECISION   SEST, SESTPR
   11:       COMPLEX*16         C, GAMMA, S
   12: *     ..
   13: *     .. Array Arguments ..
   14:       COMPLEX*16         W( J ), X( J )
   15: *     ..
   16: *
   17: *  Purpose
   18: *  =======
   19: *
   20: *  ZLAIC1 applies one step of incremental condition estimation in
   21: *  its simplest version:
   22: *
   23: *  Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
   24: *  lower triangular matrix L, such that
   25: *           twonorm(L*x) = sest
   26: *  Then ZLAIC1 computes sestpr, s, c such that
   27: *  the vector
   28: *                  [ s*x ]
   29: *           xhat = [  c  ]
   30: *  is an approximate singular vector of
   31: *                  [ L     0  ]
   32: *           Lhat = [ w' gamma ]
   33: *  in the sense that
   34: *           twonorm(Lhat*xhat) = sestpr.
   35: *
   36: *  Depending on JOB, an estimate for the largest or smallest singular
   37: *  value is computed.
   38: *
   39: *  Note that [s c]' and sestpr**2 is an eigenpair of the system
   40: *
   41: *      diag(sest*sest, 0) + [alpha  gamma] * [ conjg(alpha) ]
   42: *                                            [ conjg(gamma) ]
   43: *
   44: *  where  alpha =  conjg(x)'*w.
   45: *
   46: *  Arguments
   47: *  =========
   48: *
   49: *  JOB     (input) INTEGER
   50: *          = 1: an estimate for the largest singular value is computed.
   51: *          = 2: an estimate for the smallest singular value is computed.
   52: *
   53: *  J       (input) INTEGER
   54: *          Length of X and W
   55: *
   56: *  X       (input) COMPLEX*16 array, dimension (J)
   57: *          The j-vector x.
   58: *
   59: *  SEST    (input) DOUBLE PRECISION
   60: *          Estimated singular value of j by j matrix L
   61: *
   62: *  W       (input) COMPLEX*16 array, dimension (J)
   63: *          The j-vector w.
   64: *
   65: *  GAMMA   (input) COMPLEX*16
   66: *          The diagonal element gamma.
   67: *
   68: *  SESTPR  (output) DOUBLE PRECISION
   69: *          Estimated singular value of (j+1) by (j+1) matrix Lhat.
   70: *
   71: *  S       (output) COMPLEX*16
   72: *          Sine needed in forming xhat.
   73: *
   74: *  C       (output) COMPLEX*16
   75: *          Cosine needed in forming xhat.
   76: *
   77: *  =====================================================================
   78: *
   79: *     .. Parameters ..
   80:       DOUBLE PRECISION   ZERO, ONE, TWO
   81:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
   82:       DOUBLE PRECISION   HALF, FOUR
   83:       PARAMETER          ( HALF = 0.5D0, FOUR = 4.0D0 )
   84: *     ..
   85: *     .. Local Scalars ..
   86:       DOUBLE PRECISION   ABSALP, ABSEST, ABSGAM, B, EPS, NORMA, S1, S2,
   87:      $                   SCL, T, TEST, TMP, ZETA1, ZETA2
   88:       COMPLEX*16         ALPHA, COSINE, SINE
   89: *     ..
   90: *     .. Intrinsic Functions ..
   91:       INTRINSIC          ABS, DCONJG, MAX, SQRT
   92: *     ..
   93: *     .. External Functions ..
   94:       DOUBLE PRECISION   DLAMCH
   95:       COMPLEX*16         ZDOTC
   96:       EXTERNAL           DLAMCH, ZDOTC
   97: *     ..
   98: *     .. Executable Statements ..
   99: *
  100:       EPS = DLAMCH( 'Epsilon' )
  101:       ALPHA = ZDOTC( J, X, 1, W, 1 )
  102: *
  103:       ABSALP = ABS( ALPHA )
  104:       ABSGAM = ABS( GAMMA )
  105:       ABSEST = ABS( SEST )
  106: *
  107:       IF( JOB.EQ.1 ) THEN
  108: *
  109: *        Estimating largest singular value
  110: *
  111: *        special cases
  112: *
  113:          IF( SEST.EQ.ZERO ) THEN
  114:             S1 = MAX( ABSGAM, ABSALP )
  115:             IF( S1.EQ.ZERO ) THEN
  116:                S = ZERO
  117:                C = ONE
  118:                SESTPR = ZERO
  119:             ELSE
  120:                S = ALPHA / S1
  121:                C = GAMMA / S1
  122:                TMP = SQRT( S*DCONJG( S )+C*DCONJG( C ) )
  123:                S = S / TMP
  124:                C = C / TMP
  125:                SESTPR = S1*TMP
  126:             END IF
  127:             RETURN
  128:          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
  129:             S = ONE
  130:             C = ZERO
  131:             TMP = MAX( ABSEST, ABSALP )
  132:             S1 = ABSEST / TMP
  133:             S2 = ABSALP / TMP
  134:             SESTPR = TMP*SQRT( S1*S1+S2*S2 )
  135:             RETURN
  136:          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
  137:             S1 = ABSGAM
  138:             S2 = ABSEST
  139:             IF( S1.LE.S2 ) THEN
  140:                S = ONE
  141:                C = ZERO
  142:                SESTPR = S2
  143:             ELSE
  144:                S = ZERO
  145:                C = ONE
  146:                SESTPR = S1
  147:             END IF
  148:             RETURN
  149:          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
  150:             S1 = ABSGAM
  151:             S2 = ABSALP
  152:             IF( S1.LE.S2 ) THEN
  153:                TMP = S1 / S2
  154:                SCL = SQRT( ONE+TMP*TMP )
  155:                SESTPR = S2*SCL
  156:                S = ( ALPHA / S2 ) / SCL
  157:                C = ( GAMMA / S2 ) / SCL
  158:             ELSE
  159:                TMP = S2 / S1
  160:                SCL = SQRT( ONE+TMP*TMP )
  161:                SESTPR = S1*SCL
  162:                S = ( ALPHA / S1 ) / SCL
  163:                C = ( GAMMA / S1 ) / SCL
  164:             END IF
  165:             RETURN
  166:          ELSE
  167: *
  168: *           normal case
  169: *
  170:             ZETA1 = ABSALP / ABSEST
  171:             ZETA2 = ABSGAM / ABSEST
  172: *
  173:             B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
  174:             C = ZETA1*ZETA1
  175:             IF( B.GT.ZERO ) THEN
  176:                T = C / ( B+SQRT( B*B+C ) )
  177:             ELSE
  178:                T = SQRT( B*B+C ) - B
  179:             END IF
  180: *
  181:             SINE = -( ALPHA / ABSEST ) / T
  182:             COSINE = -( GAMMA / ABSEST ) / ( ONE+T )
  183:             TMP = SQRT( SINE*DCONJG( SINE )+COSINE*DCONJG( COSINE ) )
  184:             S = SINE / TMP
  185:             C = COSINE / TMP
  186:             SESTPR = SQRT( T+ONE )*ABSEST
  187:             RETURN
  188:          END IF
  189: *
  190:       ELSE IF( JOB.EQ.2 ) THEN
  191: *
  192: *        Estimating smallest singular value
  193: *
  194: *        special cases
  195: *
  196:          IF( SEST.EQ.ZERO ) THEN
  197:             SESTPR = ZERO
  198:             IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
  199:                SINE = ONE
  200:                COSINE = ZERO
  201:             ELSE
  202:                SINE = -DCONJG( GAMMA )
  203:                COSINE = DCONJG( ALPHA )
  204:             END IF
  205:             S1 = MAX( ABS( SINE ), ABS( COSINE ) )
  206:             S = SINE / S1
  207:             C = COSINE / S1
  208:             TMP = SQRT( S*DCONJG( S )+C*DCONJG( C ) )
  209:             S = S / TMP
  210:             C = C / TMP
  211:             RETURN
  212:          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
  213:             S = ZERO
  214:             C = ONE
  215:             SESTPR = ABSGAM
  216:             RETURN
  217:          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
  218:             S1 = ABSGAM
  219:             S2 = ABSEST
  220:             IF( S1.LE.S2 ) THEN
  221:                S = ZERO
  222:                C = ONE
  223:                SESTPR = S1
  224:             ELSE
  225:                S = ONE
  226:                C = ZERO
  227:                SESTPR = S2
  228:             END IF
  229:             RETURN
  230:          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
  231:             S1 = ABSGAM
  232:             S2 = ABSALP
  233:             IF( S1.LE.S2 ) THEN
  234:                TMP = S1 / S2
  235:                SCL = SQRT( ONE+TMP*TMP )
  236:                SESTPR = ABSEST*( TMP / SCL )
  237:                S = -( DCONJG( GAMMA ) / S2 ) / SCL
  238:                C = ( DCONJG( ALPHA ) / S2 ) / SCL
  239:             ELSE
  240:                TMP = S2 / S1
  241:                SCL = SQRT( ONE+TMP*TMP )
  242:                SESTPR = ABSEST / SCL
  243:                S = -( DCONJG( GAMMA ) / S1 ) / SCL
  244:                C = ( DCONJG( ALPHA ) / S1 ) / SCL
  245:             END IF
  246:             RETURN
  247:          ELSE
  248: *
  249: *           normal case
  250: *
  251:             ZETA1 = ABSALP / ABSEST
  252:             ZETA2 = ABSGAM / ABSEST
  253: *
  254:             NORMA = MAX( ONE+ZETA1*ZETA1+ZETA1*ZETA2,
  255:      $              ZETA1*ZETA2+ZETA2*ZETA2 )
  256: *
  257: *           See if root is closer to zero or to ONE
  258: *
  259:             TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
  260:             IF( TEST.GE.ZERO ) THEN
  261: *
  262: *              root is close to zero, compute directly
  263: *
  264:                B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
  265:                C = ZETA2*ZETA2
  266:                T = C / ( B+SQRT( ABS( B*B-C ) ) )
  267:                SINE = ( ALPHA / ABSEST ) / ( ONE-T )
  268:                COSINE = -( GAMMA / ABSEST ) / T
  269:                SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
  270:             ELSE
  271: *
  272: *              root is closer to ONE, shift by that amount
  273: *
  274:                B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
  275:                C = ZETA1*ZETA1
  276:                IF( B.GE.ZERO ) THEN
  277:                   T = -C / ( B+SQRT( B*B+C ) )
  278:                ELSE
  279:                   T = B - SQRT( B*B+C )
  280:                END IF
  281:                SINE = -( ALPHA / ABSEST ) / T
  282:                COSINE = -( GAMMA / ABSEST ) / ( ONE+T )
  283:                SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
  284:             END IF
  285:             TMP = SQRT( SINE*DCONJG( SINE )+COSINE*DCONJG( COSINE ) )
  286:             S = SINE / TMP
  287:             C = COSINE / TMP
  288:             RETURN
  289: *
  290:          END IF
  291:       END IF
  292:       RETURN
  293: *
  294: *     End of ZLAIC1
  295: *
  296:       END