1*> \brief \b CGET51
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 CGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK,
12*                          RWORK, RESULT )
13*
14*       .. Scalar Arguments ..
15*       INTEGER            ITYPE, LDA, LDB, LDU, LDV, N
16*       REAL               RESULT
17*       ..
18*       .. Array Arguments ..
19*       REAL               RWORK( * )
20*       COMPLEX            A( LDA, * ), B( LDB, * ), U( LDU, * ),
21*      $                   V( LDV, * ), WORK( * )
22*       ..
23*
24*
25*> \par Purpose:
26*  =============
27*>
28*> \verbatim
29*>
30*>      CGET51  generally checks a decomposition of the form
31*>
32*>              A = U B V**H
33*>
34*>      where **H means conjugate transpose and U and V are unitary.
35*>
36*>      Specifically, if ITYPE=1
37*>
38*>              RESULT = | A - U B V**H | / ( |A| n ulp )
39*>
40*>      If ITYPE=2, then:
41*>
42*>              RESULT = | A - B | / ( |A| n ulp )
43*>
44*>      If ITYPE=3, then:
45*>
46*>              RESULT = | I - U U**H | / ( n ulp )
47*> \endverbatim
48*
49*  Arguments:
50*  ==========
51*
52*> \param[in] ITYPE
53*> \verbatim
54*>          ITYPE is INTEGER
55*>          Specifies the type of tests to be performed.
56*>          =1: RESULT = | A - U B V**H | / ( |A| n ulp )
57*>          =2: RESULT = | A - B | / ( |A| n ulp )
58*>          =3: RESULT = | I - U U**H | / ( n ulp )
59*> \endverbatim
60*>
61*> \param[in] N
62*> \verbatim
63*>          N is INTEGER
64*>          The size of the matrix.  If it is zero, CGET51 does nothing.
65*>          It must be at least zero.
66*> \endverbatim
67*>
68*> \param[in] A
69*> \verbatim
70*>          A is COMPLEX array, dimension (LDA, N)
71*>          The original (unfactored) matrix.
72*> \endverbatim
73*>
74*> \param[in] LDA
75*> \verbatim
76*>          LDA is INTEGER
77*>          The leading dimension of A.  It must be at least 1
78*>          and at least N.
79*> \endverbatim
80*>
81*> \param[in] B
82*> \verbatim
83*>          B is COMPLEX array, dimension (LDB, N)
84*>          The factored matrix.
85*> \endverbatim
86*>
87*> \param[in] LDB
88*> \verbatim
89*>          LDB is INTEGER
90*>          The leading dimension of B.  It must be at least 1
91*>          and at least N.
92*> \endverbatim
93*>
94*> \param[in] U
95*> \verbatim
96*>          U is COMPLEX array, dimension (LDU, N)
97*>          The unitary matrix on the left-hand side in the
98*>          decomposition.
99*>          Not referenced if ITYPE=2
100*> \endverbatim
101*>
102*> \param[in] LDU
103*> \verbatim
104*>          LDU is INTEGER
105*>          The leading dimension of U.  LDU must be at least N and
106*>          at least 1.
107*> \endverbatim
108*>
109*> \param[in] V
110*> \verbatim
111*>          V is COMPLEX array, dimension (LDV, N)
112*>          The unitary matrix on the left-hand side in the
113*>          decomposition.
114*>          Not referenced if ITYPE=2
115*> \endverbatim
116*>
117*> \param[in] LDV
118*> \verbatim
119*>          LDV is INTEGER
120*>          The leading dimension of V.  LDV must be at least N and
121*>          at least 1.
122*> \endverbatim
123*>
124*> \param[out] WORK
125*> \verbatim
126*>          WORK is COMPLEX array, dimension (2*N**2)
127*> \endverbatim
128*>
129*> \param[out] RWORK
130*> \verbatim
131*>          RWORK is REAL array, dimension (N)
132*> \endverbatim
133*>
134*> \param[out] RESULT
135*> \verbatim
136*>          RESULT is REAL
137*>          The values computed by the test specified by ITYPE.  The
138*>          value is currently limited to 1/ulp, to avoid overflow.
139*>          Errors are flagged by RESULT=10/ulp.
140*> \endverbatim
141*
142*  Authors:
143*  ========
144*
145*> \author Univ. of Tennessee
146*> \author Univ. of California Berkeley
147*> \author Univ. of Colorado Denver
148*> \author NAG Ltd.
149*
150*> \ingroup complex_eig
151*
152*  =====================================================================
153      SUBROUTINE CGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK,
154     $                   RWORK, RESULT )
155*
156*  -- LAPACK test routine --
157*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
158*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
159*
160*     .. Scalar Arguments ..
161      INTEGER            ITYPE, LDA, LDB, LDU, LDV, N
162      REAL               RESULT
163*     ..
164*     .. Array Arguments ..
165      REAL               RWORK( * )
166      COMPLEX            A( LDA, * ), B( LDB, * ), U( LDU, * ),
167     $                   V( LDV, * ), WORK( * )
168*     ..
169*
170*  =====================================================================
171*
172*     .. Parameters ..
173      REAL               ZERO, ONE, TEN
174      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 10.0E+0 )
175      COMPLEX            CZERO, CONE
176      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
177     $                   CONE = ( 1.0E+0, 0.0E+0 ) )
178*     ..
179*     .. Local Scalars ..
180      INTEGER            JCOL, JDIAG, JROW
181      REAL               ANORM, ULP, UNFL, WNORM
182*     ..
183*     .. External Functions ..
184      REAL               CLANGE, SLAMCH
185      EXTERNAL           CLANGE, SLAMCH
186*     ..
187*     .. External Subroutines ..
188      EXTERNAL           CGEMM, CLACPY
189*     ..
190*     .. Intrinsic Functions ..
191      INTRINSIC          MAX, MIN, REAL
192*     ..
193*     .. Executable Statements ..
194*
195      RESULT = ZERO
196      IF( N.LE.0 )
197     $   RETURN
198*
199*     Constants
200*
201      UNFL = SLAMCH( 'Safe minimum' )
202      ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
203*
204*     Some Error Checks
205*
206      IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
207         RESULT = TEN / ULP
208         RETURN
209      END IF
210*
211      IF( ITYPE.LE.2 ) THEN
212*
213*        Tests scaled by the norm(A)
214*
215         ANORM = MAX( CLANGE( '1', N, N, A, LDA, RWORK ), UNFL )
216*
217         IF( ITYPE.EQ.1 ) THEN
218*
219*           ITYPE=1: Compute W = A - U B V**H
220*
221            CALL CLACPY( ' ', N, N, A, LDA, WORK, N )
222            CALL CGEMM( 'N', 'N', N, N, N, CONE, U, LDU, B, LDB, CZERO,
223     $                  WORK( N**2+1 ), N )
224*
225            CALL CGEMM( 'N', 'C', N, N, N, -CONE, WORK( N**2+1 ), N, V,
226     $                  LDV, CONE, WORK, N )
227*
228         ELSE
229*
230*           ITYPE=2: Compute W = A - B
231*
232            CALL CLACPY( ' ', N, N, B, LDB, WORK, N )
233*
234            DO 20 JCOL = 1, N
235               DO 10 JROW = 1, N
236                  WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
237     $                - A( JROW, JCOL )
238   10          CONTINUE
239   20       CONTINUE
240         END IF
241*
242*        Compute norm(W)/ ( ulp*norm(A) )
243*
244         WNORM = CLANGE( '1', N, N, WORK, N, RWORK )
245*
246         IF( ANORM.GT.WNORM ) THEN
247            RESULT = ( WNORM / ANORM ) / ( N*ULP )
248         ELSE
249            IF( ANORM.LT.ONE ) THEN
250               RESULT = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
251            ELSE
252               RESULT = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP )
253            END IF
254         END IF
255*
256      ELSE
257*
258*        Tests not scaled by norm(A)
259*
260*        ITYPE=3: Compute  U U**H - I
261*
262         CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,
263     $               WORK, N )
264*
265         DO 30 JDIAG = 1, N
266            WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+
267     $         1 ) - CONE
268   30    CONTINUE
269*
270         RESULT = MIN( CLANGE( '1', N, N, WORK, N, RWORK ),
271     $            REAL( N ) ) / ( N*ULP )
272      END IF
273*
274      RETURN
275*
276*     End of CGET51
277*
278      END
279