SRC/zla_syrcond_x.f

   1 *> \brief \b ZLA_SYRCOND_X computes the infinity norm condition number of op(A)*diag(x) for symmetric 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_SYRCOND_X + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zla_syrcond_x.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zla_syrcond_x.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zla_syrcond_x.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       DOUBLE PRECISION FUNCTION ZLA_SYRCOND_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_SYRCOND_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 ZSYTRF.
  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 ZSYTRF.
  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 complex16SYcomputational
 130 *
 131 *  =====================================================================
 132       DOUBLE PRECISION FUNCTION ZLA_SYRCOND_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
 155       DOUBLE PRECISION   AINVNM, ANORM, TMP
 156       INTEGER            I, J
 157       LOGICAL            UP, UPPER
 158       COMPLEX*16         ZDUM
 159 *     ..
 160 *     .. Local Arrays ..
 161       INTEGER            ISAVE( 3 )
 162 *     ..
 163 *     .. External Functions ..
 164       LOGICAL            LSAME
 165       EXTERNAL           LSAME
 166 *     ..
 167 *     .. External Subroutines ..
 168       EXTERNAL           ZLACN2, ZSYTRS, XERBLA
 169 *     ..
 170 *     .. Intrinsic Functions ..
 171       INTRINSIC          ABS, MAX
 172 *     ..
 173 *     .. Statement Functions ..
 174       DOUBLE PRECISION   CABS1
 175 *     ..
 176 *     .. Statement Function Definitions ..
 177       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 178 *     ..
 179 *     .. Executable Statements ..
 180 *
 181       ZLA_SYRCOND_X = 0.0D+0
 182 *
 183       INFO = 0
 184       UPPER = LSAME( UPLO, 'U' )
 185       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 186          INFO = -1
 187       ELSE IF( N.LT.0 ) THEN
 188          INFO = -2
 189       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 190          INFO = -4
 191       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
 192          INFO = -6
 193       END IF
 194       IF( INFO.NE.0 ) THEN
 195          CALL XERBLA( 'ZLA_SYRCOND_X', -INFO )
 196          RETURN
 197       END IF
 198       UP = .FALSE.
 199       IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
 200 *
 201 *     Compute norm of op(A)*op2(C).
 202 *
 203       ANORM = 0.0D+0
 204       IF ( UP ) THEN
 205          DO I = 1, N
 206             TMP = 0.0D+0
 207             DO J = 1, I
 208                TMP = TMP + CABS1( A( J, I ) * X( J ) )
 209             END DO
 210             DO J = I+1, N
 211                TMP = TMP + CABS1( A( I, J ) * X( J ) )
 212             END DO
 213             RWORK( I ) = TMP
 214             ANORM = MAX( ANORM, TMP )
 215          END DO
 216       ELSE
 217          DO I = 1, N
 218             TMP = 0.0D+0
 219             DO J = 1, I
 220                TMP = TMP + CABS1( A( I, J ) * X( J ) )
 221             END DO
 222             DO J = I+1, N
 223                TMP = TMP + CABS1( A( J, I ) * X( J ) )
 224             END DO
 225             RWORK( I ) = TMP
 226             ANORM = MAX( ANORM, TMP )
 227          END DO
 228       END IF
 229 *
 230 *     Quick return if possible.
 231 *
 232       IF( N.EQ.0 ) THEN
 233          ZLA_SYRCOND_X = 1.0D+0
 234          RETURN
 235       ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
 236          RETURN
 237       END IF
 238 *
 239 *     Estimate the norm of inv(op(A)).
 240 *
 241       AINVNM = 0.0D+0
 242 *
 243       KASE = 0
 244    10 CONTINUE
 245       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
 246       IF( KASE.NE.0 ) THEN
 247          IF( KASE.EQ.2 ) THEN
 248 *
 249 *           Multiply by R.
 250 *
 251             DO I = 1, N
 252                WORK( I ) = WORK( I ) * RWORK( I )
 253             END DO
 254 *
 255             IF ( UP ) THEN
 256                CALL ZSYTRS( 'U', N, 1, AF, LDAF, IPIV,
 257      $            WORK, N, INFO )
 258             ELSE
 259                CALL ZSYTRS( 'L', N, 1, AF, LDAF, IPIV,
 260      $            WORK, N, INFO )
 261             ENDIF
 262 *
 263 *           Multiply by inv(X).
 264 *
 265             DO I = 1, N
 266                WORK( I ) = WORK( I ) / X( I )
 267             END DO
 268          ELSE
 269 *
 270 *           Multiply by inv(X**T).
 271 *
 272             DO I = 1, N
 273                WORK( I ) = WORK( I ) / X( I )
 274             END DO
 275 *
 276             IF ( UP ) THEN
 277                CALL ZSYTRS( 'U', N, 1, AF, LDAF, IPIV,
 278      $            WORK, N, INFO )
 279             ELSE
 280                CALL ZSYTRS( 'L', N, 1, AF, LDAF, IPIV,
 281      $            WORK, N, INFO )
 282             END IF
 283 *
 284 *           Multiply by R.
 285 *
 286             DO I = 1, N
 287                WORK( I ) = WORK( I ) * RWORK( I )
 288             END DO
 289          END IF
 290          GO TO 10
 291       END IF
 292 *
 293 *     Compute the estimate of the reciprocal condition number.
 294 *
 295       IF( AINVNM .NE. 0.0D+0 )
 296      $   ZLA_SYRCOND_X = 1.0D+0 / AINVNM
 297 *
 298       RETURN
 299 *
 300       END