1 SUBROUTINE SGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK,
2 $ LDWORK, RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER LDWORK, N
10 REAL RESID
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 REAL D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
15 $ DU2( * ), DUF( * ), RWORK( * ),
16 $ WORK( LDWORK, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * SGTT01 reconstructs a tridiagonal matrix A from its LU factorization
23 * and computes the residual
24 * norm(L*U - A) / ( norm(A) * EPS ),
25 * where EPS is the machine epsilon.
26 *
27 * Arguments
28 * =========
29 *
30 * N (input) INTEGTER
31 * The order of the matrix A. N >= 0.
32 *
33 * DL (input) REAL array, dimension (N-1)
34 * The (n-1) sub-diagonal elements of A.
35 *
36 * D (input) REAL array, dimension (N)
37 * The diagonal elements of A.
38 *
39 * DU (input) REAL array, dimension (N-1)
40 * The (n-1) super-diagonal elements of A.
41 *
42 * DLF (input) REAL array, dimension (N-1)
43 * The (n-1) multipliers that define the matrix L from the
44 * LU factorization of A.
45 *
46 * DF (input) REAL array, dimension (N)
47 * The n diagonal elements of the upper triangular matrix U from
48 * the LU factorization of A.
49 *
50 * DUF (input) REAL array, dimension (N-1)
51 * The (n-1) elements of the first super-diagonal of U.
52 *
53 * DU2F (input) REAL array, dimension (N-2)
54 * The (n-2) elements of the second super-diagonal of U.
55 *
56 * IPIV (input) INTEGER array, dimension (N)
57 * The pivot indices; for 1 <= i <= n, row i of the matrix was
58 * interchanged with row IPIV(i). IPIV(i) will always be either
59 * i or i+1; IPIV(i) = i indicates a row interchange was not
60 * required.
61 *
62 * WORK (workspace) REAL array, dimension (LDWORK,N)
63 *
64 * LDWORK (input) INTEGER
65 * The leading dimension of the array WORK. LDWORK >= max(1,N).
66 *
67 * RWORK (workspace) REAL array, dimension (N)
68 *
69 * RESID (output) REAL
70 * The scaled residual: norm(L*U - A) / (norm(A) * EPS)
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 REAL ONE, ZERO
76 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
77 * ..
78 * .. Local Scalars ..
79 INTEGER I, IP, J, LASTJ
80 REAL ANORM, EPS, LI
81 * ..
82 * .. External Functions ..
83 REAL SLAMCH, SLANGT, SLANHS
84 EXTERNAL SLAMCH, SLANGT, SLANHS
85 * ..
86 * .. Intrinsic Functions ..
87 INTRINSIC MIN
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL SAXPY, SSWAP
91 * ..
92 * .. Executable Statements ..
93 *
94 * Quick return if possible
95 *
96 IF( N.LE.0 ) THEN
97 RESID = ZERO
98 RETURN
99 END IF
100 *
101 EPS = SLAMCH( 'Epsilon' )
102 *
103 * Copy the matrix U to WORK.
104 *
105 DO 20 J = 1, N
106 DO 10 I = 1, N
107 WORK( I, J ) = ZERO
108 10 CONTINUE
109 20 CONTINUE
110 DO 30 I = 1, N
111 IF( I.EQ.1 ) THEN
112 WORK( I, I ) = DF( I )
113 IF( N.GE.2 )
114 $ WORK( I, I+1 ) = DUF( I )
115 IF( N.GE.3 )
116 $ WORK( I, I+2 ) = DU2( I )
117 ELSE IF( I.EQ.N ) THEN
118 WORK( I, I ) = DF( I )
119 ELSE
120 WORK( I, I ) = DF( I )
121 WORK( I, I+1 ) = DUF( I )
122 IF( I.LT.N-1 )
123 $ WORK( I, I+2 ) = DU2( I )
124 END IF
125 30 CONTINUE
126 *
127 * Multiply on the left by L.
128 *
129 LASTJ = N
130 DO 40 I = N - 1, 1, -1
131 LI = DLF( I )
132 CALL SAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK,
133 $ WORK( I+1, I ), LDWORK )
134 IP = IPIV( I )
135 IF( IP.EQ.I ) THEN
136 LASTJ = MIN( I+2, N )
137 ELSE
138 CALL SSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ),
139 $ LDWORK )
140 END IF
141 40 CONTINUE
142 *
143 * Subtract the matrix A.
144 *
145 WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 )
146 IF( N.GT.1 ) THEN
147 WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 )
148 WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 )
149 WORK( N, N ) = WORK( N, N ) - D( N )
150 DO 50 I = 2, N - 1
151 WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 )
152 WORK( I, I ) = WORK( I, I ) - D( I )
153 WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I )
154 50 CONTINUE
155 END IF
156 *
157 * Compute the 1-norm of the tridiagonal matrix A.
158 *
159 ANORM = SLANGT( '1', N, DL, D, DU )
160 *
161 * Compute the 1-norm of WORK, which is only guaranteed to be
162 * upper Hessenberg.
163 *
164 RESID = SLANHS( '1', N, WORK, LDWORK, RWORK )
165 *
166 * Compute norm(L*U - A) / (norm(A) * EPS)
167 *
168 IF( ANORM.LE.ZERO ) THEN
169 IF( RESID.NE.ZERO )
170 $ RESID = ONE / EPS
171 ELSE
172 RESID = ( RESID / ANORM ) / EPS
173 END IF
174 *
175 RETURN
176 *
177 * End of SGTT01
178 *
179 END
2 $ LDWORK, RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER LDWORK, N
10 REAL RESID
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 REAL D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
15 $ DU2( * ), DUF( * ), RWORK( * ),
16 $ WORK( LDWORK, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * SGTT01 reconstructs a tridiagonal matrix A from its LU factorization
23 * and computes the residual
24 * norm(L*U - A) / ( norm(A) * EPS ),
25 * where EPS is the machine epsilon.
26 *
27 * Arguments
28 * =========
29 *
30 * N (input) INTEGTER
31 * The order of the matrix A. N >= 0.
32 *
33 * DL (input) REAL array, dimension (N-1)
34 * The (n-1) sub-diagonal elements of A.
35 *
36 * D (input) REAL array, dimension (N)
37 * The diagonal elements of A.
38 *
39 * DU (input) REAL array, dimension (N-1)
40 * The (n-1) super-diagonal elements of A.
41 *
42 * DLF (input) REAL array, dimension (N-1)
43 * The (n-1) multipliers that define the matrix L from the
44 * LU factorization of A.
45 *
46 * DF (input) REAL array, dimension (N)
47 * The n diagonal elements of the upper triangular matrix U from
48 * the LU factorization of A.
49 *
50 * DUF (input) REAL array, dimension (N-1)
51 * The (n-1) elements of the first super-diagonal of U.
52 *
53 * DU2F (input) REAL array, dimension (N-2)
54 * The (n-2) elements of the second super-diagonal of U.
55 *
56 * IPIV (input) INTEGER array, dimension (N)
57 * The pivot indices; for 1 <= i <= n, row i of the matrix was
58 * interchanged with row IPIV(i). IPIV(i) will always be either
59 * i or i+1; IPIV(i) = i indicates a row interchange was not
60 * required.
61 *
62 * WORK (workspace) REAL array, dimension (LDWORK,N)
63 *
64 * LDWORK (input) INTEGER
65 * The leading dimension of the array WORK. LDWORK >= max(1,N).
66 *
67 * RWORK (workspace) REAL array, dimension (N)
68 *
69 * RESID (output) REAL
70 * The scaled residual: norm(L*U - A) / (norm(A) * EPS)
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 REAL ONE, ZERO
76 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
77 * ..
78 * .. Local Scalars ..
79 INTEGER I, IP, J, LASTJ
80 REAL ANORM, EPS, LI
81 * ..
82 * .. External Functions ..
83 REAL SLAMCH, SLANGT, SLANHS
84 EXTERNAL SLAMCH, SLANGT, SLANHS
85 * ..
86 * .. Intrinsic Functions ..
87 INTRINSIC MIN
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL SAXPY, SSWAP
91 * ..
92 * .. Executable Statements ..
93 *
94 * Quick return if possible
95 *
96 IF( N.LE.0 ) THEN
97 RESID = ZERO
98 RETURN
99 END IF
100 *
101 EPS = SLAMCH( 'Epsilon' )
102 *
103 * Copy the matrix U to WORK.
104 *
105 DO 20 J = 1, N
106 DO 10 I = 1, N
107 WORK( I, J ) = ZERO
108 10 CONTINUE
109 20 CONTINUE
110 DO 30 I = 1, N
111 IF( I.EQ.1 ) THEN
112 WORK( I, I ) = DF( I )
113 IF( N.GE.2 )
114 $ WORK( I, I+1 ) = DUF( I )
115 IF( N.GE.3 )
116 $ WORK( I, I+2 ) = DU2( I )
117 ELSE IF( I.EQ.N ) THEN
118 WORK( I, I ) = DF( I )
119 ELSE
120 WORK( I, I ) = DF( I )
121 WORK( I, I+1 ) = DUF( I )
122 IF( I.LT.N-1 )
123 $ WORK( I, I+2 ) = DU2( I )
124 END IF
125 30 CONTINUE
126 *
127 * Multiply on the left by L.
128 *
129 LASTJ = N
130 DO 40 I = N - 1, 1, -1
131 LI = DLF( I )
132 CALL SAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK,
133 $ WORK( I+1, I ), LDWORK )
134 IP = IPIV( I )
135 IF( IP.EQ.I ) THEN
136 LASTJ = MIN( I+2, N )
137 ELSE
138 CALL SSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ),
139 $ LDWORK )
140 END IF
141 40 CONTINUE
142 *
143 * Subtract the matrix A.
144 *
145 WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 )
146 IF( N.GT.1 ) THEN
147 WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 )
148 WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 )
149 WORK( N, N ) = WORK( N, N ) - D( N )
150 DO 50 I = 2, N - 1
151 WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 )
152 WORK( I, I ) = WORK( I, I ) - D( I )
153 WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I )
154 50 CONTINUE
155 END IF
156 *
157 * Compute the 1-norm of the tridiagonal matrix A.
158 *
159 ANORM = SLANGT( '1', N, DL, D, DU )
160 *
161 * Compute the 1-norm of WORK, which is only guaranteed to be
162 * upper Hessenberg.
163 *
164 RESID = SLANHS( '1', N, WORK, LDWORK, RWORK )
165 *
166 * Compute norm(L*U - A) / (norm(A) * EPS)
167 *
168 IF( ANORM.LE.ZERO ) THEN
169 IF( RESID.NE.ZERO )
170 $ RESID = ONE / EPS
171 ELSE
172 RESID = ( RESID / ANORM ) / EPS
173 END IF
174 *
175 RETURN
176 *
177 * End of SGTT01
178 *
179 END