LAPACK: SRC/dla_gercond.f Source File

Go to the documentation of this file.
00001       DOUBLE PRECISION FUNCTION DLA_GERCOND ( TRANS, N, A, LDA, AF,
00002      $                                        LDAF, IPIV, CMODE, C,
00003      $                                        INFO, WORK, IWORK )
00004 *
00005 *     -- LAPACK routine (version 3.2.1)                                 --
00006 *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
00007 *     -- Jason Riedy of Univ. of California Berkeley.                 --
00008 *     -- April 2009                                                   --
00009 *
00010 *     -- LAPACK is a software package provided by Univ. of Tennessee, --
00011 *     -- Univ. of California Berkeley and NAG Ltd.                    --
00012 *
00013       IMPLICIT NONE
00014 *     ..
00015 *     .. Scalar Arguments ..
00016       CHARACTER          TRANS
00017       INTEGER            N, LDA, LDAF, INFO, CMODE
00018 *     ..
00019 *     .. Array Arguments ..
00020       INTEGER            IPIV( * ), IWORK( * )
00021       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * ),
00022      $                   C( * )
00023 *     ..
00024 *
00025 *  Purpose
00026 *  =======
00027 *
00028 *     DLA_GERCOND estimates the Skeel condition number of op(A) * op2(C)
00029 *     where op2 is determined by CMODE as follows
00030 *     CMODE =  1    op2(C) = C
00031 *     CMODE =  0    op2(C) = I
00032 *     CMODE = -1    op2(C) = inv(C)
00033 *     The Skeel condition number cond(A) = norminf( |inv(A)||A| )
00034 *     is computed by computing scaling factors R such that
00035 *     diag(R)*A*op2(C) is row equilibrated and computing the standard
00036 *     infinity-norm condition number.
00037 *
00038 *  Arguments
00039 *  ==========
00040 *
00041 *     TRANS   (input) CHARACTER*1
00042 *     Specifies the form of the system of equations:
00043 *       = 'N':  A * X = B     (No transpose)
00044 *       = 'T':  A**T * X = B  (Transpose)
00045 *       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose)
00046 *
00047 *     N       (input) INTEGER
00048 *     The number of linear equations, i.e., the order of the
00049 *     matrix A.  N >= 0.
00050 *
00051 *     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
00052 *     On entry, the N-by-N matrix A.
00053 *
00054 *     LDA     (input) INTEGER
00055 *     The leading dimension of the array A.  LDA >= max(1,N).
00056 *
00057 *     AF      (input) DOUBLE PRECISION array, dimension (LDAF,N)
00058 *     The factors L and U from the factorization
00059 *     A = P*L*U as computed by DGETRF.
00060 *
00061 *     LDAF    (input) INTEGER
00062 *     The leading dimension of the array AF.  LDAF >= max(1,N).
00063 *
00064 *     IPIV    (input) INTEGER array, dimension (N)
00065 *     The pivot indices from the factorization A = P*L*U
00066 *     as computed by DGETRF; row i of the matrix was interchanged
00067 *     with row IPIV(i).
00068 *
00069 *     CMODE   (input) INTEGER
00070 *     Determines op2(C) in the formula op(A) * op2(C) as follows:
00071 *     CMODE =  1    op2(C) = C
00072 *     CMODE =  0    op2(C) = I
00073 *     CMODE = -1    op2(C) = inv(C)
00074 *
00075 *     C       (input) DOUBLE PRECISION array, dimension (N)
00076 *     The vector C in the formula op(A) * op2(C).
00077 *
00078 *     INFO    (output) INTEGER
00079 *       = 0:  Successful exit.
00080 *     i > 0:  The ith argument is invalid.
00081 *
00082 *     WORK    (input) DOUBLE PRECISION array, dimension (3*N).
00083 *     Workspace.
00084 *
00085 *     IWORK   (input) INTEGER array, dimension (N).
00086 *     Workspace.
00087 *
00088 *  =====================================================================
00089 *
00090 *     .. Local Scalars ..
00091       LOGICAL            NOTRANS
00092       INTEGER            KASE, I, J
00093       DOUBLE PRECISION   AINVNM, TMP
00094 *     ..
00095 *     .. Local Arrays ..
00096       INTEGER            ISAVE( 3 )
00097 *     ..
00098 *     .. External Functions ..
00099       LOGICAL            LSAME
00100       EXTERNAL           LSAME
00101 *     ..
00102 *     .. External Subroutines ..
00103       EXTERNAL           DLACN2, DGETRS, XERBLA
00104 *     ..
00105 *     .. Intrinsic Functions ..
00106       INTRINSIC          ABS, MAX
00107 *     ..
00108 *     .. Executable Statements ..
00109 *
00110       DLA_GERCOND = 0.0D+0
00111 *
00112       INFO = 0
00113       NOTRANS = LSAME( TRANS, 'N' )
00114       IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T')
00115      $     .AND. .NOT. LSAME(TRANS, 'C') ) THEN
00116          INFO = -1
00117       ELSE IF( N.LT.0 ) THEN
00118          INFO = -2
00119       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00120          INFO = -4
00121       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
00122          INFO = -6
00123       END IF
00124       IF( INFO.NE.0 ) THEN
00125          CALL XERBLA( 'DLA_GERCOND', -INFO )
00126          RETURN
00127       END IF
00128       IF( N.EQ.0 ) THEN
00129          DLA_GERCOND = 1.0D+0
00130          RETURN
00131       END IF
00132 *
00133 *     Compute the equilibration matrix R such that
00134 *     inv(R)*A*C has unit 1-norm.
00135 *
00136       IF (NOTRANS) THEN
00137          DO I = 1, N
00138             TMP = 0.0D+0
00139             IF ( CMODE .EQ. 1 ) THEN
00140                DO J = 1, N
00141                   TMP = TMP + ABS( A( I, J ) * C( J ) )
00142                END DO
00143             ELSE IF ( CMODE .EQ. 0 ) THEN
00144                DO J = 1, N
00145                   TMP = TMP + ABS( A( I, J ) )
00146                END DO
00147             ELSE
00148                DO J = 1, N
00149                   TMP = TMP + ABS( A( I, J ) / C( J ) )
00150                END DO
00151             END IF
00152             WORK( 2*N+I ) = TMP
00153          END DO
00154       ELSE
00155          DO I = 1, N
00156             TMP = 0.0D+0
00157             IF ( CMODE .EQ. 1 ) THEN
00158                DO J = 1, N
00159                   TMP = TMP + ABS( A( J, I ) * C( J ) )
00160                END DO
00161             ELSE IF ( CMODE .EQ. 0 ) THEN
00162                DO J = 1, N
00163                   TMP = TMP + ABS( A( J, I ) )
00164                END DO
00165             ELSE
00166                DO J = 1, N
00167                   TMP = TMP + ABS( A( J, I ) / C( J ) )
00168                END DO
00169             END IF
00170             WORK( 2*N+I ) = TMP
00171          END DO
00172       END IF
00173 *
00174 *     Estimate the norm of inv(op(A)).
00175 *
00176       AINVNM = 0.0D+0
00177 
00178       KASE = 0
00179    10 CONTINUE
00180       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
00181       IF( KASE.NE.0 ) THEN
00182          IF( KASE.EQ.2 ) THEN
00183 *
00184 *           Multiply by R.
00185 *
00186             DO I = 1, N
00187                WORK(I) = WORK(I) * WORK(2*N+I)
00188             END DO
00189 
00190             IF (NOTRANS) THEN
00191                CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
00192      $            WORK, N, INFO )
00193             ELSE
00194                CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
00195      $            WORK, N, INFO )
00196             END IF
00197 *
00198 *           Multiply by inv(C).
00199 *
00200             IF ( CMODE .EQ. 1 ) THEN
00201                DO I = 1, N
00202                   WORK( I ) = WORK( I ) / C( I )
00203                END DO
00204             ELSE IF ( CMODE .EQ. -1 ) THEN
00205                DO I = 1, N
00206                   WORK( I ) = WORK( I ) * C( I )
00207                END DO
00208             END IF
00209          ELSE
00210 *
00211 *           Multiply by inv(C').
00212 *
00213             IF ( CMODE .EQ. 1 ) THEN
00214                DO I = 1, N
00215                   WORK( I ) = WORK( I ) / C( I )
00216                END DO
00217             ELSE IF ( CMODE .EQ. -1 ) THEN
00218                DO I = 1, N
00219                   WORK( I ) = WORK( I ) * C( I )
00220                END DO
00221             END IF
00222 
00223             IF (NOTRANS) THEN
00224                CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
00225      $            WORK, N, INFO )
00226             ELSE
00227                CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
00228      $            WORK, N, INFO )
00229             END IF
00230 *
00231 *           Multiply by R.
00232 *
00233             DO I = 1, N
00234                WORK( I ) = WORK( I ) * WORK( 2*N+I )
00235             END DO
00236          END IF
00237          GO TO 10
00238       END IF
00239 *
00240 *     Compute the estimate of the reciprocal condition number.
00241 *
00242       IF( AINVNM .NE. 0.0D+0 )
00243      $   DLA_GERCOND = ( 1.0D+0 / AINVNM )
00244 *
00245       RETURN
00246 *
00247       END