LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 INTEGTER
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 *> \date November 2011
89 *
90 *> \ingroup complex16_lin
91 *
92 * =====================================================================
93  SUBROUTINE zptt01( N, D, E, DF, EF, WORK, RESID )
94 *
95 * -- LAPACK test routine (version 3.4.0) --
96 * -- LAPACK is a software package provided by Univ. of Tennessee, --
97 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
98 * November 2011
99 *
100 * .. Scalar Arguments ..
101  INTEGER n
102  DOUBLE PRECISION resid
103 * ..
104 * .. Array Arguments ..
105  DOUBLE PRECISION d( * ), df( * )
106  COMPLEX*16 e( * ), ef( * ), work( * )
107 * ..
108 *
109 * =====================================================================
110 *
111 * .. Parameters ..
112  DOUBLE PRECISION one, zero
113  parameter( one = 1.0d+0, zero = 0.0d+0 )
114 * ..
115 * .. Local Scalars ..
116  INTEGER i
117  DOUBLE PRECISION anorm, eps
118  COMPLEX*16 de
119 * ..
120 * .. External Functions ..
121  DOUBLE PRECISION dlamch
122  EXTERNAL dlamch
123 * ..
124 * .. Intrinsic Functions ..
125  INTRINSIC abs, dble, dconjg, max
126 * ..
127 * .. Executable Statements ..
128 *
129 * Quick return if possible
130 *
131  IF( n.LE.0 ) THEN
132  resid = zero
133  return
134  END IF
135 *
136  eps = dlamch( 'Epsilon' )
137 *
138 * Construct the difference L*D*L' - A.
139 *
140  work( 1 ) = df( 1 ) - d( 1 )
141  DO 10 i = 1, n - 1
142  de = df( i )*ef( i )
143  work( n+i ) = de - e( i )
144  work( 1+i ) = de*dconjg( ef( i ) ) + df( i+1 ) - d( i+1 )
145  10 continue
146 *
147 * Compute the 1-norms of the tridiagonal matrices A and WORK.
148 *
149  IF( n.EQ.1 ) THEN
150  anorm = d( 1 )
151  resid = abs( work( 1 ) )
152  ELSE
153  anorm = max( d( 1 )+abs( e( 1 ) ), d( n )+abs( e( n-1 ) ) )
154  resid = max( abs( work( 1 ) )+abs( work( n+1 ) ),
155  $ abs( work( n ) )+abs( work( 2*n-1 ) ) )
156  DO 20 i = 2, n - 1
157  anorm = max( anorm, d( i )+abs( e( i ) )+abs( e( i-1 ) ) )
158  resid = max( resid, abs( work( i ) )+abs( work( n+i-1 ) )+
159  $ abs( work( n+i ) ) )
160  20 continue
161  END IF
162 *
163 * Compute norm(L*D*L' - A) / (n * norm(A) * EPS)
164 *
165  IF( anorm.LE.zero ) THEN
166  IF( resid.NE.zero )
167  $ resid = one / eps
168  ELSE
169  resid = ( ( resid / dble( n ) ) / anorm ) / eps
170  END IF
171 *
172  return
173 *
174 * End of ZPTT01
175 *
176  END