LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dgsvts.f
Go to the documentation of this file.
1 *> \brief \b DGSVTS
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 DGSVTS( 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 *> DGSVTS tests DGGSVD, 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 DGGSVD,
69 *> see DGGSVD 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 DGGSVD,
89 *> see DGGSVD 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 DGGSVD 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 November 2011
205 *
206 *> \ingroup double_eig
207 *
208 * =====================================================================
209  SUBROUTINE dgsvts( 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.4.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 * November 2011
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, dggsvd, 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 dggsvd( 'U', 'V', 'Q', m, n, p, k, l, af, lda, bf, ldb,
267  $ alpha, beta, u, ldu, v, ldv, q, ldq, work, iwork,
268  $ 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 DGSVTS
395 *
396  END