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