Annotation of rpl/lapack/lapack/dla_gercond.f, revision 1.1
1.1 ! bertrand 1: DOUBLE PRECISION FUNCTION DLA_GERCOND ( TRANS, N, A, LDA, AF,
! 2: $ LDAF, IPIV, CMODE, C,
! 3: $ INFO, WORK, IWORK )
! 4: *
! 5: * -- LAPACK routine (version 3.2.1) --
! 6: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
! 7: * -- Jason Riedy of Univ. of California Berkeley. --
! 8: * -- April 2009 --
! 9: *
! 10: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 11: * -- Univ. of California Berkeley and NAG Ltd. --
! 12: *
! 13: IMPLICIT NONE
! 14: * ..
! 15: * .. Scalar Arguments ..
! 16: CHARACTER TRANS
! 17: INTEGER N, LDA, LDAF, INFO, CMODE
! 18: * ..
! 19: * .. Array Arguments ..
! 20: INTEGER IPIV( * ), IWORK( * )
! 21: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ),
! 22: $ C( * )
! 23: * ..
! 24: *
! 25: * Purpose
! 26: * =======
! 27: *
! 28: * DLA_GERCOND estimates the Skeel condition number of op(A) * op2(C)
! 29: * where op2 is determined by CMODE as follows
! 30: * CMODE = 1 op2(C) = C
! 31: * CMODE = 0 op2(C) = I
! 32: * CMODE = -1 op2(C) = inv(C)
! 33: * The Skeel condition number cond(A) = norminf( |inv(A)||A| )
! 34: * is computed by computing scaling factors R such that
! 35: * diag(R)*A*op2(C) is row equilibrated and computing the standard
! 36: * infinity-norm condition number.
! 37: *
! 38: * Arguments
! 39: * ==========
! 40: *
! 41: * TRANS (input) CHARACTER*1
! 42: * Specifies the form of the system of equations:
! 43: * = 'N': A * X = B (No transpose)
! 44: * = 'T': A**T * X = B (Transpose)
! 45: * = 'C': A**H * X = B (Conjugate Transpose = Transpose)
! 46: *
! 47: * N (input) INTEGER
! 48: * The number of linear equations, i.e., the order of the
! 49: * matrix A. N >= 0.
! 50: *
! 51: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
! 52: * On entry, the N-by-N matrix A.
! 53: *
! 54: * LDA (input) INTEGER
! 55: * The leading dimension of the array A. LDA >= max(1,N).
! 56: *
! 57: * AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
! 58: * The factors L and U from the factorization
! 59: * A = P*L*U as computed by DGETRF.
! 60: *
! 61: * LDAF (input) INTEGER
! 62: * The leading dimension of the array AF. LDAF >= max(1,N).
! 63: *
! 64: * IPIV (input) INTEGER array, dimension (N)
! 65: * The pivot indices from the factorization A = P*L*U
! 66: * as computed by DGETRF; row i of the matrix was interchanged
! 67: * with row IPIV(i).
! 68: *
! 69: * CMODE (input) INTEGER
! 70: * Determines op2(C) in the formula op(A) * op2(C) as follows:
! 71: * CMODE = 1 op2(C) = C
! 72: * CMODE = 0 op2(C) = I
! 73: * CMODE = -1 op2(C) = inv(C)
! 74: *
! 75: * C (input) DOUBLE PRECISION array, dimension (N)
! 76: * The vector C in the formula op(A) * op2(C).
! 77: *
! 78: * INFO (output) INTEGER
! 79: * = 0: Successful exit.
! 80: * i > 0: The ith argument is invalid.
! 81: *
! 82: * WORK (input) DOUBLE PRECISION array, dimension (3*N).
! 83: * Workspace.
! 84: *
! 85: * IWORK (input) INTEGER array, dimension (N).
! 86: * Workspace.
! 87: *
! 88: * =====================================================================
! 89: *
! 90: * .. Local Scalars ..
! 91: LOGICAL NOTRANS
! 92: INTEGER KASE, I, J
! 93: DOUBLE PRECISION AINVNM, TMP
! 94: * ..
! 95: * .. Local Arrays ..
! 96: INTEGER ISAVE( 3 )
! 97: * ..
! 98: * .. External Functions ..
! 99: LOGICAL LSAME
! 100: EXTERNAL LSAME
! 101: * ..
! 102: * .. External Subroutines ..
! 103: EXTERNAL DLACN2, DGETRS, XERBLA
! 104: * ..
! 105: * .. Intrinsic Functions ..
! 106: INTRINSIC ABS, MAX
! 107: * ..
! 108: * .. Executable Statements ..
! 109: *
! 110: DLA_GERCOND = 0.0D+0
! 111: *
! 112: INFO = 0
! 113: NOTRANS = LSAME( TRANS, 'N' )
! 114: IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T')
! 115: $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN
! 116: INFO = -1
! 117: ELSE IF( N.LT.0 ) THEN
! 118: INFO = -2
! 119: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 120: INFO = -4
! 121: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
! 122: INFO = -6
! 123: END IF
! 124: IF( INFO.NE.0 ) THEN
! 125: CALL XERBLA( 'DLA_GERCOND', -INFO )
! 126: RETURN
! 127: END IF
! 128: IF( N.EQ.0 ) THEN
! 129: DLA_GERCOND = 1.0D+0
! 130: RETURN
! 131: END IF
! 132: *
! 133: * Compute the equilibration matrix R such that
! 134: * inv(R)*A*C has unit 1-norm.
! 135: *
! 136: IF (NOTRANS) THEN
! 137: DO I = 1, N
! 138: TMP = 0.0D+0
! 139: IF ( CMODE .EQ. 1 ) THEN
! 140: DO J = 1, N
! 141: TMP = TMP + ABS( A( I, J ) * C( J ) )
! 142: END DO
! 143: ELSE IF ( CMODE .EQ. 0 ) THEN
! 144: DO J = 1, N
! 145: TMP = TMP + ABS( A( I, J ) )
! 146: END DO
! 147: ELSE
! 148: DO J = 1, N
! 149: TMP = TMP + ABS( A( I, J ) / C( J ) )
! 150: END DO
! 151: END IF
! 152: WORK( 2*N+I ) = TMP
! 153: END DO
! 154: ELSE
! 155: DO I = 1, N
! 156: TMP = 0.0D+0
! 157: IF ( CMODE .EQ. 1 ) THEN
! 158: DO J = 1, N
! 159: TMP = TMP + ABS( A( J, I ) * C( J ) )
! 160: END DO
! 161: ELSE IF ( CMODE .EQ. 0 ) THEN
! 162: DO J = 1, N
! 163: TMP = TMP + ABS( A( J, I ) )
! 164: END DO
! 165: ELSE
! 166: DO J = 1, N
! 167: TMP = TMP + ABS( A( J, I ) / C( J ) )
! 168: END DO
! 169: END IF
! 170: WORK( 2*N+I ) = TMP
! 171: END DO
! 172: END IF
! 173: *
! 174: * Estimate the norm of inv(op(A)).
! 175: *
! 176: AINVNM = 0.0D+0
! 177:
! 178: KASE = 0
! 179: 10 CONTINUE
! 180: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
! 181: IF( KASE.NE.0 ) THEN
! 182: IF( KASE.EQ.2 ) THEN
! 183: *
! 184: * Multiply by R.
! 185: *
! 186: DO I = 1, N
! 187: WORK(I) = WORK(I) * WORK(2*N+I)
! 188: END DO
! 189:
! 190: IF (NOTRANS) THEN
! 191: CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
! 192: $ WORK, N, INFO )
! 193: ELSE
! 194: CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
! 195: $ WORK, N, INFO )
! 196: END IF
! 197: *
! 198: * Multiply by inv(C).
! 199: *
! 200: IF ( CMODE .EQ. 1 ) THEN
! 201: DO I = 1, N
! 202: WORK( I ) = WORK( I ) / C( I )
! 203: END DO
! 204: ELSE IF ( CMODE .EQ. -1 ) THEN
! 205: DO I = 1, N
! 206: WORK( I ) = WORK( I ) * C( I )
! 207: END DO
! 208: END IF
! 209: ELSE
! 210: *
! 211: * Multiply by inv(C').
! 212: *
! 213: IF ( CMODE .EQ. 1 ) THEN
! 214: DO I = 1, N
! 215: WORK( I ) = WORK( I ) / C( I )
! 216: END DO
! 217: ELSE IF ( CMODE .EQ. -1 ) THEN
! 218: DO I = 1, N
! 219: WORK( I ) = WORK( I ) * C( I )
! 220: END DO
! 221: END IF
! 222:
! 223: IF (NOTRANS) THEN
! 224: CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
! 225: $ WORK, N, INFO )
! 226: ELSE
! 227: CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
! 228: $ WORK, N, INFO )
! 229: END IF
! 230: *
! 231: * Multiply by R.
! 232: *
! 233: DO I = 1, N
! 234: WORK( I ) = WORK( I ) * WORK( 2*N+I )
! 235: END DO
! 236: END IF
! 237: GO TO 10
! 238: END IF
! 239: *
! 240: * Compute the estimate of the reciprocal condition number.
! 241: *
! 242: IF( AINVNM .NE. 0.0D+0 )
! 243: $ DLA_GERCOND = ( 1.0D+0 / AINVNM )
! 244: *
! 245: RETURN
! 246: *
! 247: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dla_gercond.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack