LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ cgsvts3()

subroutine cgsvts3 ( integer m,
integer p,
integer n,
complex, dimension( lda, * ) a,
complex, dimension( lda, * ) af,
integer lda,
complex, dimension( ldb, * ) b,
complex, dimension( ldb, * ) bf,
integer ldb,
complex, dimension( ldu, * ) u,
integer ldu,
complex, dimension( ldv, * ) v,
integer ldv,
complex, dimension( ldq, * ) q,
integer ldq,
real, dimension( * ) alpha,
real, dimension( * ) beta,
complex, dimension( ldr, * ) r,
integer ldr,
integer, dimension( * ) iwork,
complex, dimension( lwork ) work,
integer lwork,
real, dimension( * ) rwork,
real, dimension( 6 ) result )

CGSVTS3

Purpose:
!>
!> CGSVTS3 tests CGGSVD3, which computes the GSVD of an M-by-N matrix A
!> and a P-by-N matrix B:
!>              U'*A*Q = D1*R and V'*B*Q = D2*R.
!>
Parameters
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.
!>
[in]P
!>          P is INTEGER
!>          The number of rows of the matrix B.  P >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrices A and B.  N >= 0.
!>
[in]A
!>          A is COMPLEX array, dimension (LDA,M)
!>          The M-by-N matrix A.
!>
[out]AF
!>          AF is COMPLEX array, dimension (LDA,N)
!>          Details of the GSVD of A and B, as returned by CGGSVD3,
!>          see CGGSVD3 for further details.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the arrays A and AF.
!>          LDA >= max( 1,M ).
!>
[in]B
!>          B is COMPLEX array, dimension (LDB,P)
!>          On entry, the P-by-N matrix B.
!>
[out]BF
!>          BF is COMPLEX array, dimension (LDB,N)
!>          Details of the GSVD of A and B, as returned by CGGSVD3,
!>          see CGGSVD3 for further details.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the arrays B and BF.
!>          LDB >= max(1,P).
!>
[out]U
!>          U is COMPLEX array, dimension(LDU,M)
!>          The M by M unitary matrix U.
!>
[in]LDU
!>          LDU is INTEGER
!>          The leading dimension of the array U. LDU >= max(1,M).
!>
[out]V
!>          V is COMPLEX array, dimension(LDV,M)
!>          The P by P unitary matrix V.
!>
[in]LDV
!>          LDV is INTEGER
!>          The leading dimension of the array V. LDV >= max(1,P).
!>
[out]Q
!>          Q is COMPLEX array, dimension(LDQ,N)
!>          The N by N unitary matrix Q.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q. LDQ >= max(1,N).
!>
[out]ALPHA
!>          ALPHA is REAL array, dimension (N)
!>
[out]BETA
!>          BETA is REAL array, dimension (N)
!>
!>          The generalized singular value pairs of A and B, the
!>          ``diagonal'' matrices D1 and D2 are constructed from
!>          ALPHA and BETA, see subroutine CGGSVD3 for details.
!>
[out]R
!>          R is COMPLEX array, dimension(LDQ,N)
!>          The upper triangular matrix R.
!>
[in]LDR
!>          LDR is INTEGER
!>          The leading dimension of the array R. LDR >= max(1,N).
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (LWORK)
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK,
!>          LWORK >= max(M,P,N)*max(M,P,N).
!>
[out]RWORK
!>          RWORK is REAL array, dimension (max(M,P,N))
!>
[out]RESULT
!>          RESULT is REAL array, dimension (6)
!>          The test ratios:
!>          RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP)
!>          RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP)
!>          RESULT(3) = norm( I - U'*U ) / ( M*ULP )
!>          RESULT(4) = norm( I - V'*V ) / ( P*ULP )
!>          RESULT(5) = norm( I - Q'*Q ) / ( N*ULP )
!>          RESULT(6) = 0        if ALPHA is in decreasing order;
!>                    = ULPINV   otherwise.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 206 of file cgsvts3.f.

209*
210* -- LAPACK test routine --
211* -- LAPACK is a software package provided by Univ. of Tennessee, --
212* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213*
214* .. Scalar Arguments ..
215 INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
216* ..
217* .. Array Arguments ..
218 INTEGER IWORK( * )
219 REAL ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
220 COMPLEX A( LDA, * ), AF( LDA, * ), B( LDB, * ),
221 $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ),
222 $ U( LDU, * ), V( LDV, * ), WORK( LWORK )
223* ..
224*
225* =====================================================================
226*
227* .. Parameters ..
228 REAL ZERO, ONE
229 parameter( zero = 0.0e+0, one = 1.0e+0 )
230 COMPLEX CZERO, CONE
231 parameter( czero = ( 0.0e+0, 0.0e+0 ),
232 $ cone = ( 1.0e+0, 0.0e+0 ) )
233* ..
234* .. Local Scalars ..
235 INTEGER I, INFO, J, K, L
236 REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
237* ..
238* .. External Functions ..
239 REAL CLANGE, CLANHE, SLAMCH
240 EXTERNAL clange, clanhe, slamch
241* ..
242* .. External Subroutines ..
243 EXTERNAL cgemm, cggsvd3, cherk, clacpy, claset, scopy
244* ..
245* .. Intrinsic Functions ..
246 INTRINSIC max, min, real
247* ..
248* .. Executable Statements ..
249*
250 ulp = slamch( 'Precision' )
251 ulpinv = one / ulp
252 unfl = slamch( 'Safe minimum' )
253*
254* Copy the matrix A to the array AF.
255*
256 CALL clacpy( 'Full', m, n, a, lda, af, lda )
257 CALL clacpy( 'Full', p, n, b, ldb, bf, ldb )
258*
259 anorm = max( clange( '1', m, n, a, lda, rwork ), unfl )
260 bnorm = max( clange( '1', p, n, b, ldb, rwork ), unfl )
261*
262* Factorize the matrices A and B in the arrays AF and BF.
263*
264 CALL cggsvd3( 'U', 'V', 'Q', m, n, p, k, l, af, lda, bf, ldb,
265 $ alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork,
266 $ rwork, iwork, info )
267*
268* Copy R
269*
270 DO 20 i = 1, min( k+l, m )
271 DO 10 j = i, k + l
272 r( i, j ) = af( i, n-k-l+j )
273 10 CONTINUE
274 20 CONTINUE
275*
276 IF( m-k-l.LT.0 ) THEN
277 DO 40 i = m + 1, k + l
278 DO 30 j = i, k + l
279 r( i, j ) = bf( i-k, n-k-l+j )
280 30 CONTINUE
281 40 CONTINUE
282 END IF
283*
284* Compute A:= U'*A*Q - D1*R
285*
286 CALL cgemm( 'No transpose', 'No transpose', m, n, n, cone, a, lda,
287 $ q, ldq, czero, work, lda )
288*
289 CALL cgemm( 'Conjugate transpose', 'No transpose', m, n, m, cone,
290 $ u, ldu, work, lda, czero, a, lda )
291*
292 DO 60 i = 1, k
293 DO 50 j = i, k + l
294 a( i, n-k-l+j ) = a( i, n-k-l+j ) - r( i, j )
295 50 CONTINUE
296 60 CONTINUE
297*
298 DO 80 i = k + 1, min( k+l, m )
299 DO 70 j = i, k + l
300 a( i, n-k-l+j ) = a( i, n-k-l+j ) - alpha( i )*r( i, j )
301 70 CONTINUE
302 80 CONTINUE
303*
304* Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
305*
306 resid = clange( '1', m, n, a, lda, rwork )
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 cgemm( 'No transpose', 'No transpose', p, n, n, cone, b, ldb,
317 $ q, ldq, czero, work, ldb )
318*
319 CALL cgemm( 'Conjugate transpose', 'No transpose', p, n, p, cone,
320 $ v, ldv, work, ldb, czero, 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 = clange( '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 claset( 'Full', m, m, czero, cone, work, ldq )
341 CALL cherk( 'Upper', 'Conjugate transpose', m, m, -one, u, ldu,
342 $ one, work, ldu )
343*
344* Compute norm( I - U'*U ) / ( M * ULP ) .
345*
346 resid = clanhe( '1', 'Upper', m, work, ldu, rwork )
347 result( 3 ) = ( resid / real( max( 1, m ) ) ) / ulp
348*
349* Compute I - V'*V
350*
351 CALL claset( 'Full', p, p, czero, cone, work, ldv )
352 CALL cherk( 'Upper', 'Conjugate transpose', p, p, -one, v, ldv,
353 $ one, work, ldv )
354*
355* Compute norm( I - V'*V ) / ( P * ULP ) .
356*
357 resid = clanhe( '1', 'Upper', p, work, ldv, rwork )
358 result( 4 ) = ( resid / real( max( 1, p ) ) ) / ulp
359*
360* Compute I - Q'*Q
361*
362 CALL claset( 'Full', n, n, czero, cone, work, ldq )
363 CALL cherk( 'Upper', 'Conjugate transpose', n, n, -one, q, ldq,
364 $ one, work, ldq )
365*
366* Compute norm( I - Q'*Q ) / ( N * ULP ) .
367*
368 resid = clanhe( '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, rwork, 1 )
374 DO 110 i = k + 1, min( k+l, m )
375 j = iwork( i )
376 IF( i.NE.j ) THEN
377 temp = rwork( i )
378 rwork( i ) = rwork( j )
379 rwork( 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( rwork( i ).LT.rwork( i+1 ) )
386 $ result( 6 ) = ulpinv
387 120 CONTINUE
388*
389 RETURN
390*
391* End of CGSVTS3
392*
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
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:352
cherk
subroutine cherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
CHERK
Definition cherk.f:173
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:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clange
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:113
clanhe
real function clanhe(norm, uplo, n, a, lda, work)
CLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanhe.f:122
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:104
Here is the call graph for this function:
Here is the caller graph for this function: