SRC/zla_gercond_c.f

   1       DOUBLE PRECISION FUNCTION ZLA_GERCOND_C( TRANS, N, A, LDA, AF,
   2      $                                         LDAF, IPIV, C, CAPPLY,
   3      $                                         INFO, WORK, 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 Aguments ..
  16       CHARACTER          TRANS
  17       LOGICAL            CAPPLY
  18       INTEGER            N, LDA, LDAF, INFO
  19 *     ..
  20 *     .. Array Arguments ..
  21       INTEGER            IPIV( * )
  22       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), WORK( * )
  23       DOUBLE PRECISION   C( * ), RWORK( * )
  24 *     ..
  25 *
  26 *  Purpose
  27 *  =======
  28 *
  29 *     ZLA_GERCOND_C computes the infinity norm condition number of
  30 *     op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
  31 *
  32 *  Arguments
  33 *  =========
  34 *
  35 *     TRANS   (input) CHARACTER*1
  36 *     Specifies the form of the system of equations:
  37 *       = 'N':  A * X = B     (No transpose)
  38 *       = 'T':  A**T * X = B  (Transpose)
  39 *       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose)
  40 *
  41 *     N       (input) INTEGER
  42 *     The number of linear equations, i.e., the order of the
  43 *     matrix A.  N >= 0.
  44 *
  45 *     A       (input) COMPLEX*16 array, dimension (LDA,N)
  46 *     On entry, the N-by-N matrix A
  47 *
  48 *     LDA     (input) INTEGER
  49 *     The leading dimension of the array A.  LDA >= max(1,N).
  50 *
  51 *     AF      (input) COMPLEX*16 array, dimension (LDAF,N)
  52 *     The factors L and U from the factorization
  53 *     A = P*L*U as computed by ZGETRF.
  54 *
  55 *     LDAF    (input) INTEGER
  56 *     The leading dimension of the array AF.  LDAF >= max(1,N).
  57 *
  58 *     IPIV    (input) INTEGER array, dimension (N)
  59 *     The pivot indices from the factorization A = P*L*U
  60 *     as computed by ZGETRF; row i of the matrix was interchanged
  61 *     with row IPIV(i).
  62 *
  63 *     C       (input) DOUBLE PRECISION array, dimension (N)
  64 *     The vector C in the formula op(A) * inv(diag(C)).
  65 *
  66 *     CAPPLY  (input) LOGICAL
  67 *     If .TRUE. then access the vector C in the formula above.
  68 *
  69 *     INFO    (output) INTEGER
  70 *       = 0:  Successful exit.
  71 *     i > 0:  The ith argument is invalid.
  72 *
  73 *     WORK    (input) COMPLEX*16 array, dimension (2*N).
  74 *     Workspace.
  75 *
  76 *     RWORK   (input) DOUBLE PRECISION array, dimension (N).
  77 *     Workspace.
  78 *
  79 *  =====================================================================
  80 *
  81 *     .. Local Scalars ..
  82       LOGICAL            NOTRANS
  83       INTEGER            KASE, I, J
  84       DOUBLE PRECISION   AINVNM, ANORM, TMP
  85       COMPLEX*16         ZDUM
  86 *     ..
  87 *     .. Local Arrays ..
  88       INTEGER            ISAVE( 3 )
  89 *     ..
  90 *     .. External Functions ..
  91       LOGICAL            LSAME
  92       EXTERNAL           LSAME
  93 *     ..
  94 *     .. External Subroutines ..
  95       EXTERNAL           ZLACN2, ZGETRS, XERBLA
  96 *     ..
  97 *     .. Intrinsic Functions ..
  98       INTRINSIC          ABS, MAX, REAL, DIMAG
  99 *     ..
 100 *     .. Statement Functions ..
 101       DOUBLE PRECISION   CABS1
 102 *     ..
 103 *     .. Statement Function Definitions ..
 104       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 105 *     ..
 106 *     .. Executable Statements ..
 107       ZLA_GERCOND_C = 0.0D+0
 108 *
 109       INFO = 0
 110       NOTRANS = LSAME( TRANS, 'N' )
 111       IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
 112      $     LSAME( TRANS, 'C' ) ) THEN
 113       ELSE IF( N.LT.0 ) THEN
 114          INFO = -2
 115       END IF
 116       IF( INFO.NE.0 ) THEN
 117          CALL XERBLA( 'ZLA_GERCOND_C', -INFO )
 118          RETURN
 119       END IF
 120 *
 121 *     Compute norm of op(A)*op2(C).
 122 *
 123       ANORM = 0.0D+0
 124       IF ( NOTRANS ) THEN
 125          DO I = 1, N
 126             TMP = 0.0D+0
 127             IF ( CAPPLY ) THEN
 128                DO J = 1, N
 129                   TMP = TMP + CABS1( A( I, J ) ) / C( J )
 130                END DO
 131             ELSE
 132                DO J = 1, N
 133                   TMP = TMP + CABS1( A( I, J ) )
 134                END DO
 135             END IF
 136             RWORK( I ) = TMP
 137             ANORM = MAX( ANORM, TMP )
 138          END DO
 139       ELSE
 140          DO I = 1, N
 141             TMP = 0.0D+0
 142             IF ( CAPPLY ) THEN
 143                DO J = 1, N
 144                   TMP = TMP + CABS1( A( J, I ) ) / C( J )
 145                END DO
 146             ELSE
 147                DO J = 1, N
 148                   TMP = TMP + CABS1( A( J, I ) )
 149                END DO
 150             END IF
 151             RWORK( I ) = TMP
 152             ANORM = MAX( ANORM, TMP )
 153          END DO
 154       END IF
 155 *
 156 *     Quick return if possible.
 157 *
 158       IF( N.EQ.0 ) THEN
 159          ZLA_GERCOND_C = 1.0D+0
 160          RETURN
 161       ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
 162          RETURN
 163       END IF
 164 *
 165 *     Estimate the norm of inv(op(A)).
 166 *
 167       AINVNM = 0.0D+0
 168 *
 169       KASE = 0
 170    10 CONTINUE
 171       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
 172       IF( KASE.NE.0 ) THEN
 173          IF( KASE.EQ.2 ) THEN
 174 *
 175 *           Multiply by R.
 176 *
 177             DO I = 1, N
 178                WORK( I ) = WORK( I ) * RWORK( I )
 179             END DO
 180 *
 181             IF (NOTRANS) THEN
 182                CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
 183      $            WORK, N, INFO )
 184             ELSE
 185                CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
 186      $            WORK, N, INFO )
 187             ENDIF
 188 *
 189 *           Multiply by inv(C).
 190 *
 191             IF ( CAPPLY ) THEN
 192                DO I = 1, N
 193                   WORK( I ) = WORK( I ) * C( I )
 194                END DO
 195             END IF
 196          ELSE
 197 *
 198 *           Multiply by inv(C**H).
 199 *
 200             IF ( CAPPLY ) THEN
 201                DO I = 1, N
 202                   WORK( I ) = WORK( I ) * C( I )
 203                END DO
 204             END IF
 205 *
 206             IF ( NOTRANS ) THEN
 207                CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
 208      $            WORK, N, INFO )
 209             ELSE
 210                CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
 211      $            WORK, N, INFO )
 212             END IF
 213 *
 214 *           Multiply by R.
 215 *
 216             DO I = 1, N
 217                WORK( I ) = WORK( I ) * RWORK( I )
 218             END DO
 219          END IF
 220          GO TO 10
 221       END IF
 222 *
 223 *     Compute the estimate of the reciprocal condition number.
 224 *
 225       IF( AINVNM .NE. 0.0D+0 )
 226      $   ZLA_GERCOND_C = 1.0D+0 / AINVNM
 227 *
 228       RETURN
 229 *
 230       END