SRC/zla_gercond_c.f

   1       DOUBLE PRECISION FUNCTION ZLA_GERCOND_C( TRANS, N, A, LDA, AF,
   2      $                             LDAF, IPIV, C, CAPPLY, INFO, WORK,
   3      $     RWORK )
   4 *
   5 *     -- LAPACK routine (version 3.2)                                 --
   6 *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
   7 *     -- Jason Riedy of Univ. of California Berkeley.                 --
   8 *     -- November 2008                                                --
   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 *     ZLA_GERCOND_C computes the infinity norm condition number of
  26 *     op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
  27 *     WORK is a COMPLEX*16 workspace of size 2*N, and
  28 *     RWORK is a DOUBLE PRECISION workspace of size 3*N.
  29 *     ..
  30 *     .. Local Scalars ..
  31       LOGICAL            NOTRANS
  32       INTEGER            KASE, I, J
  33       DOUBLE PRECISION   AINVNM, ANORM, TMP
  34       COMPLEX*16         ZDUM
  35 *     ..
  36 *     .. Local Arrays ..
  37       INTEGER            ISAVE( 3 )
  38 *     ..
  39 *     .. External Functions ..
  40       LOGICAL            LSAME
  41       EXTERNAL           LSAME
  42 *     ..
  43 *     .. External Subroutines ..
  44       EXTERNAL           ZLACN2, ZGETRS, XERBLA
  45 *     ..
  46 *     .. Intrinsic Functions ..
  47       INTRINSIC          ABS, MAX, REAL, DIMAG
  48 *     ..
  49 *     .. Statement Functions ..
  50       DOUBLE PRECISION   CABS1
  51 *     ..
  52 *     .. Statement Function Definitions ..
  53       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
  54 *     ..
  55 *     .. Executable Statements ..
  56       ZLA_GERCOND_C = 0.0D+0
  57 *
  58       INFO = 0
  59       NOTRANS = LSAME( TRANS, 'N' )
  60       IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
  61      $     LSAME( TRANS, 'C' ) ) THEN
  62       ELSE IF( N.LT.0 ) THEN
  63          INFO = -2
  64       END IF
  65       IF( INFO.NE.0 ) THEN
  66          CALL XERBLA( 'ZLA_GERCOND_C', -INFO )
  67          RETURN
  68       END IF
  69 *
  70 *     Compute norm of op(A)*op2(C).
  71 *
  72       ANORM = 0.0D+0
  73       IF ( NOTRANS ) THEN
  74          DO I = 1, N
  75             TMP = 0.0D+0
  76             IF ( CAPPLY ) THEN
  77                DO J = 1, N
  78                   TMP = TMP + CABS1( A( I, J ) ) / C( J )
  79                END DO
  80             ELSE
  81                DO J = 1, N
  82                   TMP = TMP + CABS1( A( I, J ) )
  83                END DO
  84             END IF
  85             RWORK( 2*N+I ) = TMP
  86             ANORM = MAX( ANORM, TMP )
  87          END DO
  88       ELSE
  89          DO I = 1, N
  90             TMP = 0.0D+0
  91             IF ( CAPPLY ) THEN
  92                DO J = 1, N
  93                   TMP = TMP + CABS1( A( J, I ) ) / C( J )
  94                END DO
  95             ELSE
  96                DO J = 1, N
  97                   TMP = TMP + CABS1( A( J, I ) )
  98                END DO
  99             END IF
 100             RWORK( 2*N+I ) = TMP
 101             ANORM = MAX( ANORM, TMP )
 102          END DO
 103       END IF
 104 *
 105 *     Quick return if possible.
 106 *
 107       IF( N.EQ.0 ) THEN
 108          ZLA_GERCOND_C = 1.0D+0
 109          RETURN
 110       ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
 111          RETURN
 112       END IF
 113 *
 114 *     Estimate the norm of inv(op(A)).
 115 *
 116       AINVNM = 0.0D+0
 117 *
 118       KASE = 0
 119    10 CONTINUE
 120       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
 121       IF( KASE.NE.0 ) THEN
 122          IF( KASE.EQ.2 ) THEN
 123 *
 124 *           Multiply by R.
 125 *
 126             DO I = 1, N
 127                WORK( I ) = WORK( I ) * RWORK( 2*N+I )
 128             END DO
 129 *
 130             IF (NOTRANS) THEN
 131                CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
 132      $            WORK, N, INFO )
 133             ELSE
 134                CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
 135      $            WORK, N, INFO )
 136             ENDIF
 137 *
 138 *           Multiply by inv(C).
 139 *
 140             IF ( CAPPLY ) THEN
 141                DO I = 1, N
 142                   WORK( I ) = WORK( I ) * C( I )
 143                END DO
 144             END IF
 145          ELSE
 146 *
 147 *           Multiply by inv(C').
 148 *
 149             IF ( CAPPLY ) THEN
 150                DO I = 1, N
 151                   WORK( I ) = WORK( I ) * C( I )
 152                END DO
 153             END IF
 154 *
 155             IF ( NOTRANS ) THEN
 156                CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
 157      $            WORK, N, INFO )
 158             ELSE
 159                CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
 160      $            WORK, N, INFO )
 161             END IF
 162 *
 163 *           Multiply by R.
 164 *
 165             DO I = 1, N
 166                WORK( I ) = WORK( I ) * RWORK( 2*N+I )
 167             END DO
 168          END IF
 169          GO TO 10
 170       END IF
 171 *
 172 *     Compute the estimate of the reciprocal condition number.
 173 *
 174       IF( AINVNM .NE. 0.0D+0 )
 175      $   ZLA_GERCOND_C = 1.0D+0 / AINVNM
 176 *
 177       RETURN
 178 *
 179       END