SRC/dla_gercond.f

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