SRC/sla_syrcond.f

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