diff options
Diffstat (limited to 'TESTING/LIN/chpt01.f')
-rw-r--r-- | TESTING/LIN/chpt01.f | 175 |
1 files changed, 175 insertions, 0 deletions
diff --git a/TESTING/LIN/chpt01.f b/TESTING/LIN/chpt01.f new file mode 100644 index 00000000..2710562a --- /dev/null +++ b/TESTING/LIN/chpt01.f @@ -0,0 +1,175 @@ + SUBROUTINE CHPT01( UPLO, N, A, AFAC, IPIV, C, LDC, 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 LDC, N + REAL RESID +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + REAL RWORK( * ) + COMPLEX A( * ), AFAC( * ), C( LDC, * ) +* .. +* +* Purpose +* ======= +* +* CHPT01 reconstructs a Hermitian indefinite packed matrix A from its +* block L*D*L' or U*D*U' factorization and computes the residual +* norm( C - A ) / ( N * norm(A) * EPS ), +* where C is the reconstructed matrix, 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) COMPLEX array, dimension (N*(N+1)/2) +* The factored form of the matrix A, stored as a packed +* triangular matrix. AFAC contains the block diagonal matrix D +* and the multipliers used to obtain the factor L or U from the +* block L*D*L' or U*D*U' factorization as computed by CHPTRF. +* +* IPIV (input) INTEGER array, dimension (N) +* The pivot indices from CHPTRF. +* +* C (workspace) COMPLEX array, dimension (LDC,N) +* +* LDC (integer) INTEGER +* The leading dimension of the array C. LDC >= max(1,N). +* +* RWORK (workspace) REAL array, dimension (N) +* +* RESID (output) REAL +* If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) +* If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + COMPLEX CZERO, CONE + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), + $ CONE = ( 1.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, INFO, J, JC + REAL ANORM, EPS +* .. +* .. External Functions .. + LOGICAL LSAME + REAL CLANHE, CLANHP, SLAMCH + EXTERNAL LSAME, CLANHE, CLANHP, SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL CLAVHP, CLASET +* .. +* .. Intrinsic Functions .. + INTRINSIC AIMAG, REAL +* .. +* .. Executable Statements .. +* +* Quick exit if N = 0. +* + IF( N.LE.0 ) THEN + RESID = ZERO + RETURN + END IF +* +* Determine EPS and the norm of A. +* + EPS = SLAMCH( 'Epsilon' ) + ANORM = CLANHP( '1', UPLO, N, A, RWORK ) +* +* Check the imaginary parts of the diagonal elements and return with +* an error code if any are nonzero. +* + JC = 1 + IF( LSAME( UPLO, 'U' ) ) THEN + DO 10 J = 1, N + IF( AIMAG( AFAC( JC ) ).NE.ZERO ) THEN + RESID = ONE / EPS + RETURN + END IF + JC = JC + J + 1 + 10 CONTINUE + ELSE + DO 20 J = 1, N + IF( AIMAG( AFAC( JC ) ).NE.ZERO ) THEN + RESID = ONE / EPS + RETURN + END IF + JC = JC + N - J + 1 + 20 CONTINUE + END IF +* +* Initialize C to the identity matrix. +* + CALL CLASET( 'Full', N, N, CZERO, CONE, C, LDC ) +* +* Call CLAVHP to form the product D * U' (or D * L' ). +* + CALL CLAVHP( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC, IPIV, C, + $ LDC, INFO ) +* +* Call CLAVHP again to multiply by U ( or L ). +* + CALL CLAVHP( UPLO, 'No transpose', 'Unit', N, N, AFAC, IPIV, C, + $ LDC, INFO ) +* +* Compute the difference C - A . +* + IF( LSAME( UPLO, 'U' ) ) THEN + JC = 0 + DO 40 J = 1, N + DO 30 I = 1, J - 1 + C( I, J ) = C( I, J ) - A( JC+I ) + 30 CONTINUE + C( J, J ) = C( J, J ) - REAL( A( JC+J ) ) + JC = JC + J + 40 CONTINUE + ELSE + JC = 1 + DO 60 J = 1, N + C( J, J ) = C( J, J ) - REAL( A( JC ) ) + DO 50 I = J + 1, N + C( I, J ) = C( I, J ) - A( JC+I-J ) + 50 CONTINUE + JC = JC + N - J + 1 + 60 CONTINUE + END IF +* +* Compute norm( C - A ) / ( N * norm(A) * EPS ) +* + RESID = CLANHE( '1', UPLO, N, C, LDC, RWORK ) +* + IF( ANORM.LE.ZERO ) THEN + IF( RESID.NE.ZERO ) + $ RESID = ONE / EPS + ELSE + RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS + END IF +* + RETURN +* +* End of CHPT01 +* + END |