1*> \brief \b CGTT01
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 CGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK,
12*                          LDWORK, RWORK, RESID )
13*
14*       .. Scalar Arguments ..
15*       INTEGER            LDWORK, N
16*       REAL               RESID
17*       ..
18*       .. Array Arguments ..
19*       INTEGER            IPIV( * )
20*       REAL               RWORK( * )
21*       COMPLEX            D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
22*      $                   DU2( * ), DUF( * ), WORK( LDWORK, * )
23*       ..
24*
25*
26*> \par Purpose:
27*  =============
28*>
29*> \verbatim
30*>
31*> CGTT01 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 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 array, dimension (N)
55*>          The diagonal elements of A.
56*> \endverbatim
57*>
58*> \param[in] DU
59*> \verbatim
60*>          DU is COMPLEX 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 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 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 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 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 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 REAL array, dimension (N)
113*> \endverbatim
114*>
115*> \param[out] RESID
116*> \verbatim
117*>          RESID is REAL
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 complex_lin
130*
131*  =====================================================================
132      SUBROUTINE CGTT01( 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      REAL               RESID
142*     ..
143*     .. Array Arguments ..
144      INTEGER            IPIV( * )
145      REAL               RWORK( * )
146      COMPLEX            D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
147     $                   DU2( * ), DUF( * ), WORK( LDWORK, * )
148*     ..
149*
150*  =====================================================================
151*
152*     .. Parameters ..
153      REAL               ONE, ZERO
154      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
155*     ..
156*     .. Local Scalars ..
157      INTEGER            I, IP, J, LASTJ
158      REAL               ANORM, EPS
159      COMPLEX            LI
160*     ..
161*     .. External Functions ..
162      REAL               CLANGT, CLANHS, SLAMCH
163      EXTERNAL           CLANGT, CLANHS, SLAMCH
164*     ..
165*     .. Intrinsic Functions ..
166      INTRINSIC          MIN
167*     ..
168*     .. External Subroutines ..
169      EXTERNAL           CAXPY, CSWAP
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 = SLAMCH( '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 CAXPY( 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 CSWAP( 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 = CLANGT( '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 = CLANHS( '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 CGTT01
257*
258      END
259