SRC/zla_gercond_c.f

   1 *> \brief \b ZLA_GERCOND_C
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition
   9 *  ==========
  10 *
  11 *       DOUBLE PRECISION FUNCTION ZLA_GERCOND_C( TRANS, N, A, LDA, AF,
  12 *                                                LDAF, IPIV, C, CAPPLY,
  13 *                                                INFO, WORK, RWORK )
  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 *>\details \b Purpose:
  30 *>\verbatim
  31 *>
  32 *>    ZLA_GERCOND_C computes the infinity norm condition number of
  33 *>    op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
  34 *>
  35 *>\endverbatim
  36 *
  37 *  Arguments
  38 *  =========
  39 *
  40 *> \param[in] TRANS
  41 *> \verbatim
  42 *>          TRANS is CHARACTER*1
  43 *>     Specifies the form of the system of equations:
  44 *>       = 'N':  A * X = B     (No transpose)
  45 *>       = 'T':  A**T * X = B  (Transpose)
  46 *>       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose)
  47 *> \endverbatim
  48 *>
  49 *> \param[in] N
  50 *> \verbatim
  51 *>          N is INTEGER
  52 *>     The number of linear equations, i.e., the order of the
  53 *>     matrix A.  N >= 0.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] A
  57 *> \verbatim
  58 *>          A is COMPLEX*16 array, dimension (LDA,N)
  59 *>     On entry, the N-by-N matrix A
  60 *> \endverbatim
  61 *>
  62 *> \param[in] LDA
  63 *> \verbatim
  64 *>          LDA is INTEGER
  65 *>     The leading dimension of the array A.  LDA >= max(1,N).
  66 *> \endverbatim
  67 *>
  68 *> \param[in] AF
  69 *> \verbatim
  70 *>          AF is COMPLEX*16 array, dimension (LDAF,N)
  71 *>     The factors L and U from the factorization
  72 *>     A = P*L*U as computed by ZGETRF.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] LDAF
  76 *> \verbatim
  77 *>          LDAF is INTEGER
  78 *>     The leading dimension of the array AF.  LDAF >= max(1,N).
  79 *> \endverbatim
  80 *>
  81 *> \param[in] IPIV
  82 *> \verbatim
  83 *>          IPIV is INTEGER array, dimension (N)
  84 *>     The pivot indices from the factorization A = P*L*U
  85 *>     as computed by ZGETRF; row i of the matrix was interchanged
  86 *>     with row IPIV(i).
  87 *> \endverbatim
  88 *>
  89 *> \param[in] C
  90 *> \verbatim
  91 *>          C is DOUBLE PRECISION array, dimension (N)
  92 *>     The vector C in the formula op(A) * inv(diag(C)).
  93 *> \endverbatim
  94 *>
  95 *> \param[in] CAPPLY
  96 *> \verbatim
  97 *>          CAPPLY is LOGICAL
  98 *>     If .TRUE. then access the vector C in the formula above.
  99 *> \endverbatim
 100 *>
 101 *> \param[out] INFO
 102 *> \verbatim
 103 *>          INFO is INTEGER
 104 *>       = 0:  Successful exit.
 105 *>     i > 0:  The ith argument is invalid.
 106 *> \endverbatim
 107 *>
 108 *> \param[in] WORK
 109 *> \verbatim
 110 *>          WORK is COMPLEX*16 array, dimension (2*N).
 111 *>     Workspace.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] RWORK
 115 *> \verbatim
 116 *>          RWORK is DOUBLE PRECISION array, dimension (N).
 117 *>     Workspace.
 118 *> \endverbatim
 119 *>
 120 *
 121 *  Authors
 122 *  =======
 123 *
 124 *> \author Univ. of Tennessee
 125 *> \author Univ. of California Berkeley
 126 *> \author Univ. of Colorado Denver
 127 *> \author NAG Ltd.
 128 *
 129 *> \date November 2011
 130 *
 131 *> \ingroup complex16GEcomputational
 132 *
 133 *  =====================================================================
 134       DOUBLE PRECISION FUNCTION ZLA_GERCOND_C( TRANS, N, A, LDA, AF,
 135      $                                         LDAF, IPIV, C, CAPPLY,
 136      $                                         INFO, WORK, RWORK )
 137 *
 138 *  -- LAPACK computational routine (version 3.2.1) --
 139 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 140 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 141 *     November 2011
 142 *
 143 *     .. Scalar Aguments ..
 144       CHARACTER          TRANS
 145       LOGICAL            CAPPLY
 146       INTEGER            N, LDA, LDAF, INFO
 147 *     ..
 148 *     .. Array Arguments ..
 149       INTEGER            IPIV( * )
 150       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), WORK( * )
 151       DOUBLE PRECISION   C( * ), RWORK( * )
 152 *     ..
 153 *
 154 *  =====================================================================
 155 *
 156 *     .. Local Scalars ..
 157       LOGICAL            NOTRANS
 158       INTEGER            KASE, I, J
 159       DOUBLE PRECISION   AINVNM, ANORM, TMP
 160       COMPLEX*16         ZDUM
 161 *     ..
 162 *     .. Local Arrays ..
 163       INTEGER            ISAVE( 3 )
 164 *     ..
 165 *     .. External Functions ..
 166       LOGICAL            LSAME
 167       EXTERNAL           LSAME
 168 *     ..
 169 *     .. External Subroutines ..
 170       EXTERNAL           ZLACN2, ZGETRS, XERBLA
 171 *     ..
 172 *     .. Intrinsic Functions ..
 173       INTRINSIC          ABS, MAX, REAL, DIMAG
 174 *     ..
 175 *     .. Statement Functions ..
 176       DOUBLE PRECISION   CABS1
 177 *     ..
 178 *     .. Statement Function Definitions ..
 179       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 180 *     ..
 181 *     .. Executable Statements ..
 182       ZLA_GERCOND_C = 0.0D+0
 183 *
 184       INFO = 0
 185       NOTRANS = LSAME( TRANS, 'N' )
 186       IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
 187      $     LSAME( TRANS, 'C' ) ) THEN
 188       ELSE IF( N.LT.0 ) THEN
 189          INFO = -2
 190       END IF
 191       IF( INFO.NE.0 ) THEN
 192          CALL XERBLA( 'ZLA_GERCOND_C', -INFO )
 193          RETURN
 194       END IF
 195 *
 196 *     Compute norm of op(A)*op2(C).
 197 *
 198       ANORM = 0.0D+0
 199       IF ( NOTRANS ) THEN
 200          DO I = 1, N
 201             TMP = 0.0D+0
 202             IF ( CAPPLY ) THEN
 203                DO J = 1, N
 204                   TMP = TMP + CABS1( A( I, J ) ) / C( J )
 205                END DO
 206             ELSE
 207                DO J = 1, N
 208                   TMP = TMP + CABS1( A( I, J ) )
 209                END DO
 210             END IF
 211             RWORK( I ) = TMP
 212             ANORM = MAX( ANORM, TMP )
 213          END DO
 214       ELSE
 215          DO I = 1, N
 216             TMP = 0.0D+0
 217             IF ( CAPPLY ) THEN
 218                DO J = 1, N
 219                   TMP = TMP + CABS1( A( J, I ) ) / C( J )
 220                END DO
 221             ELSE
 222                DO J = 1, N
 223                   TMP = TMP + CABS1( A( J, I ) )
 224                END DO
 225             END IF
 226             RWORK( I ) = TMP
 227             ANORM = MAX( ANORM, TMP )
 228          END DO
 229       END IF
 230 *
 231 *     Quick return if possible.
 232 *
 233       IF( N.EQ.0 ) THEN
 234          ZLA_GERCOND_C = 1.0D+0
 235          RETURN
 236       ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
 237          RETURN
 238       END IF
 239 *
 240 *     Estimate the norm of inv(op(A)).
 241 *
 242       AINVNM = 0.0D+0
 243 *
 244       KASE = 0
 245    10 CONTINUE
 246       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
 247       IF( KASE.NE.0 ) THEN
 248          IF( KASE.EQ.2 ) THEN
 249 *
 250 *           Multiply by R.
 251 *
 252             DO I = 1, N
 253                WORK( I ) = WORK( I ) * RWORK( I )
 254             END DO
 255 *
 256             IF (NOTRANS) THEN
 257                CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
 258      $            WORK, N, INFO )
 259             ELSE
 260                CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
 261      $            WORK, N, INFO )
 262             ENDIF
 263 *
 264 *           Multiply by inv(C).
 265 *
 266             IF ( CAPPLY ) THEN
 267                DO I = 1, N
 268                   WORK( I ) = WORK( I ) * C( I )
 269                END DO
 270             END IF
 271          ELSE
 272 *
 273 *           Multiply by inv(C**H).
 274 *
 275             IF ( CAPPLY ) THEN
 276                DO I = 1, N
 277                   WORK( I ) = WORK( I ) * C( I )
 278                END DO
 279             END IF
 280 *
 281             IF ( NOTRANS ) THEN
 282                CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
 283      $            WORK, N, INFO )
 284             ELSE
 285                CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
 286      $            WORK, N, INFO )
 287             END IF
 288 *
 289 *           Multiply by R.
 290 *
 291             DO I = 1, N
 292                WORK( I ) = WORK( I ) * RWORK( I )
 293             END DO
 294          END IF
 295          GO TO 10
 296       END IF
 297 *
 298 *     Compute the estimate of the reciprocal condition number.
 299 *
 300       IF( AINVNM .NE. 0.0D+0 )
 301      $   ZLA_GERCOND_C = 1.0D+0 / AINVNM
 302 *
 303       RETURN
 304 *
 305       END