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