LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sla_gercond.f
Go to the documentation of this file.
00001 *> \brief \b SLA_GERCOND
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLA_GERCOND + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sla_gercond.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sla_gercond.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sla_gercond.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       REAL FUNCTION SLA_GERCOND ( TRANS, N, A, LDA, AF, LDAF, IPIV,
00022 *                                   CMODE, C, INFO, WORK, IWORK )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          TRANS
00026 *       INTEGER            N, LDA, LDAF, INFO, CMODE
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IPIV( * ), IWORK( * )
00030 *       REAL               A( LDA, * ), AF( LDAF, * ), WORK( * ),
00031 *      $                   C( * )
00032 *      ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *>    SLA_GERCOND estimates the Skeel condition number of op(A) * op2(C)
00041 *>    where op2 is determined by CMODE as follows
00042 *>    CMODE =  1    op2(C) = C
00043 *>    CMODE =  0    op2(C) = I
00044 *>    CMODE = -1    op2(C) = inv(C)
00045 *>    The Skeel condition number cond(A) = norminf( |inv(A)||A| )
00046 *>    is computed by computing scaling factors R such that
00047 *>    diag(R)*A*op2(C) is row equilibrated and computing the standard
00048 *>    infinity-norm condition number.
00049 *> \endverbatim
00050 *
00051 *  Arguments:
00052 *  ==========
00053 *
00054 *> \param[in] TRANS
00055 *> \verbatim
00056 *>          TRANS is CHARACTER*1
00057 *>     Specifies the form of the system of equations:
00058 *>       = 'N':  A * X = B     (No transpose)
00059 *>       = 'T':  A**T * X = B  (Transpose)
00060 *>       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose)
00061 *> \endverbatim
00062 *>
00063 *> \param[in] N
00064 *> \verbatim
00065 *>          N is INTEGER
00066 *>     The number of linear equations, i.e., the order of the
00067 *>     matrix A.  N >= 0.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] A
00071 *> \verbatim
00072 *>          A is REAL array, dimension (LDA,N)
00073 *>     On entry, the N-by-N matrix A.
00074 *> \endverbatim
00075 *>
00076 *> \param[in] LDA
00077 *> \verbatim
00078 *>          LDA is INTEGER
00079 *>     The leading dimension of the array A.  LDA >= max(1,N).
00080 *> \endverbatim
00081 *>
00082 *> \param[in] AF
00083 *> \verbatim
00084 *>          AF is REAL array, dimension (LDAF,N)
00085 *>     The factors L and U from the factorization
00086 *>     A = P*L*U as computed by SGETRF.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] LDAF
00090 *> \verbatim
00091 *>          LDAF is INTEGER
00092 *>     The leading dimension of the array AF.  LDAF >= max(1,N).
00093 *> \endverbatim
00094 *>
00095 *> \param[in] IPIV
00096 *> \verbatim
00097 *>          IPIV is INTEGER array, dimension (N)
00098 *>     The pivot indices from the factorization A = P*L*U
00099 *>     as computed by SGETRF; row i of the matrix was interchanged
00100 *>     with row IPIV(i).
00101 *> \endverbatim
00102 *>
00103 *> \param[in] CMODE
00104 *> \verbatim
00105 *>          CMODE is INTEGER
00106 *>     Determines op2(C) in the formula op(A) * op2(C) as follows:
00107 *>     CMODE =  1    op2(C) = C
00108 *>     CMODE =  0    op2(C) = I
00109 *>     CMODE = -1    op2(C) = inv(C)
00110 *> \endverbatim
00111 *>
00112 *> \param[in] C
00113 *> \verbatim
00114 *>          C is REAL array, dimension (N)
00115 *>     The vector C in the formula op(A) * op2(C).
00116 *> \endverbatim
00117 *>
00118 *> \param[out] INFO
00119 *> \verbatim
00120 *>          INFO is INTEGER
00121 *>       = 0:  Successful exit.
00122 *>     i > 0:  The ith argument is invalid.
00123 *> \endverbatim
00124 *>
00125 *> \param[in] WORK
00126 *> \verbatim
00127 *>          WORK is REAL array, dimension (3*N).
00128 *>     Workspace.
00129 *> \endverbatim
00130 *>
00131 *> \param[in] IWORK
00132 *> \verbatim
00133 *>          IWORK is INTEGER array, dimension (N).
00134 *>     Workspace.2
00135 *> \endverbatim
00136 *
00137 *  Authors:
00138 *  ========
00139 *
00140 *> \author Univ. of Tennessee 
00141 *> \author Univ. of California Berkeley 
00142 *> \author Univ. of Colorado Denver 
00143 *> \author NAG Ltd. 
00144 *
00145 *> \date November 2011
00146 *
00147 *> \ingroup realGEcomputational
00148 *
00149 *  =====================================================================
00150       REAL FUNCTION SLA_GERCOND ( TRANS, N, A, LDA, AF, LDAF, IPIV,
00151      $                            CMODE, C, INFO, WORK, IWORK )
00152 *
00153 *  -- LAPACK computational routine (version 3.4.0) --
00154 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00155 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00156 *     November 2011
00157 *
00158 *     .. Scalar Arguments ..
00159       CHARACTER          TRANS
00160       INTEGER            N, LDA, LDAF, INFO, CMODE
00161 *     ..
00162 *     .. Array Arguments ..
00163       INTEGER            IPIV( * ), IWORK( * )
00164       REAL               A( LDA, * ), AF( LDAF, * ), WORK( * ),
00165      $                   C( * )
00166 *    ..
00167 *
00168 *  =====================================================================
00169 *
00170 *     .. Local Scalars ..
00171       LOGICAL            NOTRANS
00172       INTEGER            KASE, I, J
00173       REAL               AINVNM, TMP
00174 *     ..
00175 *     .. Local Arrays ..
00176       INTEGER            ISAVE( 3 )
00177 *     ..
00178 *     .. External Functions ..
00179       LOGICAL            LSAME
00180       EXTERNAL           LSAME
00181 *     ..
00182 *     .. External Subroutines ..
00183       EXTERNAL           SLACN2, SGETRS, XERBLA
00184 *     ..
00185 *     .. Intrinsic Functions ..
00186       INTRINSIC          ABS, MAX
00187 *     ..
00188 *     .. Executable Statements ..
00189 *
00190       SLA_GERCOND = 0.0
00191 *
00192       INFO = 0
00193       NOTRANS = LSAME( TRANS, 'N' )
00194       IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T')
00195      $     .AND. .NOT. LSAME(TRANS, 'C') ) THEN
00196          INFO = -1
00197       ELSE IF( N.LT.0 ) THEN
00198          INFO = -2
00199       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00200          INFO = -4
00201       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
00202          INFO = -6
00203       END IF
00204       IF( INFO.NE.0 ) THEN
00205          CALL XERBLA( 'SLA_GERCOND', -INFO )
00206          RETURN
00207       END IF
00208       IF( N.EQ.0 ) THEN
00209          SLA_GERCOND = 1.0
00210          RETURN
00211       END IF
00212 *
00213 *     Compute the equilibration matrix R such that
00214 *     inv(R)*A*C has unit 1-norm.
00215 *
00216       IF (NOTRANS) THEN
00217          DO I = 1, N
00218             TMP = 0.0
00219             IF ( CMODE .EQ. 1 ) THEN
00220                DO J = 1, N
00221                   TMP = TMP + ABS( A( I, J ) * C( J ) )
00222                END DO
00223             ELSE IF ( CMODE .EQ. 0 ) THEN
00224                DO J = 1, N
00225                   TMP = TMP + ABS( A( I, J ) )
00226                END DO
00227             ELSE
00228                DO J = 1, N
00229                   TMP = TMP + ABS( A( I, J ) / C( J ) )
00230                END DO
00231             END IF
00232             WORK( 2*N+I ) = TMP
00233          END DO
00234       ELSE
00235          DO I = 1, N
00236             TMP = 0.0
00237             IF ( CMODE .EQ. 1 ) THEN
00238                DO J = 1, N
00239                   TMP = TMP + ABS( A( J, I ) * C( J ) )
00240                END DO
00241             ELSE IF ( CMODE .EQ. 0 ) THEN
00242                DO J = 1, N
00243                   TMP = TMP + ABS( A( J, I ) )
00244                END DO
00245             ELSE
00246                DO J = 1, N
00247                   TMP = TMP + ABS( A( J, I ) / C( J ) )
00248                END DO
00249             END IF
00250             WORK( 2*N+I ) = TMP
00251          END DO
00252       END IF
00253 *
00254 *     Estimate the norm of inv(op(A)).
00255 *
00256       AINVNM = 0.0
00257 
00258       KASE = 0
00259    10 CONTINUE
00260       CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
00261       IF( KASE.NE.0 ) THEN
00262          IF( KASE.EQ.2 ) THEN
00263 *
00264 *           Multiply by R.
00265 *
00266             DO I = 1, N
00267                WORK(I) = WORK(I) * WORK(2*N+I)
00268             END DO
00269 
00270             IF (NOTRANS) THEN
00271                CALL SGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
00272      $            WORK, N, INFO )
00273             ELSE
00274                CALL SGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
00275      $            WORK, N, INFO )
00276             END IF
00277 *
00278 *           Multiply by inv(C).
00279 *
00280             IF ( CMODE .EQ. 1 ) THEN
00281                DO I = 1, N
00282                   WORK( I ) = WORK( I ) / C( I )
00283                END DO
00284             ELSE IF ( CMODE .EQ. -1 ) THEN
00285                DO I = 1, N
00286                   WORK( I ) = WORK( I ) * C( I )
00287                END DO
00288             END IF
00289          ELSE
00290 *
00291 *           Multiply by inv(C**T).
00292 *
00293             IF ( CMODE .EQ. 1 ) THEN
00294                DO I = 1, N
00295                   WORK( I ) = WORK( I ) / C( I )
00296                END DO
00297             ELSE IF ( CMODE .EQ. -1 ) THEN
00298                DO I = 1, N
00299                   WORK( I ) = WORK( I ) * C( I )
00300                END DO
00301             END IF
00302 
00303             IF (NOTRANS) THEN
00304                CALL SGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
00305      $            WORK, N, INFO )
00306             ELSE
00307                CALL SGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
00308      $            WORK, N, INFO )
00309             END IF
00310 *
00311 *           Multiply by R.
00312 *
00313             DO I = 1, N
00314                WORK( I ) = WORK( I ) * WORK( 2*N+I )
00315             END DO
00316          END IF
00317          GO TO 10
00318       END IF
00319 *
00320 *     Compute the estimate of the reciprocal condition number.
00321 *
00322       IF( AINVNM .NE. 0.0 )
00323      $   SLA_GERCOND = ( 1.0 / AINVNM )
00324 *
00325       RETURN
00326 *
00327       END
 All Files Functions