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