SRC/sla_gercond.f

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