1 SUBROUTINE CGTT01( 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 RWORK( * )
15 COMPLEX D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
16 $ DU2( * ), DUF( * ), WORK( LDWORK, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * CGTT01 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) COMPLEX array, dimension (N-1)
34 * The (n-1) sub-diagonal elements of A.
35 *
36 * D (input) COMPLEX array, dimension (N)
37 * The diagonal elements of A.
38 *
39 * DU (input) COMPLEX array, dimension (N-1)
40 * The (n-1) super-diagonal elements of A.
41 *
42 * DLF (input) COMPLEX 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) COMPLEX 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) COMPLEX array, dimension (N-1)
51 * The (n-1) elements of the first super-diagonal of U.
52 *
53 * DU2 (input) COMPLEX 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) COMPLEX 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
81 COMPLEX LI
82 * ..
83 * .. External Functions ..
84 REAL CLANGT, CLANHS, SLAMCH
85 EXTERNAL CLANGT, CLANHS, SLAMCH
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC MIN
89 * ..
90 * .. External Subroutines ..
91 EXTERNAL CAXPY, CSWAP
92 * ..
93 * .. Executable Statements ..
94 *
95 * Quick return if possible
96 *
97 IF( N.LE.0 ) THEN
98 RESID = ZERO
99 RETURN
100 END IF
101 *
102 EPS = SLAMCH( 'Epsilon' )
103 *
104 * Copy the matrix U to WORK.
105 *
106 DO 20 J = 1, N
107 DO 10 I = 1, N
108 WORK( I, J ) = ZERO
109 10 CONTINUE
110 20 CONTINUE
111 DO 30 I = 1, N
112 IF( I.EQ.1 ) THEN
113 WORK( I, I ) = DF( I )
114 IF( N.GE.2 )
115 $ WORK( I, I+1 ) = DUF( I )
116 IF( N.GE.3 )
117 $ WORK( I, I+2 ) = DU2( I )
118 ELSE IF( I.EQ.N ) THEN
119 WORK( I, I ) = DF( I )
120 ELSE
121 WORK( I, I ) = DF( I )
122 WORK( I, I+1 ) = DUF( I )
123 IF( I.LT.N-1 )
124 $ WORK( I, I+2 ) = DU2( I )
125 END IF
126 30 CONTINUE
127 *
128 * Multiply on the left by L.
129 *
130 LASTJ = N
131 DO 40 I = N - 1, 1, -1
132 LI = DLF( I )
133 CALL CAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK,
134 $ WORK( I+1, I ), LDWORK )
135 IP = IPIV( I )
136 IF( IP.EQ.I ) THEN
137 LASTJ = MIN( I+2, N )
138 ELSE
139 CALL CSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ),
140 $ LDWORK )
141 END IF
142 40 CONTINUE
143 *
144 * Subtract the matrix A.
145 *
146 WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 )
147 IF( N.GT.1 ) THEN
148 WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 )
149 WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 )
150 WORK( N, N ) = WORK( N, N ) - D( N )
151 DO 50 I = 2, N - 1
152 WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 )
153 WORK( I, I ) = WORK( I, I ) - D( I )
154 WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I )
155 50 CONTINUE
156 END IF
157 *
158 * Compute the 1-norm of the tridiagonal matrix A.
159 *
160 ANORM = CLANGT( '1', N, DL, D, DU )
161 *
162 * Compute the 1-norm of WORK, which is only guaranteed to be
163 * upper Hessenberg.
164 *
165 RESID = CLANHS( '1', N, WORK, LDWORK, RWORK )
166 *
167 * Compute norm(L*U - A) / (norm(A) * EPS)
168 *
169 IF( ANORM.LE.ZERO ) THEN
170 IF( RESID.NE.ZERO )
171 $ RESID = ONE / EPS
172 ELSE
173 RESID = ( RESID / ANORM ) / EPS
174 END IF
175 *
176 RETURN
177 *
178 * End of CGTT01
179 *
180 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 RWORK( * )
15 COMPLEX D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
16 $ DU2( * ), DUF( * ), WORK( LDWORK, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * CGTT01 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) COMPLEX array, dimension (N-1)
34 * The (n-1) sub-diagonal elements of A.
35 *
36 * D (input) COMPLEX array, dimension (N)
37 * The diagonal elements of A.
38 *
39 * DU (input) COMPLEX array, dimension (N-1)
40 * The (n-1) super-diagonal elements of A.
41 *
42 * DLF (input) COMPLEX 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) COMPLEX 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) COMPLEX array, dimension (N-1)
51 * The (n-1) elements of the first super-diagonal of U.
52 *
53 * DU2 (input) COMPLEX 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) COMPLEX 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
81 COMPLEX LI
82 * ..
83 * .. External Functions ..
84 REAL CLANGT, CLANHS, SLAMCH
85 EXTERNAL CLANGT, CLANHS, SLAMCH
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC MIN
89 * ..
90 * .. External Subroutines ..
91 EXTERNAL CAXPY, CSWAP
92 * ..
93 * .. Executable Statements ..
94 *
95 * Quick return if possible
96 *
97 IF( N.LE.0 ) THEN
98 RESID = ZERO
99 RETURN
100 END IF
101 *
102 EPS = SLAMCH( 'Epsilon' )
103 *
104 * Copy the matrix U to WORK.
105 *
106 DO 20 J = 1, N
107 DO 10 I = 1, N
108 WORK( I, J ) = ZERO
109 10 CONTINUE
110 20 CONTINUE
111 DO 30 I = 1, N
112 IF( I.EQ.1 ) THEN
113 WORK( I, I ) = DF( I )
114 IF( N.GE.2 )
115 $ WORK( I, I+1 ) = DUF( I )
116 IF( N.GE.3 )
117 $ WORK( I, I+2 ) = DU2( I )
118 ELSE IF( I.EQ.N ) THEN
119 WORK( I, I ) = DF( I )
120 ELSE
121 WORK( I, I ) = DF( I )
122 WORK( I, I+1 ) = DUF( I )
123 IF( I.LT.N-1 )
124 $ WORK( I, I+2 ) = DU2( I )
125 END IF
126 30 CONTINUE
127 *
128 * Multiply on the left by L.
129 *
130 LASTJ = N
131 DO 40 I = N - 1, 1, -1
132 LI = DLF( I )
133 CALL CAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK,
134 $ WORK( I+1, I ), LDWORK )
135 IP = IPIV( I )
136 IF( IP.EQ.I ) THEN
137 LASTJ = MIN( I+2, N )
138 ELSE
139 CALL CSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ),
140 $ LDWORK )
141 END IF
142 40 CONTINUE
143 *
144 * Subtract the matrix A.
145 *
146 WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 )
147 IF( N.GT.1 ) THEN
148 WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 )
149 WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 )
150 WORK( N, N ) = WORK( N, N ) - D( N )
151 DO 50 I = 2, N - 1
152 WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 )
153 WORK( I, I ) = WORK( I, I ) - D( I )
154 WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I )
155 50 CONTINUE
156 END IF
157 *
158 * Compute the 1-norm of the tridiagonal matrix A.
159 *
160 ANORM = CLANGT( '1', N, DL, D, DU )
161 *
162 * Compute the 1-norm of WORK, which is only guaranteed to be
163 * upper Hessenberg.
164 *
165 RESID = CLANHS( '1', N, WORK, LDWORK, RWORK )
166 *
167 * Compute norm(L*U - A) / (norm(A) * EPS)
168 *
169 IF( ANORM.LE.ZERO ) THEN
170 IF( RESID.NE.ZERO )
171 $ RESID = ONE / EPS
172 ELSE
173 RESID = ( RESID / ANORM ) / EPS
174 END IF
175 *
176 RETURN
177 *
178 * End of CGTT01
179 *
180 END