LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
zstt21.f
Go to the documentation of this file.
1 *> \brief \b ZSTT21
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 ZSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RWORK,
12 * RESULT )
13 *
14 * .. Scalar Arguments ..
15 * INTEGER KBAND, LDU, N
16 * ..
17 * .. Array Arguments ..
18 * DOUBLE PRECISION AD( * ), AE( * ), RESULT( 2 ), RWORK( * ),
19 * $ SD( * ), SE( * )
20 * COMPLEX*16 U( LDU, * ), WORK( * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> ZSTT21 checks a decomposition of the form
30 *>
31 *> A = U S UC>
32 *> where * means conjugate transpose, A is real symmetric tridiagonal,
33 *> U is unitary, and S is real and diagonal (if KBAND=0) or symmetric
34 *> tridiagonal (if KBAND=1). Two tests are performed:
35 *>
36 *> RESULT(1) = | A - U S U* | / ( |A| n ulp )
37 *>
38 *> RESULT(2) = | I - UU* | / ( n ulp )
39 *> \endverbatim
40 *
41 * Arguments:
42 * ==========
43 *
44 *> \param[in] N
45 *> \verbatim
46 *> N is INTEGER
47 *> The size of the matrix. If it is zero, ZSTT21 does nothing.
48 *> It must be at least zero.
49 *> \endverbatim
50 *>
51 *> \param[in] KBAND
52 *> \verbatim
53 *> KBAND is INTEGER
54 *> The bandwidth of the matrix S. It may only be zero or one.
55 *> If zero, then S is diagonal, and SE is not referenced. If
56 *> one, then S is symmetric tri-diagonal.
57 *> \endverbatim
58 *>
59 *> \param[in] AD
60 *> \verbatim
61 *> AD is DOUBLE PRECISION array, dimension (N)
62 *> The diagonal of the original (unfactored) matrix A. A is
63 *> assumed to be real symmetric tridiagonal.
64 *> \endverbatim
65 *>
66 *> \param[in] AE
67 *> \verbatim
68 *> AE is DOUBLE PRECISION array, dimension (N-1)
69 *> The off-diagonal of the original (unfactored) matrix A. A
70 *> is assumed to be symmetric tridiagonal. AE(1) is the (1,2)
71 *> and (2,1) element, AE(2) is the (2,3) and (3,2) element, etc.
72 *> \endverbatim
73 *>
74 *> \param[in] SD
75 *> \verbatim
76 *> SD is DOUBLE PRECISION array, dimension (N)
77 *> The diagonal of the real (symmetric tri-) diagonal matrix S.
78 *> \endverbatim
79 *>
80 *> \param[in] SE
81 *> \verbatim
82 *> SE is DOUBLE PRECISION array, dimension (N-1)
83 *> The off-diagonal of the (symmetric tri-) diagonal matrix S.
84 *> Not referenced if KBSND=0. If KBAND=1, then AE(1) is the
85 *> (1,2) and (2,1) element, SE(2) is the (2,3) and (3,2)
86 *> element, etc.
87 *> \endverbatim
88 *>
89 *> \param[in] U
90 *> \verbatim
91 *> U is COMPLEX*16 array, dimension (LDU, N)
92 *> The unitary matrix in the decomposition.
93 *> \endverbatim
94 *>
95 *> \param[in] LDU
96 *> \verbatim
97 *> LDU is INTEGER
98 *> The leading dimension of U. LDU must be at least N.
99 *> \endverbatim
100 *>
101 *> \param[out] WORK
102 *> \verbatim
103 *> WORK is COMPLEX*16 array, dimension (N**2)
104 *> \endverbatim
105 *>
106 *> \param[out] RWORK
107 *> \verbatim
108 *> RWORK is DOUBLE PRECISION array, dimension (N)
109 *> \endverbatim
110 *>
111 *> \param[out] RESULT
112 *> \verbatim
113 *> RESULT is DOUBLE PRECISION array, dimension (2)
114 *> The values computed by the two tests described above. The
115 *> values are currently limited to 1/ulp, to avoid overflow.
116 *> RESULT(1) is always modified.
117 *> \endverbatim
118 *
119 * Authors:
120 * ========
121 *
122 *> \author Univ. of Tennessee
123 *> \author Univ. of California Berkeley
124 *> \author Univ. of Colorado Denver
125 *> \author NAG Ltd.
126 *
127 *> \date November 2011
128 *
129 *> \ingroup complex16_eig
130 *
131 * =====================================================================
132  SUBROUTINE zstt21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RWORK,
133  $ result )
134 *
135 * -- LAPACK test routine (version 3.4.0) --
136 * -- LAPACK is a software package provided by Univ. of Tennessee, --
137 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138 * November 2011
139 *
140 * .. Scalar Arguments ..
141  INTEGER kband, ldu, n
142 * ..
143 * .. Array Arguments ..
144  DOUBLE PRECISION ad( * ), ae( * ), result( 2 ), rwork( * ),
145  $ sd( * ), se( * )
146  COMPLEX*16 u( ldu, * ), work( * )
147 * ..
148 *
149 * =====================================================================
150 *
151 * .. Parameters ..
152  DOUBLE PRECISION zero, one
153  parameter( zero = 0.0d+0, one = 1.0d+0 )
154  COMPLEX*16 czero, cone
155  parameter( czero = ( 0.0d+0, 0.0d+0 ),
156  $ cone = ( 1.0d+0, 0.0d+0 ) )
157 * ..
158 * .. Local Scalars ..
159  INTEGER j
160  DOUBLE PRECISION anorm, temp1, temp2, ulp, unfl, wnorm
161 * ..
162 * .. External Functions ..
163  DOUBLE PRECISION dlamch, zlange, zlanhe
164  EXTERNAL dlamch, zlange, zlanhe
165 * ..
166 * .. External Subroutines ..
167  EXTERNAL zgemm, zher, zher2, zlaset
168 * ..
169 * .. Intrinsic Functions ..
170  INTRINSIC abs, dble, dcmplx, max, min
171 * ..
172 * .. Executable Statements ..
173 *
174 * 1) Constants
175 *
176  result( 1 ) = zero
177  result( 2 ) = zero
178  IF( n.LE.0 )
179  $ RETURN
180 *
181  unfl = dlamch( 'Safe minimum' )
182  ulp = dlamch( 'Precision' )
183 *
184 * Do Test 1
185 *
186 * Copy A & Compute its 1-Norm:
187 *
188  CALL zlaset( 'Full', n, n, czero, czero, work, n )
189 *
190  anorm = zero
191  temp1 = zero
192 *
193  DO 10 j = 1, n - 1
194  work( ( n+1 )*( j-1 )+1 ) = ad( j )
195  work( ( n+1 )*( j-1 )+2 ) = ae( j )
196  temp2 = abs( ae( j ) )
197  anorm = max( anorm, abs( ad( j ) )+temp1+temp2 )
198  temp1 = temp2
199  10 CONTINUE
200 *
201  work( n**2 ) = ad( n )
202  anorm = max( anorm, abs( ad( n ) )+temp1, unfl )
203 *
204 * Norm of A - USU*
205 *
206  DO 20 j = 1, n
207  CALL zher( 'L', n, -sd( j ), u( 1, j ), 1, work, n )
208  20 CONTINUE
209 *
210  IF( n.GT.1 .AND. kband.EQ.1 ) THEN
211  DO 30 j = 1, n - 1
212  CALL zher2( 'L', n, -dcmplx( se( j ) ), u( 1, j ), 1,
213  $ u( 1, j+1 ), 1, work, n )
214  30 CONTINUE
215  END IF
216 *
217  wnorm = zlanhe( '1', 'L', n, work, n, rwork )
218 *
219  IF( anorm.GT.wnorm ) THEN
220  result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
221  ELSE
222  IF( anorm.LT.one ) THEN
223  result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
224  ELSE
225  result( 1 ) = min( wnorm / anorm, dble( n ) ) / ( n*ulp )
226  END IF
227  END IF
228 *
229 * Do Test 2
230 *
231 * Compute UU* - I
232 *
233  CALL zgemm( 'N', 'C', n, n, n, cone, u, ldu, u, ldu, czero, work,
234  $ n )
235 *
236  DO 40 j = 1, n
237  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
238  40 CONTINUE
239 *
240  result( 2 ) = min( dble( n ), zlange( '1', n, n, work, n,
241  $ rwork ) ) / ( n*ulp )
242 *
243  RETURN
244 *
245 * End of ZSTT21
246 *
247  END