SRC/cla_syrcond_c.f

   1 *> \brief \b CLA_SYRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) 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 CLA_SYRCOND_C + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cla_syrcond_c.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cla_syrcond_c.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cla_syrcond_c.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       REAL FUNCTION CLA_SYRCOND_C( UPLO, N, A, LDA, AF, LDAF, IPIV, C,
  22 *                                    CAPPLY, INFO, WORK, RWORK )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       LOGICAL            CAPPLY
  27 *       INTEGER            N, LDA, LDAF, INFO
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       INTEGER            IPIV( * )
  31 *       COMPLEX            A( LDA, * ), AF( LDAF, * ), WORK( * )
  32 *       REAL               C( * ), RWORK( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *>    CLA_SYRCOND_C Computes the infinity norm condition number of
  42 *>    op(A) * inv(diag(C)) where C is a REAL 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 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 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 CSYTRF.
  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 CSYTRF.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] C
  95 *> \verbatim
  96 *>          C is REAL array, dimension (N)
  97 *>     The vector C in the formula op(A) * inv(diag(C)).
  98 *> \endverbatim
  99 *>
 100 *> \param[in] CAPPLY
 101 *> \verbatim
 102 *>          CAPPLY is LOGICAL
 103 *>     If .TRUE. then access the vector C in the formula above.
 104 *> \endverbatim
 105 *>
 106 *> \param[out] INFO
 107 *> \verbatim
 108 *>          INFO is INTEGER
 109 *>       = 0:  Successful exit.
 110 *>     i > 0:  The ith argument is invalid.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] WORK
 114 *> \verbatim
 115 *>          WORK is COMPLEX array, dimension (2*N).
 116 *>     Workspace.
 117 *> \endverbatim
 118 *>
 119 *> \param[in] RWORK
 120 *> \verbatim
 121 *>          RWORK is REAL array, dimension (N).
 122 *>     Workspace.
 123 *> \endverbatim
 124 *
 125 *  Authors:
 126 *  ========
 127 *
 128 *> \author Univ. of Tennessee
 129 *> \author Univ. of California Berkeley
 130 *> \author Univ. of Colorado Denver
 131 *> \author NAG Ltd.
 132 *
 133 *> \date September 2012
 134 *
 135 *> \ingroup complexSYcomputational
 136 *
 137 *  =====================================================================
 138       REAL FUNCTION CLA_SYRCOND_C( UPLO, N, A, LDA, AF, LDAF, IPIV, C,
 139      $                             CAPPLY, INFO, WORK, RWORK )
 140 *
 141 *  -- LAPACK computational routine (version 3.4.2) --
 142 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 143 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 144 *     September 2012
 145 *
 146 *     .. Scalar Arguments ..
 147       CHARACTER          UPLO
 148       LOGICAL            CAPPLY
 149       INTEGER            N, LDA, LDAF, INFO
 150 *     ..
 151 *     .. Array Arguments ..
 152       INTEGER            IPIV( * )
 153       COMPLEX            A( LDA, * ), AF( LDAF, * ), WORK( * )
 154       REAL               C( * ), RWORK( * )
 155 *     ..
 156 *
 157 *  =====================================================================
 158 *
 159 *     .. Local Scalars ..
 160       INTEGER            KASE
 161       REAL               AINVNM, ANORM, TMP
 162       INTEGER            I, J
 163       LOGICAL            UP, UPPER
 164       COMPLEX            ZDUM
 165 *     ..
 166 *     .. Local Arrays ..
 167       INTEGER            ISAVE( 3 )
 168 *     ..
 169 *     .. External Functions ..
 170       LOGICAL            LSAME
 171       EXTERNAL           LSAME
 172 *     ..
 173 *     .. External Subroutines ..
 174       EXTERNAL           CLACN2, CSYTRS, XERBLA
 175 *     ..
 176 *     .. Intrinsic Functions ..
 177       INTRINSIC          ABS, MAX
 178 *     ..
 179 *     .. Statement Functions ..
 180       REAL CABS1
 181 *     ..
 182 *     .. Statement Function Definitions ..
 183       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
 184 *     ..
 185 *     .. Executable Statements ..
 186 *
 187       CLA_SYRCOND_C = 0.0E+0
 188 *
 189       INFO = 0
 190       UPPER = LSAME( UPLO, 'U' )
 191       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 192          INFO = -1
 193       ELSE IF( N.LT.0 ) THEN
 194          INFO = -2
 195       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 196          INFO = -4
 197       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
 198          INFO = -6
 199       END IF
 200       IF( INFO.NE.0 ) THEN
 201          CALL XERBLA( 'CLA_SYRCOND_C', -INFO )
 202          RETURN
 203       END IF
 204       UP = .FALSE.
 205       IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
 206 *
 207 *     Compute norm of op(A)*op2(C).
 208 *
 209       ANORM = 0.0E+0
 210       IF ( UP ) THEN
 211          DO I = 1, N
 212             TMP = 0.0E+0
 213             IF ( CAPPLY ) THEN
 214                DO J = 1, I
 215                   TMP = TMP + CABS1( A( J, I ) ) / C( J )
 216                END DO
 217                DO J = I+1, N
 218                   TMP = TMP + CABS1( A( I, J ) ) / C( J )
 219                END DO
 220             ELSE
 221                DO J = 1, I
 222                   TMP = TMP + CABS1( A( J, I ) )
 223                END DO
 224                DO J = I+1, N
 225                   TMP = TMP + CABS1( A( I, J ) )
 226                END DO
 227             END IF
 228             RWORK( I ) = TMP
 229             ANORM = MAX( ANORM, TMP )
 230          END DO
 231       ELSE
 232          DO I = 1, N
 233             TMP = 0.0E+0
 234             IF ( CAPPLY ) THEN
 235                DO J = 1, I
 236                   TMP = TMP + CABS1( A( I, J ) ) / C( J )
 237                END DO
 238                DO J = I+1, N
 239                   TMP = TMP + CABS1( A( J, I ) ) / C( J )
 240                END DO
 241             ELSE
 242                DO J = 1, I
 243                   TMP = TMP + CABS1( A( I, J ) )
 244                END DO
 245                DO J = I+1, N
 246                   TMP = TMP + CABS1( A( J, I ) )
 247                END DO
 248             END IF
 249             RWORK( I ) = TMP
 250             ANORM = MAX( ANORM, TMP )
 251          END DO
 252       END IF
 253 *
 254 *     Quick return if possible.
 255 *
 256       IF( N.EQ.0 ) THEN
 257          CLA_SYRCOND_C = 1.0E+0
 258          RETURN
 259       ELSE IF( ANORM .EQ. 0.0E+0 ) THEN
 260          RETURN
 261       END IF
 262 *
 263 *     Estimate the norm of inv(op(A)).
 264 *
 265       AINVNM = 0.0E+0
 266 *
 267       KASE = 0
 268    10 CONTINUE
 269       CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
 270       IF( KASE.NE.0 ) THEN
 271          IF( KASE.EQ.2 ) THEN
 272 *
 273 *           Multiply by R.
 274 *
 275             DO I = 1, N
 276                WORK( I ) = WORK( I ) * RWORK( I )
 277             END DO
 278 *
 279             IF ( UP ) THEN
 280                CALL CSYTRS( 'U', N, 1, AF, LDAF, IPIV,
 281      $            WORK, N, INFO )
 282             ELSE
 283                CALL CSYTRS( 'L', N, 1, AF, LDAF, IPIV,
 284      $            WORK, N, INFO )
 285             ENDIF
 286 *
 287 *           Multiply by inv(C).
 288 *
 289             IF ( CAPPLY ) THEN
 290                DO I = 1, N
 291                   WORK( I ) = WORK( I ) * C( I )
 292                END DO
 293             END IF
 294          ELSE
 295 *
 296 *           Multiply by inv(C**T).
 297 *
 298             IF ( CAPPLY ) THEN
 299                DO I = 1, N
 300                   WORK( I ) = WORK( I ) * C( I )
 301                END DO
 302             END IF
 303 *
 304             IF ( UP ) THEN
 305                CALL CSYTRS( 'U', N, 1, AF, LDAF, IPIV,
 306      $            WORK, N, INFO )
 307             ELSE
 308                CALL CSYTRS( 'L', N, 1, AF, LDAF, IPIV,
 309      $            WORK, N, INFO )
 310             END IF
 311 *
 312 *           Multiply by R.
 313 *
 314             DO I = 1, N
 315                WORK( I ) = WORK( I ) * RWORK( I )
 316             END DO
 317          END IF
 318          GO TO 10
 319       END IF
 320 *
 321 *     Compute the estimate of the reciprocal condition number.
 322 *
 323       IF( AINVNM .NE. 0.0E+0 )
 324      $   CLA_SYRCOND_C = 1.0E+0 / AINVNM
 325 *
 326       RETURN
 327 *
 328       END