LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
cgtt01.f
Go to the documentation of this file.
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 *> \date December 2016
130 *
131 *> \ingroup complex_lin
132 *
133 * =====================================================================
134  SUBROUTINE cgtt01( 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  REAL RESID
145 * ..
146 * .. Array Arguments ..
147  INTEGER IPIV( * )
148  REAL RWORK( * )
149  COMPLEX D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
150  $ du2( * ), duf( * ), 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
162  COMPLEX LI
163 * ..
164 * .. External Functions ..
165  REAL CLANGT, CLANHS, SLAMCH
166  EXTERNAL clangt, clanhs, slamch
167 * ..
168 * .. Intrinsic Functions ..
169  INTRINSIC min
170 * ..
171 * .. External Subroutines ..
172  EXTERNAL caxpy, cswap
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 = slamch( '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 caxpy( 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 cswap( 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 = clangt( '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 = clanhs( '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 CGTT01
260 *
261  END
subroutine cgtt01(N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK, LDWORK, RWORK, RESID)
CGTT01
Definition: cgtt01.f:136
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:90