SRC/zla_gercond_c.f

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