SRC/cla_porcond_x.f

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