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