1*> \brief \b CHBT21
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 CHBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK,
12*                          RWORK, RESULT )
13*
14*       .. Scalar Arguments ..
15*       CHARACTER          UPLO
16*       INTEGER            KA, KS, LDA, LDU, N
17*       ..
18*       .. Array Arguments ..
19*       REAL               D( * ), E( * ), RESULT( 2 ), RWORK( * )
20*       COMPLEX            A( LDA, * ), U( LDU, * ), WORK( * )
21*       ..
22*
23*
24*> \par Purpose:
25*  =============
26*>
27*> \verbatim
28*>
29*> CHBT21  generally checks a decomposition of the form
30*>
31*>         A = U S UC>
32*> where * means conjugate transpose, A is hermitian banded, U is
33*> unitary, and S is diagonal (if KS=0) or symmetric
34*> tridiagonal (if KS=1).
35*>
36*> Specifically:
37*>
38*>         RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC>         RESULT(2) = | I - UU* | / ( n ulp )
39*> \endverbatim
40*
41*  Arguments:
42*  ==========
43*
44*> \param[in] UPLO
45*> \verbatim
46*>          UPLO is CHARACTER
47*>          If UPLO='U', the upper triangle of A and V will be used and
48*>          the (strictly) lower triangle will not be referenced.
49*>          If UPLO='L', the lower triangle of A and V will be used and
50*>          the (strictly) upper triangle will not be referenced.
51*> \endverbatim
52*>
53*> \param[in] N
54*> \verbatim
55*>          N is INTEGER
56*>          The size of the matrix.  If it is zero, CHBT21 does nothing.
57*>          It must be at least zero.
58*> \endverbatim
59*>
60*> \param[in] KA
61*> \verbatim
62*>          KA is INTEGER
63*>          The bandwidth of the matrix A.  It must be at least zero.  If
64*>          it is larger than N-1, then max( 0, N-1 ) will be used.
65*> \endverbatim
66*>
67*> \param[in] KS
68*> \verbatim
69*>          KS is INTEGER
70*>          The bandwidth of the matrix S.  It may only be zero or one.
71*>          If zero, then S is diagonal, and E is not referenced.  If
72*>          one, then S is symmetric tri-diagonal.
73*> \endverbatim
74*>
75*> \param[in] A
76*> \verbatim
77*>          A is COMPLEX array, dimension (LDA, N)
78*>          The original (unfactored) matrix.  It is assumed to be
79*>          hermitian, and only the upper (UPLO='U') or only the lower
80*>          (UPLO='L') will be referenced.
81*> \endverbatim
82*>
83*> \param[in] LDA
84*> \verbatim
85*>          LDA is INTEGER
86*>          The leading dimension of A.  It must be at least 1
87*>          and at least min( KA, N-1 ).
88*> \endverbatim
89*>
90*> \param[in] D
91*> \verbatim
92*>          D is REAL array, dimension (N)
93*>          The diagonal of the (symmetric tri-) diagonal matrix S.
94*> \endverbatim
95*>
96*> \param[in] E
97*> \verbatim
98*>          E is REAL array, dimension (N-1)
99*>          The off-diagonal of the (symmetric tri-) diagonal matrix S.
100*>          E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
101*>          (3,2) element, etc.
102*>          Not referenced if KS=0.
103*> \endverbatim
104*>
105*> \param[in] U
106*> \verbatim
107*>          U is COMPLEX array, dimension (LDU, N)
108*>          The unitary matrix in the decomposition, expressed as a
109*>          dense matrix (i.e., not as a product of Householder
110*>          transformations, Givens transformations, etc.)
111*> \endverbatim
112*>
113*> \param[in] LDU
114*> \verbatim
115*>          LDU is INTEGER
116*>          The leading dimension of U.  LDU must be at least N and
117*>          at least 1.
118*> \endverbatim
119*>
120*> \param[out] WORK
121*> \verbatim
122*>          WORK is COMPLEX array, dimension (N**2)
123*> \endverbatim
124*>
125*> \param[out] RWORK
126*> \verbatim
127*>          RWORK is REAL array, dimension (N)
128*> \endverbatim
129*>
130*> \param[out] RESULT
131*> \verbatim
132*>          RESULT is REAL array, dimension (2)
133*>          The values computed by the two tests described above.  The
134*>          values are currently limited to 1/ulp, to avoid overflow.
135*> \endverbatim
136*
137*  Authors:
138*  ========
139*
140*> \author Univ. of Tennessee
141*> \author Univ. of California Berkeley
142*> \author Univ. of Colorado Denver
143*> \author NAG Ltd.
144*
145*> \date November 2011
146*
147*> \ingroup complex_eig
148*
149*  =====================================================================
150      SUBROUTINE CHBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK,
151     $                   RWORK, RESULT )
152*
153*  -- LAPACK test routine (version 3.4.0) --
154*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
155*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
156*     November 2011
157*
158*     .. Scalar Arguments ..
159      CHARACTER          UPLO
160      INTEGER            KA, KS, LDA, LDU, N
161*     ..
162*     .. Array Arguments ..
163      REAL               D( * ), E( * ), RESULT( 2 ), RWORK( * )
164      COMPLEX            A( LDA, * ), U( LDU, * ), WORK( * )
165*     ..
166*
167*  =====================================================================
168*
169*     .. Parameters ..
170      COMPLEX            CZERO, CONE
171      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
172     $                   CONE = ( 1.0E+0, 0.0E+0 ) )
173      REAL               ZERO, ONE
174      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
175*     ..
176*     .. Local Scalars ..
177      LOGICAL            LOWER
178      CHARACTER          CUPLO
179      INTEGER            IKA, J, JC, JR
180      REAL               ANORM, ULP, UNFL, WNORM
181*     ..
182*     .. External Functions ..
183      LOGICAL            LSAME
184      REAL               CLANGE, CLANHB, CLANHP, SLAMCH
185      EXTERNAL           LSAME, CLANGE, CLANHB, CLANHP, SLAMCH
186*     ..
187*     .. External Subroutines ..
188      EXTERNAL           CGEMM, CHPR, CHPR2
189*     ..
190*     .. Intrinsic Functions ..
191      INTRINSIC          CMPLX, MAX, MIN, REAL
192*     ..
193*     .. Executable Statements ..
194*
195*     Constants
196*
197      RESULT( 1 ) = ZERO
198      RESULT( 2 ) = ZERO
199      IF( N.LE.0 )
200     $   RETURN
201*
202      IKA = MAX( 0, MIN( N-1, KA ) )
203*
204      IF( LSAME( UPLO, 'U' ) ) THEN
205         LOWER = .FALSE.
206         CUPLO = 'U'
207      ELSE
208         LOWER = .TRUE.
209         CUPLO = 'L'
210      END IF
211*
212      UNFL = SLAMCH( 'Safe minimum' )
213      ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
214*
215*     Some Error Checks
216*
217*     Do Test 1
218*
219*     Norm of A:
220*
221      ANORM = MAX( CLANHB( '1', CUPLO, N, IKA, A, LDA, RWORK ), UNFL )
222*
223*     Compute error matrix:    Error = A - U S U*
224*
225*     Copy A from SB to SP storage format.
226*
227      J = 0
228      DO 50 JC = 1, N
229         IF( LOWER ) THEN
230            DO 10 JR = 1, MIN( IKA+1, N+1-JC )
231               J = J + 1
232               WORK( J ) = A( JR, JC )
233   10       CONTINUE
234            DO 20 JR = IKA + 2, N + 1 - JC
235               J = J + 1
236               WORK( J ) = ZERO
237   20       CONTINUE
238         ELSE
239            DO 30 JR = IKA + 2, JC
240               J = J + 1
241               WORK( J ) = ZERO
242   30       CONTINUE
243            DO 40 JR = MIN( IKA, JC-1 ), 0, -1
244               J = J + 1
245               WORK( J ) = A( IKA+1-JR, JC )
246   40       CONTINUE
247         END IF
248   50 CONTINUE
249*
250      DO 60 J = 1, N
251         CALL CHPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK )
252   60 CONTINUE
253*
254      IF( N.GT.1 .AND. KS.EQ.1 ) THEN
255         DO 70 J = 1, N - 1
256            CALL CHPR2( CUPLO, N, -CMPLX( E( J ) ), U( 1, J ), 1,
257     $                  U( 1, J+1 ), 1, WORK )
258   70    CONTINUE
259      END IF
260      WNORM = CLANHP( '1', CUPLO, N, WORK, RWORK )
261*
262      IF( ANORM.GT.WNORM ) THEN
263         RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
264      ELSE
265         IF( ANORM.LT.ONE ) THEN
266            RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
267         ELSE
268            RESULT( 1 ) = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP )
269         END IF
270      END IF
271*
272*     Do Test 2
273*
274*     Compute  UU* - I
275*
276      CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK,
277     $            N )
278*
279      DO 80 J = 1, N
280         WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE
281   80 CONTINUE
282*
283      RESULT( 2 ) = MIN( CLANGE( '1', N, N, WORK, N, RWORK ),
284     $              REAL( N ) ) / ( N*ULP )
285*
286      RETURN
287*
288*     End of CHBT21
289*
290      END
291