LAPACK  3.9.0 LAPACK: Linear Algebra PACKage
cgsvts3.f
Go to the documentation of this file.
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 *> \date August 2015
204 *
205 *> \ingroup complex_eig
206 *
207 * =====================================================================
208  SUBROUTINE cgsvts3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
209  \$ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
210  \$ LWORK, RWORK, RESULT )
211 *
212 * -- LAPACK test routine (version 3.7.0) --
213 * -- LAPACK is a software package provided by Univ. of Tennessee, --
214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215 * August 2015
216 *
217 * .. Scalar Arguments ..
218  INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
219 * ..
220 * .. Array Arguments ..
221  INTEGER IWORK( * )
222  REAL ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
223  COMPLEX A( LDA, * ), AF( LDA, * ), B( LDB, * ),
224  \$ bf( ldb, * ), q( ldq, * ), r( ldr, * ),
225  \$ u( ldu, * ), v( ldv, * ), work( lwork )
226 * ..
227 *
228 * =====================================================================
229 *
230 * .. Parameters ..
231  REAL ZERO, ONE
232  PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
233  COMPLEX CZERO, CONE
234  parameter( czero = ( 0.0e+0, 0.0e+0 ),
235  \$ cone = ( 1.0e+0, 0.0e+0 ) )
236 * ..
237 * .. Local Scalars ..
238  INTEGER I, INFO, J, K, L
239  REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
240 * ..
241 * .. External Functions ..
242  REAL CLANGE, CLANHE, SLAMCH
243  EXTERNAL CLANGE, CLANHE, SLAMCH
244 * ..
245 * .. External Subroutines ..
246  EXTERNAL cgemm, cggsvd3, cherk, clacpy, claset, scopy
247 * ..
248 * .. Intrinsic Functions ..
249  INTRINSIC max, min, real
250 * ..
251 * .. Executable Statements ..
252 *
253  ulp = slamch( 'Precision' )
254  ulpinv = one / ulp
255  unfl = slamch( 'Safe minimum' )
256 *
257 * Copy the matrix A to the array AF.
258 *
259  CALL clacpy( 'Full', m, n, a, lda, af, lda )
260  CALL clacpy( 'Full', p, n, b, ldb, bf, ldb )
261 *
262  anorm = max( clange( '1', m, n, a, lda, rwork ), unfl )
263  bnorm = max( clange( '1', p, n, b, ldb, rwork ), unfl )
264 *
265 * Factorize the matrices A and B in the arrays AF and BF.
266 *
267  CALL cggsvd3( 'U', 'V', 'Q', m, n, p, k, l, af, lda, bf, ldb,
268  \$ alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork,
269  \$ rwork, iwork, info )
270 *
271 * Copy R
272 *
273  DO 20 i = 1, min( k+l, m )
274  DO 10 j = i, k + l
275  r( i, j ) = af( i, n-k-l+j )
276  10 CONTINUE
277  20 CONTINUE
278 *
279  IF( m-k-l.LT.0 ) THEN
280  DO 40 i = m + 1, k + l
281  DO 30 j = i, k + l
282  r( i, j ) = bf( i-k, n-k-l+j )
283  30 CONTINUE
284  40 CONTINUE
285  END IF
286 *
287 * Compute A:= U'*A*Q - D1*R
288 *
289  CALL cgemm( 'No transpose', 'No transpose', m, n, n, cone, a, lda,
290  \$ q, ldq, czero, work, lda )
291 *
292  CALL cgemm( 'Conjugate transpose', 'No transpose', m, n, m, cone,
293  \$ u, ldu, work, lda, czero, a, lda )
294 *
295  DO 60 i = 1, k
296  DO 50 j = i, k + l
297  a( i, n-k-l+j ) = a( i, n-k-l+j ) - r( i, j )
298  50 CONTINUE
299  60 CONTINUE
300 *
301  DO 80 i = k + 1, min( k+l, m )
302  DO 70 j = i, k + l
303  a( i, n-k-l+j ) = a( i, n-k-l+j ) - alpha( i )*r( i, j )
304  70 CONTINUE
305  80 CONTINUE
306 *
307 * Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
308 *
309  resid = clange( '1', m, n, a, lda, rwork )
310  IF( anorm.GT.zero ) THEN
311  result( 1 ) = ( ( resid / real( 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 cgemm( 'No transpose', 'No transpose', p, n, n, cone, b, ldb,
320  \$ q, ldq, czero, work, ldb )
321 *
322  CALL cgemm( 'Conjugate transpose', 'No transpose', p, n, p, cone,
323  \$ v, ldv, work, ldb, czero, 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 = clange( '1', p, n, b, ldb, rwork )
334  IF( bnorm.GT.zero ) THEN
335  result( 2 ) = ( ( resid / real( 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 claset( 'Full', m, m, czero, cone, work, ldq )
344  CALL cherk( 'Upper', 'Conjugate transpose', m, m, -one, u, ldu,
345  \$ one, work, ldu )
346 *
347 * Compute norm( I - U'*U ) / ( M * ULP ) .
348 *
349  resid = clanhe( '1', 'Upper', m, work, ldu, rwork )
350  result( 3 ) = ( resid / real( max( 1, m ) ) ) / ulp
351 *
352 * Compute I - V'*V
353 *
354  CALL claset( 'Full', p, p, czero, cone, work, ldv )
355  CALL cherk( 'Upper', 'Conjugate transpose', p, p, -one, v, ldv,
356  \$ one, work, ldv )
357 *
358 * Compute norm( I - V'*V ) / ( P * ULP ) .
359 *
360  resid = clanhe( '1', 'Upper', p, work, ldv, rwork )
361  result( 4 ) = ( resid / real( max( 1, p ) ) ) / ulp
362 *
363 * Compute I - Q'*Q
364 *
365  CALL claset( 'Full', n, n, czero, cone, work, ldq )
366  CALL cherk( 'Upper', 'Conjugate transpose', n, n, -one, q, ldq,
367  \$ one, work, ldq )
368 *
369 * Compute norm( I - Q'*Q ) / ( N * ULP ) .
370 *
371  resid = clanhe( '1', 'Upper', n, work, ldq, rwork )
372  result( 5 ) = ( resid / real( max( 1, n ) ) ) / ulp
373 *
374 * Check sorting
375 *
376  CALL scopy( n, alpha, 1, rwork, 1 )
377  DO 110 i = k + 1, min( k+l, m )
378  j = iwork( i )
379  IF( i.NE.j ) THEN
380  temp = rwork( i )
381  rwork( i ) = rwork( j )
382  rwork( 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( rwork( i ).LT.rwork( i+1 ) )
389  \$ result( 6 ) = ulpinv
390  120 CONTINUE
391 *
392  RETURN
393 *
394 * End of CGSVTS3
395 *
396  END
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
cggsvd3
subroutine cggsvd3(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, RWORK, IWORK, INFO)
CGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
Definition: cggsvd3.f:356
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
cherk
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:175
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
cgsvts3
subroutine cgsvts3(M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, LWORK, RWORK, RESULT)
CGSVTS3
Definition: cgsvts3.f:211
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:108