1*> \brief \b SGTT01
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 SGTT01( 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               D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
21*      $                   DU2( * ), DUF( * ), RWORK( * ),
22*      $                   WORK( LDWORK, * )
23*       ..
24*
25*
26*> \par Purpose:
27*  =============
28*>
29*> \verbatim
30*>
31*> SGTT01 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 REAL 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 REAL array, dimension (N)
55*>          The diagonal elements of A.
56*> \endverbatim
57*>
58*> \param[in] DU
59*> \verbatim
60*>          DU is REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 single_lin
130*
131*  =====================================================================
132      SUBROUTINE SGTT01( 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               D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
146     $                   DU2( * ), DUF( * ), RWORK( * ),
147     $                   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, LI
159*     ..
160*     .. External Functions ..
161      REAL               SLAMCH, SLANGT, SLANHS
162      EXTERNAL           SLAMCH, SLANGT, SLANHS
163*     ..
164*     .. Intrinsic Functions ..
165      INTRINSIC          MIN
166*     ..
167*     .. External Subroutines ..
168      EXTERNAL           SAXPY, SSWAP
169*     ..
170*     .. Executable Statements ..
171*
172*     Quick return if possible
173*
174      IF( N.LE.0 ) THEN
175         RESID = ZERO
176         RETURN
177      END IF
178*
179      EPS = SLAMCH( 'Epsilon' )
180*
181*     Copy the matrix U to WORK.
182*
183      DO 20 J = 1, N
184         DO 10 I = 1, N
185            WORK( I, J ) = ZERO
186   10    CONTINUE
187   20 CONTINUE
188      DO 30 I = 1, N
189         IF( I.EQ.1 ) THEN
190            WORK( I, I ) = DF( I )
191            IF( N.GE.2 )
192     $         WORK( I, I+1 ) = DUF( I )
193            IF( N.GE.3 )
194     $         WORK( I, I+2 ) = DU2( I )
195         ELSE IF( I.EQ.N ) THEN
196            WORK( I, I ) = DF( I )
197         ELSE
198            WORK( I, I ) = DF( I )
199            WORK( I, I+1 ) = DUF( I )
200            IF( I.LT.N-1 )
201     $         WORK( I, I+2 ) = DU2( I )
202         END IF
203   30 CONTINUE
204*
205*     Multiply on the left by L.
206*
207      LASTJ = N
208      DO 40 I = N - 1, 1, -1
209         LI = DLF( I )
210         CALL SAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK,
211     $               WORK( I+1, I ), LDWORK )
212         IP = IPIV( I )
213         IF( IP.EQ.I ) THEN
214            LASTJ = MIN( I+2, N )
215         ELSE
216            CALL SSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ),
217     $                  LDWORK )
218         END IF
219   40 CONTINUE
220*
221*     Subtract the matrix A.
222*
223      WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 )
224      IF( N.GT.1 ) THEN
225         WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 )
226         WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 )
227         WORK( N, N ) = WORK( N, N ) - D( N )
228         DO 50 I = 2, N - 1
229            WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 )
230            WORK( I, I ) = WORK( I, I ) - D( I )
231            WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I )
232   50    CONTINUE
233      END IF
234*
235*     Compute the 1-norm of the tridiagonal matrix A.
236*
237      ANORM = SLANGT( '1', N, DL, D, DU )
238*
239*     Compute the 1-norm of WORK, which is only guaranteed to be
240*     upper Hessenberg.
241*
242      RESID = SLANHS( '1', N, WORK, LDWORK, RWORK )
243*
244*     Compute norm(L*U - A) / (norm(A) * EPS)
245*
246      IF( ANORM.LE.ZERO ) THEN
247         IF( RESID.NE.ZERO )
248     $      RESID = ONE / EPS
249      ELSE
250         RESID = ( RESID / ANORM ) / EPS
251      END IF
252*
253      RETURN
254*
255*     End of SGTT01
256*
257      END
258