Return to zla_porcond_x.f CVS log

Up to [local] / rpl / lapack / lapack

Patches pour OS/2

    1:       DOUBLE PRECISION FUNCTION ZLA_PORCOND_X( UPLO, N, A, LDA, AF,
    2:      $                                         LDAF, X, INFO, WORK,
    3:      $                                         RWORK )
    4: *
    5: *     -- LAPACK routine (version 3.2.1)                                 --
    6: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
    7: *     -- Jason Riedy of Univ. of California Berkeley.                 --
    8: *     -- April 2009                                                   --
    9: *
   10: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
   11: *     -- Univ. of California Berkeley and NAG Ltd.                    --
   12: *
   13:       IMPLICIT NONE
   14: *     ..
   15: *     .. Scalar Arguments ..
   16:       CHARACTER          UPLO
   17:       INTEGER            N, LDA, LDAF, INFO
   18: *     ..
   19: *     .. Array Arguments ..
   20:       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * )
   21:       DOUBLE PRECISION   RWORK( * )
   22: *     ..
   23: *
   24: *  Purpose
   25: *  =======
   26: *
   27: *     ZLA_PORCOND_X Computes the infinity norm condition number of
   28: *     op(A) * diag(X) where X is a COMPLEX*16 vector.
   29: *
   30: *  Arguments
   31: *  =========
   32: *
   33: *     UPLO    (input) CHARACTER*1
   34: *       = 'U':  Upper triangle of A is stored;
   35: *       = 'L':  Lower triangle of A is stored.
   36: *
   37: *     N       (input) INTEGER
   38: *     The number of linear equations, i.e., the order of the
   39: *     matrix A.  N >= 0.
   40: *
   41: *     A       (input) COMPLEX*16 array, dimension (LDA,N)
   42: *     On entry, the N-by-N matrix A.
   43: *
   44: *     LDA     (input) INTEGER
   45: *     The leading dimension of the array A.  LDA >= max(1,N).
   46: *
   47: *     AF      (input) COMPLEX*16 array, dimension (LDAF,N)
   48: *     The triangular factor U or L from the Cholesky factorization
   49: *     A = U**T*U or A = L*L**T, as computed by ZPOTRF.
   50: *
   51: *     LDAF    (input) INTEGER
   52: *     The leading dimension of the array AF.  LDAF >= max(1,N).
   53: *
   54: *     X       (input) COMPLEX*16 array, dimension (N)
   55: *     The vector X in the formula op(A) * diag(X).
   56: *
   57: *     INFO    (output) INTEGER
   58: *       = 0:  Successful exit.
   59: *     i > 0:  The ith argument is invalid.
   60: *
   61: *     WORK    (input) COMPLEX*16 array, dimension (2*N).
   62: *     Workspace.
   63: *
   64: *     RWORK   (input) DOUBLE PRECISION array, dimension (N).
   65: *     Workspace.
   66: *
   67: *  =====================================================================
   68: *
   69: *     .. Local Scalars ..
   70:       INTEGER            KASE, I, J
   71:       DOUBLE PRECISION   AINVNM, ANORM, TMP
   72:       LOGICAL            UP
   73:       COMPLEX*16         ZDUM
   74: *     ..
   75: *     .. Local Arrays ..
   76:       INTEGER            ISAVE( 3 )
   77: *     ..
   78: *     .. External Functions ..
   79:       LOGICAL            LSAME
   80:       EXTERNAL           LSAME
   81: *     ..
   82: *     .. External Subroutines ..
   83:       EXTERNAL           ZLACN2, ZPOTRS, XERBLA
   84: *     ..
   85: *     .. Intrinsic Functions ..
   86:       INTRINSIC          ABS, MAX, REAL, DIMAG
   87: *     ..
   88: *     .. Statement Functions ..
   89:       DOUBLE PRECISION CABS1
   90: *     ..
   91: *     .. Statement Function Definitions ..
   92:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
   93: *     ..
   94: *     .. Executable Statements ..
   95: *
   96:       ZLA_PORCOND_X = 0.0D+0
   97: *
   98:       INFO = 0
   99:       IF( N.LT.0 ) THEN
  100:          INFO = -2
  101:       END IF
  102:       IF( INFO.NE.0 ) THEN
  103:          CALL XERBLA( 'ZLA_PORCOND_X', -INFO )
  104:          RETURN
  105:       END IF
  106:       UP = .FALSE.
  107:       IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
  108: *
  109: *     Compute norm of op(A)*op2(C).
  110: *
  111:       ANORM = 0.0D+0
  112:       IF ( UP ) THEN
  113:          DO I = 1, N
  114:             TMP = 0.0D+0
  115:             DO J = 1, I
  116:                TMP = TMP + CABS1( A( J, I ) * X( J ) )
  117:             END DO
  118:             DO J = I+1, N
  119:                TMP = TMP + CABS1( A( I, J ) * X( J ) )
  120:             END DO
  121:             RWORK( I ) = TMP
  122:             ANORM = MAX( ANORM, TMP )
  123:          END DO
  124:       ELSE
  125:          DO I = 1, N
  126:             TMP = 0.0D+0
  127:             DO J = 1, I
  128:                TMP = TMP + CABS1( A( I, J ) * X( J ) )
  129:             END DO
  130:             DO J = I+1, N
  131:                TMP = TMP + CABS1( A( J, I ) * X( J ) )
  132:             END DO
  133:             RWORK( I ) = TMP
  134:             ANORM = MAX( ANORM, TMP )
  135:          END DO
  136:       END IF
  137: *
  138: *     Quick return if possible.
  139: *
  140:       IF( N.EQ.0 ) THEN
  141:          ZLA_PORCOND_X = 1.0D+0
  142:          RETURN
  143:       ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
  144:          RETURN
  145:       END IF
  146: *
  147: *     Estimate the norm of inv(op(A)).
  148: *
  149:       AINVNM = 0.0D+0
  150: *
  151:       KASE = 0
  152:    10 CONTINUE
  153:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
  154:       IF( KASE.NE.0 ) THEN
  155:          IF( KASE.EQ.2 ) THEN
  156: *
  157: *           Multiply by R.
  158: *
  159:             DO I = 1, N
  160:                WORK( I ) = WORK( I ) * RWORK( I )
  161:             END DO
  162: *
  163:             IF ( UP ) THEN
  164:                CALL ZPOTRS( 'U', N, 1, AF, LDAF,
  165:      $            WORK, N, INFO )
  166:             ELSE
  167:                CALL ZPOTRS( 'L', N, 1, AF, LDAF,
  168:      $            WORK, N, INFO )
  169:             ENDIF
  170: *
  171: *           Multiply by inv(X).
  172: *
  173:             DO I = 1, N
  174:                WORK( I ) = WORK( I ) / X( I )
  175:             END DO
  176:          ELSE
  177: *
  178: *           Multiply by inv(X').
  179: *
  180:             DO I = 1, N
  181:                WORK( I ) = WORK( I ) / X( I )
  182:             END DO
  183: *
  184:             IF ( UP ) THEN
  185:                CALL ZPOTRS( 'U', N, 1, AF, LDAF,
  186:      $            WORK, N, INFO )
  187:             ELSE
  188:                CALL ZPOTRS( 'L', N, 1, AF, LDAF,
  189:      $            WORK, N, INFO )
  190:             END IF
  191: *
  192: *           Multiply by R.
  193: *
  194:             DO I = 1, N
  195:                WORK( I ) = WORK( I ) * RWORK( I )
  196:             END DO
  197:          END IF
  198:          GO TO 10
  199:       END IF
  200: *
  201: *     Compute the estimate of the reciprocal condition number.
  202: *
  203:       IF( AINVNM .NE. 0.0D+0 )
  204:      $   ZLA_PORCOND_X = 1.0D+0 / AINVNM
  205: *
  206:       RETURN
  207: *
  208:       END