Mercurial > octave-nkf
diff libcruft/lapack/cpocon.f @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author | David Bateman <dbateman@free.fr> |
---|---|
date | Sun, 27 Apr 2008 22:34:17 +0200 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/cpocon.f Sun Apr 27 22:34:17 2008 +0200 @@ -0,0 +1,184 @@ + SUBROUTINE CPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, + $ INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH. +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N + REAL ANORM, RCOND +* .. +* .. Array Arguments .. + REAL RWORK( * ) + COMPLEX A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* CPOCON estimates the reciprocal of the condition number (in the +* 1-norm) of a complex Hermitian positive definite matrix using the +* Cholesky factorization A = U**H*U or A = L*L**H computed by CPOTRF. +* +* An estimate is obtained for norm(inv(A)), and the reciprocal of the +* condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* = 'U': Upper triangle of A is stored; +* = 'L': Lower triangle of A is stored. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input) COMPLEX array, dimension (LDA,N) +* The triangular factor U or L from the Cholesky factorization +* A = U**H*U or A = L*L**H, as computed by CPOTRF. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* ANORM (input) REAL +* The 1-norm (or infinity-norm) of the Hermitian matrix A. +* +* RCOND (output) REAL +* The reciprocal of the condition number of the matrix A, +* computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an +* estimate of the 1-norm of inv(A) computed in this routine. +* +* WORK (workspace) COMPLEX array, dimension (2*N) +* +* RWORK (workspace) REAL array, dimension (N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + CHARACTER NORMIN + INTEGER IX, KASE + REAL AINVNM, SCALE, SCALEL, SCALEU, SMLNUM + COMPLEX ZDUM +* .. +* .. Local Arrays .. + INTEGER ISAVE( 3 ) +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ICAMAX + REAL SLAMCH + EXTERNAL LSAME, ICAMAX, SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL CLACN2, CLATRS, CSRSCL, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, AIMAG, MAX, REAL +* .. +* .. Statement Functions .. + REAL CABS1 +* .. +* .. Statement Function definitions .. + CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( ANORM.LT.ZERO ) THEN + INFO = -5 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CPOCON', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + RCOND = ZERO + IF( N.EQ.0 ) THEN + RCOND = ONE + RETURN + ELSE IF( ANORM.EQ.ZERO ) THEN + RETURN + END IF +* + SMLNUM = SLAMCH( 'Safe minimum' ) +* +* Estimate the 1-norm of inv(A). +* + KASE = 0 + NORMIN = 'N' + 10 CONTINUE + CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) + IF( KASE.NE.0 ) THEN + IF( UPPER ) THEN +* +* Multiply by inv(U'). +* + CALL CLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', + $ NORMIN, N, A, LDA, WORK, SCALEL, RWORK, INFO ) + NORMIN = 'Y' +* +* Multiply by inv(U). +* + CALL CLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, + $ A, LDA, WORK, SCALEU, RWORK, INFO ) + ELSE +* +* Multiply by inv(L). +* + CALL CLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, + $ A, LDA, WORK, SCALEL, RWORK, INFO ) + NORMIN = 'Y' +* +* Multiply by inv(L'). +* + CALL CLATRS( 'Lower', 'Conjugate transpose', 'Non-unit', + $ NORMIN, N, A, LDA, WORK, SCALEU, RWORK, INFO ) + END IF +* +* Multiply by 1/SCALE if doing so will not cause overflow. +* + SCALE = SCALEL*SCALEU + IF( SCALE.NE.ONE ) THEN + IX = ICAMAX( N, WORK, 1 ) + IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) + $ GO TO 20 + CALL CSRSCL( N, SCALE, WORK, 1 ) + END IF + GO TO 10 + END IF +* +* Compute the estimate of the reciprocal condition number. +* + IF( AINVNM.NE.ZERO ) + $ RCOND = ( ONE / AINVNM ) / ANORM +* + 20 CONTINUE + RETURN +* +* End of CPOCON +* + END