LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zptt01.f
Go to the documentation of this file.
1*> \brief \b ZPTT01
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 ZPTT01( N, D, E, DF, EF, WORK, RESID )
12*
13* .. Scalar Arguments ..
14* INTEGER N
15* DOUBLE PRECISION RESID
16* ..
17* .. Array Arguments ..
18* DOUBLE PRECISION D( * ), DF( * )
19* COMPLEX*16 E( * ), EF( * ), WORK( * )
20* ..
21*
22*
23*> \par Purpose:
24* =============
25*>
26*> \verbatim
27*>
28*> ZPTT01 reconstructs a tridiagonal matrix A from its L*D*L'
29*> factorization and computes the residual
30*> norm(L*D*L' - A) / ( n * norm(A) * EPS ),
31*> where EPS is the machine epsilon.
32*> \endverbatim
33*
34* Arguments:
35* ==========
36*
37*> \param[in] N
38*> \verbatim
39*> N is INTEGER
40*> The order of the matrix A.
41*> \endverbatim
42*>
43*> \param[in] D
44*> \verbatim
45*> D is DOUBLE PRECISION array, dimension (N)
46*> The n diagonal elements of the tridiagonal matrix A.
47*> \endverbatim
48*>
49*> \param[in] E
50*> \verbatim
51*> E is COMPLEX*16 array, dimension (N-1)
52*> The (n-1) subdiagonal elements of the tridiagonal matrix A.
53*> \endverbatim
54*>
55*> \param[in] DF
56*> \verbatim
57*> DF is DOUBLE PRECISION array, dimension (N)
58*> The n diagonal elements of the factor L from the L*D*L'
59*> factorization of A.
60*> \endverbatim
61*>
62*> \param[in] EF
63*> \verbatim
64*> EF is COMPLEX*16 array, dimension (N-1)
65*> The (n-1) subdiagonal elements of the factor L from the
66*> L*D*L' factorization of A.
67*> \endverbatim
68*>
69*> \param[out] WORK
70*> \verbatim
71*> WORK is COMPLEX*16 array, dimension (2*N)
72*> \endverbatim
73*>
74*> \param[out] RESID
75*> \verbatim
76*> RESID is DOUBLE PRECISION
77*> norm(L*D*L' - A) / (n * norm(A) * EPS)
78*> \endverbatim
79*
80* Authors:
81* ========
82*
83*> \author Univ. of Tennessee
84*> \author Univ. of California Berkeley
85*> \author Univ. of Colorado Denver
86*> \author NAG Ltd.
87*
88*> \ingroup complex16_lin
89*
90* =====================================================================
91 SUBROUTINE zptt01( N, D, E, DF, EF, WORK, RESID )
92*
93* -- LAPACK test routine --
94* -- LAPACK is a software package provided by Univ. of Tennessee, --
95* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
96*
97* .. Scalar Arguments ..
98 INTEGER N
99 DOUBLE PRECISION RESID
100* ..
101* .. Array Arguments ..
102 DOUBLE PRECISION D( * ), DF( * )
103 COMPLEX*16 E( * ), EF( * ), WORK( * )
104* ..
105*
106* =====================================================================
107*
108* .. Parameters ..
109 DOUBLE PRECISION ONE, ZERO
110 parameter( one = 1.0d+0, zero = 0.0d+0 )
111* ..
112* .. Local Scalars ..
113 INTEGER I
114 DOUBLE PRECISION ANORM, EPS
115 COMPLEX*16 DE
116* ..
117* .. External Functions ..
118 DOUBLE PRECISION DLAMCH
119 EXTERNAL dlamch
120* ..
121* .. Intrinsic Functions ..
122 INTRINSIC abs, dble, dconjg, max
123* ..
124* .. Executable Statements ..
125*
126* Quick return if possible
127*
128 IF( n.LE.0 ) THEN
129 resid = zero
130 RETURN
131 END IF
132*
133 eps = dlamch( 'Epsilon' )
134*
135* Construct the difference L*D*L' - A.
136*
137 work( 1 ) = df( 1 ) - d( 1 )
138 DO 10 i = 1, n - 1
139 de = df( i )*ef( i )
140 work( n+i ) = de - e( i )
141 work( 1+i ) = de*dconjg( ef( i ) ) + df( i+1 ) - d( i+1 )
142 10 CONTINUE
143*
144* Compute the 1-norms of the tridiagonal matrices A and WORK.
145*
146 IF( n.EQ.1 ) THEN
147 anorm = d( 1 )
148 resid = abs( work( 1 ) )
149 ELSE
150 anorm = max( d( 1 )+abs( e( 1 ) ), d( n )+abs( e( n-1 ) ) )
151 resid = max( abs( work( 1 ) )+abs( work( n+1 ) ),
152 $ abs( work( n ) )+abs( work( 2*n-1 ) ) )
153 DO 20 i = 2, n - 1
154 anorm = max( anorm, d( i )+abs( e( i ) )+abs( e( i-1 ) ) )
155 resid = max( resid, abs( work( i ) )+abs( work( n+i-1 ) )+
156 $ abs( work( n+i ) ) )
157 20 CONTINUE
158 END IF
159*
160* Compute norm(L*D*L' - A) / (n * norm(A) * EPS)
161*
162 IF( anorm.LE.zero ) THEN
163 IF( resid.NE.zero )
164 $ resid = one / eps
165 ELSE
166 resid = ( ( resid / dble( n ) ) / anorm ) / eps
167 END IF
168*
169 RETURN
170*
171* End of ZPTT01
172*
173 END
subroutine zptt01(n, d, e, df, ef, work, resid)
ZPTT01
Definition zptt01.f:92