diff options
author | julie <julielangou@users.noreply.github.com> | 2008-12-16 17:06:58 +0000 |
---|---|---|
committer | julie <julielangou@users.noreply.github.com> | 2008-12-16 17:06:58 +0000 |
commit | ff981f106bde4ce6a74aa4f4a572c943f5a395b2 (patch) | |
tree | a386cad907bcaefd6893535c31d67ec9468e693e /SRC/dla_gercond.f | |
parent | e58b61578b55644f6391f3333262b72c1dc88437 (diff) | |
download | lapack-ff981f106bde4ce6a74aa4f4a572c943f5a395b2.tar.gz lapack-ff981f106bde4ce6a74aa4f4a572c943f5a395b2.tar.bz2 lapack-ff981f106bde4ce6a74aa4f4a572c943f5a395b2.zip |
Diffstat (limited to 'SRC/dla_gercond.f')
-rw-r--r-- | SRC/dla_gercond.f | 189 |
1 files changed, 189 insertions, 0 deletions
diff --git a/SRC/dla_gercond.f b/SRC/dla_gercond.f new file mode 100644 index 00000000..cb75a97e --- /dev/null +++ b/SRC/dla_gercond.f @@ -0,0 +1,189 @@ + DOUBLE PRECISION FUNCTION DLA_GERCOND ( TRANS, N, A, LDA, AF, + $ LDAF, IPIV, CMODE, C, INFO, WORK, + $ IWORK ) +* +* -- LAPACK routine (version 3.2) -- +* -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- +* -- Jason Riedy of Univ. of California Berkeley. -- +* -- November 2008 -- +* +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley and NAG Ltd. -- +* + IMPLICIT NONE +* .. +* .. Scalar Arguments .. + CHARACTER TRANS + INTEGER N, LDA, LDAF, INFO, CMODE +* .. +* .. Array Arguments .. + INTEGER IPIV( * ), IWORK( * ) + DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ), + $ C( * ) +* +* DLA_GERCOND estimates the Skeel condition number of op(A) * op2(C) +* where op2 is determined by CMODE as follows +* CMODE = 1 op2(C) = C +* CMODE = 0 op2(C) = I +* CMODE = -1 op2(C) = inv(C) +* The Skeel condition number cond(A) = norminf( |inv(A)||A| ) +* is computed by computing scaling factors R such that +* diag(R)*A*op2(C) is row equilibrated and computing the standard +* infinity-norm condition number. +* WORK is a DOUBLE PRECISION workspace of size 3*N, and +* IWORK is an INTEGER workspace of size N. +* .. +* .. Local Scalars .. + LOGICAL NOTRANS + INTEGER KASE, I, J + DOUBLE PRECISION AINVNM, TMP +* .. +* .. Local Arrays .. + INTEGER ISAVE( 3 ) +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL DLACN2, DGETRS, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX +* .. +* .. Executable Statements .. +* + DLA_GERCOND = 0.0D+0 +* + INFO = 0 + NOTRANS = LSAME( TRANS, 'N' ) + IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T') + $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLA_GERCOND', -INFO ) + RETURN + END IF + IF( N.EQ.0 ) THEN + DLA_GERCOND = 1.0D+0 + RETURN + END IF +* +* Compute the equilibration matrix R such that +* inv(R)*A*C has unit 1-norm. +* + IF (NOTRANS) THEN + DO I = 1, N + TMP = 0.0D+0 + IF ( CMODE .EQ. 1 ) THEN + DO J = 1, N + TMP = TMP + ABS( A( I, J ) * C( J ) ) + END DO + ELSE IF ( CMODE .EQ. 0 ) THEN + DO J = 1, N + TMP = TMP + ABS( A( I, J ) ) + END DO + ELSE + DO J = 1, N + TMP = TMP + ABS( A( I, J ) / C( J ) ) + END DO + END IF + WORK( 2*N+I ) = TMP + END DO + ELSE + DO I = 1, N + TMP = 0.0D+0 + IF ( CMODE .EQ. 1 ) THEN + DO J = 1, N + TMP = TMP + ABS( A( J, I ) * C( J ) ) + END DO + ELSE IF ( CMODE .EQ. 0 ) THEN + DO J = 1, N + TMP = TMP + ABS( A( J, I ) ) + END DO + ELSE + DO J = 1, N + TMP = TMP + ABS( A( J, I ) / C( J ) ) + END DO + END IF + WORK( 2*N+I ) = TMP + END DO + END IF +* +* Estimate the norm of inv(op(A)). +* + AINVNM = 0.0D+0 + + KASE = 0 + 10 CONTINUE + CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) + IF( KASE.NE.0 ) THEN + IF( KASE.EQ.2 ) THEN +* +* Multiply by R. +* + DO I = 1, N + WORK(I) = WORK(I) * WORK(2*N+I) + END DO + + IF (NOTRANS) THEN + CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, + $ WORK, N, INFO ) + ELSE + CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV, + $ WORK, N, INFO ) + END IF +* +* Multiply by inv(C). +* + IF ( CMODE .EQ. 1 ) THEN + DO I = 1, N + WORK( I ) = WORK( I ) / C( I ) + END DO + ELSE IF ( CMODE .EQ. -1 ) THEN + DO I = 1, N + WORK( I ) = WORK( I ) * C( I ) + END DO + END IF + ELSE +* +* Multiply by inv(C'). +* + IF ( CMODE .EQ. 1 ) THEN + DO I = 1, N + WORK( I ) = WORK( I ) / C( I ) + END DO + ELSE IF ( CMODE .EQ. -1 ) THEN + DO I = 1, N + WORK( I ) = WORK( I ) * C( I ) + END DO + END IF + + IF (NOTRANS) THEN + CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV, + $ WORK, N, INFO ) + ELSE + CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, + $ WORK, N, INFO ) + END IF +* +* Multiply by R. +* + DO I = 1, N + WORK( I ) = WORK( I ) * WORK( 2*N+I ) + END DO + END IF + GO TO 10 + END IF +* +* Compute the estimate of the reciprocal condition number. +* + IF( AINVNM .NE. 0.0D+0 ) + $ DLA_GERCOND = ( 1.0D+0 / AINVNM ) +* + RETURN +* + END |