SRC/zla_gercond_x.f

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