LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
August 2015

Definition at line 211 of file cgsvts3.f.

211 *
212 * -- LAPACK test routine (version 3.6.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 *
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, or the element of largest absolute value of a complex Hermitian matrix.
Definition: clanhe.f:126
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
cherk
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:175
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:117
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
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
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: