SRC/cla_hercond_x.f

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