SRC/zla_hercond_x.f

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