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