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 INTEGTER
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*> \date December 2016
130*
131*> \ingroup complex16_lin
132*
133*  =====================================================================
134      SUBROUTINE ZGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK,
135     $                   LDWORK, RWORK, RESID )
136*
137*  -- LAPACK test routine (version 3.7.0) --
138*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
139*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140*     December 2016
141*
142*     .. Scalar Arguments ..
143      INTEGER            LDWORK, N
144      DOUBLE PRECISION   RESID
145*     ..
146*     .. Array Arguments ..
147      INTEGER            IPIV( * )
148      DOUBLE PRECISION   RWORK( * )
149      COMPLEX*16         D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
150     $                   DU2( * ), DUF( * ), WORK( LDWORK, * )
151*     ..
152*
153*  =====================================================================
154*
155*     .. Parameters ..
156      DOUBLE PRECISION   ONE, ZERO
157      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
158*     ..
159*     .. Local Scalars ..
160      INTEGER            I, IP, J, LASTJ
161      DOUBLE PRECISION   ANORM, EPS
162      COMPLEX*16         LI
163*     ..
164*     .. External Functions ..
165      DOUBLE PRECISION   DLAMCH, ZLANGT, ZLANHS
166      EXTERNAL           DLAMCH, ZLANGT, ZLANHS
167*     ..
168*     .. Intrinsic Functions ..
169      INTRINSIC          MIN
170*     ..
171*     .. External Subroutines ..
172      EXTERNAL           ZAXPY, ZSWAP
173*     ..
174*     .. Executable Statements ..
175*
176*     Quick return if possible
177*
178      IF( N.LE.0 ) THEN
179         RESID = ZERO
180         RETURN
181      END IF
182*
183      EPS = DLAMCH( 'Epsilon' )
184*
185*     Copy the matrix U to WORK.
186*
187      DO 20 J = 1, N
188         DO 10 I = 1, N
189            WORK( I, J ) = ZERO
190   10    CONTINUE
191   20 CONTINUE
192      DO 30 I = 1, N
193         IF( I.EQ.1 ) THEN
194            WORK( I, I ) = DF( I )
195            IF( N.GE.2 )
196     $         WORK( I, I+1 ) = DUF( I )
197            IF( N.GE.3 )
198     $         WORK( I, I+2 ) = DU2( I )
199         ELSE IF( I.EQ.N ) THEN
200            WORK( I, I ) = DF( I )
201         ELSE
202            WORK( I, I ) = DF( I )
203            WORK( I, I+1 ) = DUF( I )
204            IF( I.LT.N-1 )
205     $         WORK( I, I+2 ) = DU2( I )
206         END IF
207   30 CONTINUE
208*
209*     Multiply on the left by L.
210*
211      LASTJ = N
212      DO 40 I = N - 1, 1, -1
213         LI = DLF( I )
214         CALL ZAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK,
215     $               WORK( I+1, I ), LDWORK )
216         IP = IPIV( I )
217         IF( IP.EQ.I ) THEN
218            LASTJ = MIN( I+2, N )
219         ELSE
220            CALL ZSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ),
221     $                  LDWORK )
222         END IF
223   40 CONTINUE
224*
225*     Subtract the matrix A.
226*
227      WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 )
228      IF( N.GT.1 ) THEN
229         WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 )
230         WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 )
231         WORK( N, N ) = WORK( N, N ) - D( N )
232         DO 50 I = 2, N - 1
233            WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 )
234            WORK( I, I ) = WORK( I, I ) - D( I )
235            WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I )
236   50    CONTINUE
237      END IF
238*
239*     Compute the 1-norm of the tridiagonal matrix A.
240*
241      ANORM = ZLANGT( '1', N, DL, D, DU )
242*
243*     Compute the 1-norm of WORK, which is only guaranteed to be
244*     upper Hessenberg.
245*
246      RESID = ZLANHS( '1', N, WORK, LDWORK, RWORK )
247*
248*     Compute norm(L*U - A) / (norm(A) * EPS)
249*
250      IF( ANORM.LE.ZERO ) THEN
251         IF( RESID.NE.ZERO )
252     $      RESID = ONE / EPS
253      ELSE
254         RESID = ( RESID / ANORM ) / EPS
255      END IF
256*
257      RETURN
258*
259*     End of ZGTT01
260*
261      END
262