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/sgtt01.f |
Move LAPACK trunk into position.
Diffstat (limited to 'TESTING/LIN/sgtt01.f')
-rw-r--r-- | TESTING/LIN/sgtt01.f | 179 |
1 files changed, 179 insertions, 0 deletions
diff --git a/TESTING/LIN/sgtt01.f b/TESTING/LIN/sgtt01.f new file mode 100644 index 00000000..fc124ad0 --- /dev/null +++ b/TESTING/LIN/sgtt01.f @@ -0,0 +1,179 @@ + SUBROUTINE SGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK, + $ LDWORK, RWORK, RESID ) +* +* -- LAPACK test routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER LDWORK, N + REAL RESID +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + REAL D( * ), DF( * ), DL( * ), DLF( * ), DU( * ), + $ DU2( * ), DUF( * ), RWORK( * ), + $ WORK( LDWORK, * ) +* .. +* +* Purpose +* ======= +* +* SGTT01 reconstructs a tridiagonal matrix A from its LU factorization +* and computes the residual +* norm(L*U - A) / ( norm(A) * EPS ), +* where EPS is the machine epsilon. +* +* Arguments +* ========= +* +* N (input) INTEGTER +* The order of the matrix A. N >= 0. +* +* DL (input) REAL array, dimension (N-1) +* The (n-1) sub-diagonal elements of A. +* +* D (input) REAL array, dimension (N) +* The diagonal elements of A. +* +* DU (input) REAL array, dimension (N-1) +* The (n-1) super-diagonal elements of A. +* +* DLF (input) REAL array, dimension (N-1) +* The (n-1) multipliers that define the matrix L from the +* LU factorization of A. +* +* DF (input) REAL array, dimension (N) +* The n diagonal elements of the upper triangular matrix U from +* the LU factorization of A. +* +* DUF (input) REAL array, dimension (N-1) +* The (n-1) elements of the first super-diagonal of U. +* +* DU2F (input) REAL array, dimension (N-2) +* The (n-2) elements of the second super-diagonal of U. +* +* IPIV (input) INTEGER array, dimension (N) +* The pivot indices; for 1 <= i <= n, row i of the matrix was +* interchanged with row IPIV(i). IPIV(i) will always be either +* i or i+1; IPIV(i) = i indicates a row interchange was not +* required. +* +* WORK (workspace) REAL array, dimension (LDWORK,N) +* +* LDWORK (input) INTEGER +* The leading dimension of the array WORK. LDWORK >= max(1,N). +* +* RWORK (workspace) REAL array, dimension (N) +* +* RESID (output) REAL +* The scaled residual: norm(L*U - A) / (norm(A) * EPS) +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER I, IP, J, LASTJ + REAL ANORM, EPS, LI +* .. +* .. External Functions .. + REAL SLAMCH, SLANGT, SLANHS + EXTERNAL SLAMCH, SLANGT, SLANHS +* .. +* .. Intrinsic Functions .. + INTRINSIC MIN +* .. +* .. External Subroutines .. + EXTERNAL SAXPY, SSWAP +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.LE.0 ) THEN + RESID = ZERO + RETURN + END IF +* + EPS = SLAMCH( 'Epsilon' ) +* +* Copy the matrix U to WORK. +* + DO 20 J = 1, N + DO 10 I = 1, N + WORK( I, J ) = ZERO + 10 CONTINUE + 20 CONTINUE + DO 30 I = 1, N + IF( I.EQ.1 ) THEN + WORK( I, I ) = DF( I ) + IF( N.GE.2 ) + $ WORK( I, I+1 ) = DUF( I ) + IF( N.GE.3 ) + $ WORK( I, I+2 ) = DU2( I ) + ELSE IF( I.EQ.N ) THEN + WORK( I, I ) = DF( I ) + ELSE + WORK( I, I ) = DF( I ) + WORK( I, I+1 ) = DUF( I ) + IF( I.LT.N-1 ) + $ WORK( I, I+2 ) = DU2( I ) + END IF + 30 CONTINUE +* +* Multiply on the left by L. +* + LASTJ = N + DO 40 I = N - 1, 1, -1 + LI = DLF( I ) + CALL SAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK, + $ WORK( I+1, I ), LDWORK ) + IP = IPIV( I ) + IF( IP.EQ.I ) THEN + LASTJ = MIN( I+2, N ) + ELSE + CALL SSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ), + $ LDWORK ) + END IF + 40 CONTINUE +* +* Subtract the matrix A. +* + WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 ) + IF( N.GT.1 ) THEN + WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 ) + WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 ) + WORK( N, N ) = WORK( N, N ) - D( N ) + DO 50 I = 2, N - 1 + WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 ) + WORK( I, I ) = WORK( I, I ) - D( I ) + WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I ) + 50 CONTINUE + END IF +* +* Compute the 1-norm of the tridiagonal matrix A. +* + ANORM = SLANGT( '1', N, DL, D, DU ) +* +* Compute the 1-norm of WORK, which is only guaranteed to be +* upper Hessenberg. +* + RESID = SLANHS( '1', N, WORK, LDWORK, RWORK ) +* +* Compute norm(L*U - A) / (norm(A) * EPS) +* + IF( ANORM.LE.ZERO ) THEN + IF( RESID.NE.ZERO ) + $ RESID = ONE / EPS + ELSE + RESID = ( RESID / ANORM ) / EPS + END IF +* + RETURN +* +* End of SGTT01 +* + END |