diff options
author | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
---|---|---|
committer | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
commit | baba851215b44ac3b60b9248eb02bcce7eb76247 (patch) | |
tree | 8c0f5c006875532a30d4409f5e94b0f310ff00a7 /TESTING/LIN/cppt01.f | |
download | lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.gz lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.bz2 lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.zip |
Move LAPACK trunk into position.
Diffstat (limited to 'TESTING/LIN/cppt01.f')
-rw-r--r-- | TESTING/LIN/cppt01.f | 192 |
1 files changed, 192 insertions, 0 deletions
diff --git a/TESTING/LIN/cppt01.f b/TESTING/LIN/cppt01.f new file mode 100644 index 00000000..26f432be --- /dev/null +++ b/TESTING/LIN/cppt01.f @@ -0,0 +1,192 @@ + SUBROUTINE CPPT01( UPLO, N, A, AFAC, RWORK, RESID ) +* +* -- LAPACK test routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER N + REAL RESID +* .. +* .. Array Arguments .. + REAL RWORK( * ) + COMPLEX A( * ), AFAC( * ) +* .. +* +* Purpose +* ======= +* +* CPPT01 reconstructs a Hermitian positive definite packed matrix A +* from its L*L' or U'*U factorization and computes the residual +* norm( L*L' - A ) / ( N * norm(A) * EPS ) or +* norm( U'*U - A ) / ( N * norm(A) * EPS ), +* where EPS is the machine epsilon, L' is the conjugate transpose of +* L, and U' is the conjugate transpose of U. +* +* Arguments +* ========== +* +* UPLO (input) CHARACTER*1 +* Specifies whether the upper or lower triangular part of the +* Hermitian matrix A is stored: +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* N (input) INTEGER +* The number of rows and columns of the matrix A. N >= 0. +* +* A (input) COMPLEX array, dimension (N*(N+1)/2) +* The original Hermitian matrix A, stored as a packed +* triangular matrix. +* +* AFAC (input/output) COMPLEX array, dimension (N*(N+1)/2) +* On entry, the factor L or U from the L*L' or U'*U +* factorization of A, stored as a packed triangular matrix. +* Overwritten with the reconstructed matrix, and then with the +* difference L*L' - A (or U'*U - A). +* +* RWORK (workspace) REAL array, dimension (N) +* +* RESID (output) REAL +* If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS ) +* If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS ) +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER I, K, KC + REAL ANORM, EPS, TR + COMPLEX TC +* .. +* .. External Functions .. + LOGICAL LSAME + REAL CLANHP, SLAMCH + COMPLEX CDOTC + EXTERNAL LSAME, CLANHP, SLAMCH, CDOTC +* .. +* .. External Subroutines .. + EXTERNAL CHPR, CSCAL, CTPMV +* .. +* .. Intrinsic Functions .. + INTRINSIC AIMAG, REAL +* .. +* .. Executable Statements .. +* +* Quick exit if N = 0 +* + IF( N.LE.0 ) THEN + RESID = ZERO + RETURN + END IF +* +* Exit with RESID = 1/EPS if ANORM = 0. +* + EPS = SLAMCH( 'Epsilon' ) + ANORM = CLANHP( '1', UPLO, N, A, RWORK ) + IF( ANORM.LE.ZERO ) THEN + RESID = ONE / EPS + RETURN + END IF +* +* Check the imaginary parts of the diagonal elements and return with +* an error code if any are nonzero. +* + KC = 1 + IF( LSAME( UPLO, 'U' ) ) THEN + DO 10 K = 1, N + IF( AIMAG( AFAC( KC ) ).NE.ZERO ) THEN + RESID = ONE / EPS + RETURN + END IF + KC = KC + K + 1 + 10 CONTINUE + ELSE + DO 20 K = 1, N + IF( AIMAG( AFAC( KC ) ).NE.ZERO ) THEN + RESID = ONE / EPS + RETURN + END IF + KC = KC + N - K + 1 + 20 CONTINUE + END IF +* +* Compute the product U'*U, overwriting U. +* + IF( LSAME( UPLO, 'U' ) ) THEN + KC = ( N*( N-1 ) ) / 2 + 1 + DO 30 K = N, 1, -1 +* +* Compute the (K,K) element of the result. +* + TR = CDOTC( K, AFAC( KC ), 1, AFAC( KC ), 1 ) + AFAC( KC+K-1 ) = TR +* +* Compute the rest of column K. +* + IF( K.GT.1 ) THEN + CALL CTPMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC, + $ AFAC( KC ), 1 ) + KC = KC - ( K-1 ) + END IF + 30 CONTINUE +* +* Compute the difference L*L' - A +* + KC = 1 + DO 50 K = 1, N + DO 40 I = 1, K - 1 + AFAC( KC+I-1 ) = AFAC( KC+I-1 ) - A( KC+I-1 ) + 40 CONTINUE + AFAC( KC+K-1 ) = AFAC( KC+K-1 ) - REAL( A( KC+K-1 ) ) + KC = KC + K + 50 CONTINUE +* +* Compute the product L*L', overwriting L. +* + ELSE + KC = ( N*( N+1 ) ) / 2 + DO 60 K = N, 1, -1 +* +* Add a multiple of column K of the factor L to each of +* columns K+1 through N. +* + IF( K.LT.N ) + $ CALL CHPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1, + $ AFAC( KC+N-K+1 ) ) +* +* Scale column K by the diagonal element. +* + TC = AFAC( KC ) + CALL CSCAL( N-K+1, TC, AFAC( KC ), 1 ) +* + KC = KC - ( N-K+2 ) + 60 CONTINUE +* +* Compute the difference U'*U - A +* + KC = 1 + DO 80 K = 1, N + AFAC( KC ) = AFAC( KC ) - REAL( A( KC ) ) + DO 70 I = K + 1, N + AFAC( KC+I-K ) = AFAC( KC+I-K ) - A( KC+I-K ) + 70 CONTINUE + KC = KC + N - K + 1 + 80 CONTINUE + END IF +* +* Compute norm( L*U - A ) / ( N * norm(A) * EPS ) +* + RESID = CLANHP( '1', UPLO, N, AFAC, RWORK ) +* + RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS +* + RETURN +* +* End of CPPT01 +* + END |