LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zgtt01.f
Go to the documentation of this file.
1*> \brief \b ZGTT01
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE ZGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK,
12* LDWORK, RWORK, RESID )
13*
14* .. Scalar Arguments ..
15* INTEGER LDWORK, N
16* DOUBLE PRECISION RESID
17* ..
18* .. Array Arguments ..
19* INTEGER IPIV( * )
20* DOUBLE PRECISION RWORK( * )
21* COMPLEX*16 D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
22* $ DU2( * ), DUF( * ), WORK( LDWORK, * )
23* ..
24*
25*
26*> \par Purpose:
27* =============
28*>
29*> \verbatim
30*>
31*> ZGTT01 reconstructs a tridiagonal matrix A from its LU factorization
32*> and computes the residual
33*> norm(L*U - A) / ( norm(A) * EPS ),
34*> where EPS is the machine epsilon.
35*> \endverbatim
36*
37* Arguments:
38* ==========
39*
40*> \param[in] N
41*> \verbatim
42*> N is INTEGER
43*> The order of the matrix A. N >= 0.
44*> \endverbatim
45*>
46*> \param[in] DL
47*> \verbatim
48*> DL is COMPLEX*16 array, dimension (N-1)
49*> The (n-1) sub-diagonal elements of A.
50*> \endverbatim
51*>
52*> \param[in] D
53*> \verbatim
54*> D is COMPLEX*16 array, dimension (N)
55*> The diagonal elements of A.
56*> \endverbatim
57*>
58*> \param[in] DU
59*> \verbatim
60*> DU is COMPLEX*16 array, dimension (N-1)
61*> The (n-1) super-diagonal elements of A.
62*> \endverbatim
63*>
64*> \param[in] DLF
65*> \verbatim
66*> DLF is COMPLEX*16 array, dimension (N-1)
67*> The (n-1) multipliers that define the matrix L from the
68*> LU factorization of A.
69*> \endverbatim
70*>
71*> \param[in] DF
72*> \verbatim
73*> DF is COMPLEX*16 array, dimension (N)
74*> The n diagonal elements of the upper triangular matrix U from
75*> the LU factorization of A.
76*> \endverbatim
77*>
78*> \param[in] DUF
79*> \verbatim
80*> DUF is COMPLEX*16 array, dimension (N-1)
81*> The (n-1) elements of the first super-diagonal of U.
82*> \endverbatim
83*>
84*> \param[in] DU2
85*> \verbatim
86*> DU2 is COMPLEX*16 array, dimension (N-2)
87*> The (n-2) elements of the second super-diagonal of U.
88*> \endverbatim
89*>
90*> \param[in] IPIV
91*> \verbatim
92*> IPIV is INTEGER array, dimension (N)
93*> The pivot indices; for 1 <= i <= n, row i of the matrix was
94*> interchanged with row IPIV(i). IPIV(i) will always be either
95*> i or i+1; IPIV(i) = i indicates a row interchange was not
96*> required.
97*> \endverbatim
98*>
99*> \param[out] WORK
100*> \verbatim
101*> WORK is COMPLEX*16 array, dimension (LDWORK,N)
102*> \endverbatim
103*>
104*> \param[in] LDWORK
105*> \verbatim
106*> LDWORK is INTEGER
107*> The leading dimension of the array WORK. LDWORK >= max(1,N).
108*> \endverbatim
109*>
110*> \param[out] RWORK
111*> \verbatim
112*> RWORK is DOUBLE PRECISION array, dimension (N)
113*> \endverbatim
114*>
115*> \param[out] RESID
116*> \verbatim
117*> RESID is DOUBLE PRECISION
118*> The scaled residual: norm(L*U - A) / (norm(A) * EPS)
119*> \endverbatim
120*
121* Authors:
122* ========
123*
124*> \author Univ. of Tennessee
125*> \author Univ. of California Berkeley
126*> \author Univ. of Colorado Denver
127*> \author NAG Ltd.
128*
129*> \ingroup complex16_lin
130*
131* =====================================================================
132 SUBROUTINE zgtt01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK,
133 $ LDWORK, RWORK, RESID )
134*
135* -- LAPACK test routine --
136* -- LAPACK is a software package provided by Univ. of Tennessee, --
137* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138*
139* .. Scalar Arguments ..
140 INTEGER LDWORK, N
141 DOUBLE PRECISION RESID
142* ..
143* .. Array Arguments ..
144 INTEGER IPIV( * )
145 DOUBLE PRECISION RWORK( * )
146 COMPLEX*16 D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
147 $ du2( * ), duf( * ), work( ldwork, * )
148* ..
149*
150* =====================================================================
151*
152* .. Parameters ..
153 DOUBLE PRECISION ONE, ZERO
154 parameter( one = 1.0d+0, zero = 0.0d+0 )
155* ..
156* .. Local Scalars ..
157 INTEGER I, IP, J, LASTJ
158 DOUBLE PRECISION ANORM, EPS
159 COMPLEX*16 LI
160* ..
161* .. External Functions ..
162 DOUBLE PRECISION DLAMCH, ZLANGT, ZLANHS
163 EXTERNAL dlamch, zlangt, zlanhs
164* ..
165* .. Intrinsic Functions ..
166 INTRINSIC min
167* ..
168* .. External Subroutines ..
169 EXTERNAL zaxpy, zswap
170* ..
171* .. Executable Statements ..
172*
173* Quick return if possible
174*
175 IF( n.LE.0 ) THEN
176 resid = zero
177 RETURN
178 END IF
179*
180 eps = dlamch( 'Epsilon' )
181*
182* Copy the matrix U to WORK.
183*
184 DO 20 j = 1, n
185 DO 10 i = 1, n
186 work( i, j ) = zero
187 10 CONTINUE
188 20 CONTINUE
189 DO 30 i = 1, n
190 IF( i.EQ.1 ) THEN
191 work( i, i ) = df( i )
192 IF( n.GE.2 )
193 $ work( i, i+1 ) = duf( i )
194 IF( n.GE.3 )
195 $ work( i, i+2 ) = du2( i )
196 ELSE IF( i.EQ.n ) THEN
197 work( i, i ) = df( i )
198 ELSE
199 work( i, i ) = df( i )
200 work( i, i+1 ) = duf( i )
201 IF( i.LT.n-1 )
202 $ work( i, i+2 ) = du2( i )
203 END IF
204 30 CONTINUE
205*
206* Multiply on the left by L.
207*
208 lastj = n
209 DO 40 i = n - 1, 1, -1
210 li = dlf( i )
211 CALL zaxpy( lastj-i+1, li, work( i, i ), ldwork,
212 $ work( i+1, i ), ldwork )
213 ip = ipiv( i )
214 IF( ip.EQ.i ) THEN
215 lastj = min( i+2, n )
216 ELSE
217 CALL zswap( lastj-i+1, work( i, i ), ldwork, work( i+1, i ),
218 $ ldwork )
219 END IF
220 40 CONTINUE
221*
222* Subtract the matrix A.
223*
224 work( 1, 1 ) = work( 1, 1 ) - d( 1 )
225 IF( n.GT.1 ) THEN
226 work( 1, 2 ) = work( 1, 2 ) - du( 1 )
227 work( n, n-1 ) = work( n, n-1 ) - dl( n-1 )
228 work( n, n ) = work( n, n ) - d( n )
229 DO 50 i = 2, n - 1
230 work( i, i-1 ) = work( i, i-1 ) - dl( i-1 )
231 work( i, i ) = work( i, i ) - d( i )
232 work( i, i+1 ) = work( i, i+1 ) - du( i )
233 50 CONTINUE
234 END IF
235*
236* Compute the 1-norm of the tridiagonal matrix A.
237*
238 anorm = zlangt( '1', n, dl, d, du )
239*
240* Compute the 1-norm of WORK, which is only guaranteed to be
241* upper Hessenberg.
242*
243 resid = zlanhs( '1', n, work, ldwork, rwork )
244*
245* Compute norm(L*U - A) / (norm(A) * EPS)
246*
247 IF( anorm.LE.zero ) THEN
248 IF( resid.NE.zero )
249 $ resid = one / eps
250 ELSE
251 resid = ( resid / anorm ) / eps
252 END IF
253*
254 RETURN
255*
256* End of ZGTT01
257*
258 END
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
subroutine zgtt01(n, dl, d, du, dlf, df, duf, du2, ipiv, work, ldwork, rwork, resid)
ZGTT01
Definition zgtt01.f:134