1*> \brief \b SGSVTS3
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 SGSVTS3( 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               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*> SGSVTS3 tests SGGSVD3, 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 REAL array, dimension (LDA,M)
62*>          The M-by-N matrix A.
63*> \endverbatim
64*>
65*> \param[out] AF
66*> \verbatim
67*>          AF is REAL array, dimension (LDA,N)
68*>          Details of the GSVD of A and B, as returned by SGGSVD3,
69*>          see SGGSVD3 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 REAL 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 REAL array, dimension (LDB,N)
88*>          Details of the GSVD of A and B, as returned by SGGSVD3,
89*>          see SGGSVD3 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (N)
138*> \endverbatim
139*>
140*> \param[out] BETA
141*> \verbatim
142*>          BETA is REAL 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 SGGSVD3 for details.
147*> \endverbatim
148*>
149*> \param[out] R
150*> \verbatim
151*>          R is REAL 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 REAL 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 REAL array, dimension (max(M,P,N))
181*> \endverbatim
182*>
183*> \param[out] RESULT
184*> \verbatim
185*>          RESULT is REAL 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*> \ingroup single_eig
205*
206*  =====================================================================
207      SUBROUTINE SGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
208     $                    LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
209     $                    LWORK, RWORK, RESULT )
210*
211*  -- LAPACK test routine --
212*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
213*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
214*
215*     .. Scalar Arguments ..
216      INTEGER            LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
217*     ..
218*     .. Array Arguments ..
219      INTEGER            IWORK( * )
220      REAL               A( LDA, * ), AF( LDA, * ), ALPHA( * ),
221     $                   B( LDB, * ), BETA( * ), BF( LDB, * ),
222     $                   Q( LDQ, * ), R( LDR, * ), RESULT( 6 ),
223     $                   RWORK( * ), U( LDU, * ), V( LDV, * ),
224     $                   WORK( LWORK )
225*     ..
226*
227*  =====================================================================
228*
229*     .. Parameters ..
230      REAL               ZERO, ONE
231      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
232*     ..
233*     .. Local Scalars ..
234      INTEGER            I, INFO, J, K, L
235      REAL               ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
236*     ..
237*     .. External Functions ..
238      REAL               SLAMCH, SLANGE, SLANSY
239      EXTERNAL           SLAMCH, SLANGE, SLANSY
240*     ..
241*     .. External Subroutines ..
242      EXTERNAL           SCOPY, SGEMM, SGGSVD3, SLACPY, SLASET, SSYRK
243*     ..
244*     .. Intrinsic Functions ..
245      INTRINSIC          MAX, MIN, REAL
246*     ..
247*     .. Executable Statements ..
248*
249      ULP = SLAMCH( 'Precision' )
250      ULPINV = ONE / ULP
251      UNFL = SLAMCH( 'Safe minimum' )
252*
253*     Copy the matrix A to the array AF.
254*
255      CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA )
256      CALL SLACPY( 'Full', P, N, B, LDB, BF, LDB )
257*
258      ANORM = MAX( SLANGE( '1', M, N, A, LDA, RWORK ), UNFL )
259      BNORM = MAX( SLANGE( '1', P, N, B, LDB, RWORK ), UNFL )
260*
261*     Factorize the matrices A and B in the arrays AF and BF.
262*
263      CALL SGGSVD3( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB,
264     $              ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK,
265     $              IWORK, INFO )
266*
267*     Copy R
268*
269      DO 20 I = 1, MIN( K+L, M )
270         DO 10 J = I, K + L
271            R( I, J ) = AF( I, N-K-L+J )
272   10    CONTINUE
273   20 CONTINUE
274*
275      IF( M-K-L.LT.0 ) THEN
276         DO 40 I = M + 1, K + L
277            DO 30 J = I, K + L
278               R( I, J ) = BF( I-K, N-K-L+J )
279   30       CONTINUE
280   40    CONTINUE
281      END IF
282*
283*     Compute A:= U'*A*Q - D1*R
284*
285      CALL SGEMM( 'No transpose', 'No transpose', M, N, N, ONE, A, LDA,
286     $            Q, LDQ, ZERO, WORK, LDA )
287*
288      CALL SGEMM( 'Transpose', 'No transpose', M, N, M, ONE, U, LDU,
289     $            WORK, LDA, ZERO, A, LDA )
290*
291      DO 60 I = 1, K
292         DO 50 J = I, K + L
293            A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J )
294   50    CONTINUE
295   60 CONTINUE
296*
297      DO 80 I = K + 1, MIN( K+L, M )
298         DO 70 J = I, K + L
299            A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J )
300   70    CONTINUE
301   80 CONTINUE
302*
303*     Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
304*
305      RESID = SLANGE( '1', M, N, A, LDA, RWORK )
306*
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 SGEMM( 'No transpose', 'No transpose', P, N, N, ONE, B, LDB,
317     $            Q, LDQ, ZERO, WORK, LDB )
318*
319      CALL SGEMM( 'Transpose', 'No transpose', P, N, P, ONE, V, LDV,
320     $            WORK, LDB, ZERO, 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 = SLANGE( '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 SLASET( 'Full', M, M, ZERO, ONE, WORK, LDQ )
341      CALL SSYRK( 'Upper', 'Transpose', M, M, -ONE, U, LDU, ONE, WORK,
342     $            LDU )
343*
344*     Compute norm( I - U'*U ) / ( M * ULP ) .
345*
346      RESID = SLANSY( '1', 'Upper', M, WORK, LDU, RWORK )
347      RESULT( 3 ) = ( RESID / REAL( MAX( 1, M ) ) ) / ULP
348*
349*     Compute I - V'*V
350*
351      CALL SLASET( 'Full', P, P, ZERO, ONE, WORK, LDV )
352      CALL SSYRK( 'Upper', 'Transpose', P, P, -ONE, V, LDV, ONE, WORK,
353     $            LDV )
354*
355*     Compute norm( I - V'*V ) / ( P * ULP ) .
356*
357      RESID = SLANSY( '1', 'Upper', P, WORK, LDV, RWORK )
358      RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP
359*
360*     Compute I - Q'*Q
361*
362      CALL SLASET( 'Full', N, N, ZERO, ONE, WORK, LDQ )
363      CALL SSYRK( 'Upper', 'Transpose', N, N, -ONE, Q, LDQ, ONE, WORK,
364     $            LDQ )
365*
366*     Compute norm( I - Q'*Q ) / ( N * ULP ) .
367*
368      RESID = SLANSY( '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, WORK, 1 )
374      DO 110 I = K + 1, MIN( K+L, M )
375         J = IWORK( I )
376         IF( I.NE.J ) THEN
377            TEMP = WORK( I )
378            WORK( I ) = WORK( J )
379            WORK( 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( WORK( I ).LT.WORK( I+1 ) )
386     $      RESULT( 6 ) = ULPINV
387  120 CONTINUE
388*
389      RETURN
390*
391*     End of SGSVTS3
392*
393      END
394