LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
dgsvts3.f
Go to the documentation of this file.
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.6.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
dlaset
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
dggsvd3
subroutine dggsvd3(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, IWORK, INFO)
DGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
Definition: dggsvd3.f:351
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
dgsvts3
subroutine dgsvts3(M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, LWORK, RWORK, RESULT)
DGSVTS3
Definition: dgsvts3.f:212
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
dsyrk
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
Definition: dsyrk.f:171