1*> \brief \b CGSVTS3
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 CGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
12*                           LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
13*                           LWORK, RWORK, RESULT )
14*
15*       .. Scalar Arguments ..
16*       INTEGER            LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
17*       ..
18*       .. Array Arguments ..
19*       INTEGER            IWORK( * )
20*       REAL               ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
21*       COMPLEX            A( LDA, * ), AF( LDA, * ), B( LDB, * ),
22*      $                   BF( LDB, * ), Q( LDQ, * ), R( LDR, * ),
23*      $                   U( LDU, * ), V( LDV, * ), WORK( LWORK )
24*       ..
25*
26*
27*> \par Purpose:
28*  =============
29*>
30*> \verbatim
31*>
32*> CGSVTS3 tests CGGSVD3, which computes the GSVD of an M-by-N matrix A
33*> and a P-by-N matrix B:
34*>              U'*A*Q = D1*R and V'*B*Q = D2*R.
35*> \endverbatim
36*
37*  Arguments:
38*  ==========
39*
40*> \param[in] M
41*> \verbatim
42*>          M is INTEGER
43*>          The number of rows of the matrix A.  M >= 0.
44*> \endverbatim
45*>
46*> \param[in] P
47*> \verbatim
48*>          P is INTEGER
49*>          The number of rows of the matrix B.  P >= 0.
50*> \endverbatim
51*>
52*> \param[in] N
53*> \verbatim
54*>          N is INTEGER
55*>          The number of columns of the matrices A and B.  N >= 0.
56*> \endverbatim
57*>
58*> \param[in] A
59*> \verbatim
60*>          A is COMPLEX array, dimension (LDA,M)
61*>          The M-by-N matrix A.
62*> \endverbatim
63*>
64*> \param[out] AF
65*> \verbatim
66*>          AF is COMPLEX array, dimension (LDA,N)
67*>          Details of the GSVD of A and B, as returned by CGGSVD3,
68*>          see CGGSVD3 for further details.
69*> \endverbatim
70*>
71*> \param[in] LDA
72*> \verbatim
73*>          LDA is INTEGER
74*>          The leading dimension of the arrays A and AF.
75*>          LDA >= max( 1,M ).
76*> \endverbatim
77*>
78*> \param[in] B
79*> \verbatim
80*>          B is COMPLEX array, dimension (LDB,P)
81*>          On entry, the P-by-N matrix B.
82*> \endverbatim
83*>
84*> \param[out] BF
85*> \verbatim
86*>          BF is COMPLEX array, dimension (LDB,N)
87*>          Details of the GSVD of A and B, as returned by CGGSVD3,
88*>          see CGGSVD3 for further details.
89*> \endverbatim
90*>
91*> \param[in] LDB
92*> \verbatim
93*>          LDB is INTEGER
94*>          The leading dimension of the arrays B and BF.
95*>          LDB >= max(1,P).
96*> \endverbatim
97*>
98*> \param[out] U
99*> \verbatim
100*>          U is COMPLEX array, dimension(LDU,M)
101*>          The M by M unitary matrix U.
102*> \endverbatim
103*>
104*> \param[in] LDU
105*> \verbatim
106*>          LDU is INTEGER
107*>          The leading dimension of the array U. LDU >= max(1,M).
108*> \endverbatim
109*>
110*> \param[out] V
111*> \verbatim
112*>          V is COMPLEX array, dimension(LDV,M)
113*>          The P by P unitary matrix V.
114*> \endverbatim
115*>
116*> \param[in] LDV
117*> \verbatim
118*>          LDV is INTEGER
119*>          The leading dimension of the array V. LDV >= max(1,P).
120*> \endverbatim
121*>
122*> \param[out] Q
123*> \verbatim
124*>          Q is COMPLEX array, dimension(LDQ,N)
125*>          The N by N unitary matrix Q.
126*> \endverbatim
127*>
128*> \param[in] LDQ
129*> \verbatim
130*>          LDQ is INTEGER
131*>          The leading dimension of the array Q. LDQ >= max(1,N).
132*> \endverbatim
133*>
134*> \param[out] ALPHA
135*> \verbatim
136*>          ALPHA is REAL array, dimension (N)
137*> \endverbatim
138*>
139*> \param[out] BETA
140*> \verbatim
141*>          BETA is REAL array, dimension (N)
142*>
143*>          The generalized singular value pairs of A and B, the
144*>          ``diagonal'' matrices D1 and D2 are constructed from
145*>          ALPHA and BETA, see subroutine CGGSVD3 for details.
146*> \endverbatim
147*>
148*> \param[out] R
149*> \verbatim
150*>          R is COMPLEX array, dimension(LDQ,N)
151*>          The upper triangular matrix R.
152*> \endverbatim
153*>
154*> \param[in] LDR
155*> \verbatim
156*>          LDR is INTEGER
157*>          The leading dimension of the array R. LDR >= max(1,N).
158*> \endverbatim
159*>
160*> \param[out] IWORK
161*> \verbatim
162*>          IWORK is INTEGER array, dimension (N)
163*> \endverbatim
164*>
165*> \param[out] WORK
166*> \verbatim
167*>          WORK is COMPLEX array, dimension (LWORK)
168*> \endverbatim
169*>
170*> \param[in] LWORK
171*> \verbatim
172*>          LWORK is INTEGER
173*>          The dimension of the array WORK,
174*>          LWORK >= max(M,P,N)*max(M,P,N).
175*> \endverbatim
176*>
177*> \param[out] RWORK
178*> \verbatim
179*>          RWORK is REAL array, dimension (max(M,P,N))
180*> \endverbatim
181*>
182*> \param[out] RESULT
183*> \verbatim
184*>          RESULT is REAL array, dimension (6)
185*>          The test ratios:
186*>          RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP)
187*>          RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP)
188*>          RESULT(3) = norm( I - U'*U ) / ( M*ULP )
189*>          RESULT(4) = norm( I - V'*V ) / ( P*ULP )
190*>          RESULT(5) = norm( I - Q'*Q ) / ( N*ULP )
191*>          RESULT(6) = 0        if ALPHA is in decreasing order;
192*>                    = ULPINV   otherwise.
193*> \endverbatim
194*
195*  Authors:
196*  ========
197*
198*> \author Univ. of Tennessee
199*> \author Univ. of California Berkeley
200*> \author Univ. of Colorado Denver
201*> \author NAG Ltd.
202*
203*> \ingroup complex_eig
204*
205*  =====================================================================
206      SUBROUTINE CGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
207     $                    LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
208     $                    LWORK, RWORK, RESULT )
209*
210*  -- LAPACK test routine --
211*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
212*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213*
214*     .. Scalar Arguments ..
215      INTEGER            LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
216*     ..
217*     .. Array Arguments ..
218      INTEGER            IWORK( * )
219      REAL               ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
220      COMPLEX            A( LDA, * ), AF( LDA, * ), B( LDB, * ),
221     $                   BF( LDB, * ), Q( LDQ, * ), R( LDR, * ),
222     $                   U( LDU, * ), V( LDV, * ), WORK( LWORK )
223*     ..
224*
225*  =====================================================================
226*
227*     .. Parameters ..
228      REAL               ZERO, ONE
229      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
230      COMPLEX            CZERO, CONE
231      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
232     $                   CONE = ( 1.0E+0, 0.0E+0 ) )
233*     ..
234*     .. Local Scalars ..
235      INTEGER            I, INFO, J, K, L
236      REAL               ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
237*     ..
238*     .. External Functions ..
239      REAL               CLANGE, CLANHE, SLAMCH
240      EXTERNAL           CLANGE, CLANHE, SLAMCH
241*     ..
242*     .. External Subroutines ..
243      EXTERNAL           CGEMM, CGGSVD3, CHERK, CLACPY, CLASET, SCOPY
244*     ..
245*     .. Intrinsic Functions ..
246      INTRINSIC          MAX, MIN, REAL
247*     ..
248*     .. Executable Statements ..
249*
250      ULP = SLAMCH( 'Precision' )
251      ULPINV = ONE / ULP
252      UNFL = SLAMCH( 'Safe minimum' )
253*
254*     Copy the matrix A to the array AF.
255*
256      CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA )
257      CALL CLACPY( 'Full', P, N, B, LDB, BF, LDB )
258*
259      ANORM = MAX( CLANGE( '1', M, N, A, LDA, RWORK ), UNFL )
260      BNORM = MAX( CLANGE( '1', P, N, B, LDB, RWORK ), UNFL )
261*
262*     Factorize the matrices A and B in the arrays AF and BF.
263*
264      CALL CGGSVD3( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB,
265     $              ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK,
266     $              RWORK, IWORK, INFO )
267*
268*     Copy R
269*
270      DO 20 I = 1, MIN( K+L, M )
271         DO 10 J = I, K + L
272            R( I, J ) = AF( I, N-K-L+J )
273   10    CONTINUE
274   20 CONTINUE
275*
276      IF( M-K-L.LT.0 ) THEN
277         DO 40 I = M + 1, K + L
278            DO 30 J = I, K + L
279               R( I, J ) = BF( I-K, N-K-L+J )
280   30       CONTINUE
281   40    CONTINUE
282      END IF
283*
284*     Compute A:= U'*A*Q - D1*R
285*
286      CALL CGEMM( 'No transpose', 'No transpose', M, N, N, CONE, A, LDA,
287     $            Q, LDQ, CZERO, WORK, LDA )
288*
289      CALL CGEMM( 'Conjugate transpose', 'No transpose', M, N, M, CONE,
290     $            U, LDU, WORK, LDA, CZERO, A, LDA )
291*
292      DO 60 I = 1, K
293         DO 50 J = I, K + L
294            A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J )
295   50    CONTINUE
296   60 CONTINUE
297*
298      DO 80 I = K + 1, MIN( K+L, M )
299         DO 70 J = I, K + L
300            A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J )
301   70    CONTINUE
302   80 CONTINUE
303*
304*     Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
305*
306      RESID = CLANGE( '1', M, N, A, LDA, RWORK )
307      IF( ANORM.GT.ZERO ) THEN
308         RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M, N ) ) ) / ANORM ) /
309     $                 ULP
310      ELSE
311         RESULT( 1 ) = ZERO
312      END IF
313*
314*     Compute B := V'*B*Q - D2*R
315*
316      CALL CGEMM( 'No transpose', 'No transpose', P, N, N, CONE, B, LDB,
317     $            Q, LDQ, CZERO, WORK, LDB )
318*
319      CALL CGEMM( 'Conjugate transpose', 'No transpose', P, N, P, CONE,
320     $            V, LDV, WORK, LDB, CZERO, B, LDB )
321*
322      DO 100 I = 1, L
323         DO 90 J = I, L
324            B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J )
325   90    CONTINUE
326  100 CONTINUE
327*
328*     Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) .
329*
330      RESID = CLANGE( '1', P, N, B, LDB, RWORK )
331      IF( BNORM.GT.ZERO ) THEN
332         RESULT( 2 ) = ( ( RESID / REAL( MAX( 1, P, N ) ) ) / BNORM ) /
333     $                 ULP
334      ELSE
335         RESULT( 2 ) = ZERO
336      END IF
337*
338*     Compute I - U'*U
339*
340      CALL CLASET( 'Full', M, M, CZERO, CONE, WORK, LDQ )
341      CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, U, LDU,
342     $            ONE, WORK, LDU )
343*
344*     Compute norm( I - U'*U ) / ( M * ULP ) .
345*
346      RESID = CLANHE( '1', 'Upper', M, WORK, LDU, RWORK )
347      RESULT( 3 ) = ( RESID / REAL( MAX( 1, M ) ) ) / ULP
348*
349*     Compute I - V'*V
350*
351      CALL CLASET( 'Full', P, P, CZERO, CONE, WORK, LDV )
352      CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, V, LDV,
353     $            ONE, WORK, LDV )
354*
355*     Compute norm( I - V'*V ) / ( P * ULP ) .
356*
357      RESID = CLANHE( '1', 'Upper', P, WORK, LDV, RWORK )
358      RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP
359*
360*     Compute I - Q'*Q
361*
362      CALL CLASET( 'Full', N, N, CZERO, CONE, WORK, LDQ )
363      CALL CHERK( 'Upper', 'Conjugate transpose', N, N, -ONE, Q, LDQ,
364     $            ONE, WORK, LDQ )
365*
366*     Compute norm( I - Q'*Q ) / ( N * ULP ) .
367*
368      RESID = CLANHE( '1', 'Upper', N, WORK, LDQ, RWORK )
369      RESULT( 5 ) = ( RESID / REAL( MAX( 1, N ) ) ) / ULP
370*
371*     Check sorting
372*
373      CALL SCOPY( N, ALPHA, 1, RWORK, 1 )
374      DO 110 I = K + 1, MIN( K+L, M )
375         J = IWORK( I )
376         IF( I.NE.J ) THEN
377            TEMP = RWORK( I )
378            RWORK( I ) = RWORK( J )
379            RWORK( J ) = TEMP
380         END IF
381  110 CONTINUE
382*
383      RESULT( 6 ) = ZERO
384      DO 120 I = K + 1, MIN( K+L, M ) - 1
385         IF( RWORK( I ).LT.RWORK( I+1 ) )
386     $      RESULT( 6 ) = ULPINV
387  120 CONTINUE
388*
389      RETURN
390*
391*     End of CGSVTS3
392*
393      END
394