LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zgsvts3 ( integer  M,
integer  P,
integer  N,
complex*16, dimension( lda, * )  A,
complex*16, dimension( lda, * )  AF,
integer  LDA,
complex*16, dimension( ldb, * )  B,
complex*16, dimension( ldb, * )  BF,
integer  LDB,
complex*16, dimension( ldu, * )  U,
integer  LDU,
complex*16, dimension( ldv, * )  V,
integer  LDV,
complex*16, dimension( ldq, * )  Q,
integer  LDQ,
double precision, dimension( * )  ALPHA,
double precision, dimension( * )  BETA,
complex*16, dimension( ldr, * )  R,
integer  LDR,
integer, dimension( * )  IWORK,
complex*16, dimension( lwork )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
double precision, dimension( 6 )  RESULT 
)

ZGSVTS3

Purpose: 
 ZGSVTS3 tests ZGGSVD3, 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*16 array, dimension (LDA,M)
          The M-by-N matrix A.
[out]AF
          AF is COMPLEX*16 array, dimension (LDA,N)
          Details of the GSVD of A and B, as returned by ZGGSVD3,
          see ZGGSVD3 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*16 array, dimension (LDB,P)
          On entry, the P-by-N matrix B.
[out]BF
          BF is COMPLEX*16 array, dimension (LDB,N)
          Details of the GSVD of A and B, as returned by ZGGSVD3,
          see ZGGSVD3 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*16 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*16 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*16 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 DOUBLE PRECISION array, dimension (N)
[out]BETA
          BETA is DOUBLE PRECISION 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 ZGGSVD3 for details.
[out]R
          R is COMPLEX*16 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*16 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 DOUBLE PRECISION array, dimension (max(M,P,N))
[out]RESULT
          RESULT is DOUBLE PRECISION 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 zgsvts3.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  DOUBLE PRECISION alpha( * ), beta( * ), result( 6 ), rwork( * )
223  COMPLEX*16 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  DOUBLE PRECISION zero, one
232  parameter ( zero = 0.0d+0, one = 1.0d+0 )
233  COMPLEX*16 czero, cone
234  parameter ( czero = ( 0.0d+0, 0.0d+0 ),
235  $ cone = ( 1.0d+0, 0.0d+0 ) )
236 * ..
237 * .. Local Scalars ..
238  INTEGER i, info, j, k, l
239  DOUBLE PRECISION anorm, bnorm, resid, temp, ulp, ulpinv, unfl
240 * ..
241 * .. External Functions ..
242  DOUBLE PRECISION dlamch, zlange, zlanhe
243  EXTERNAL dlamch, zlange, zlanhe
244 * ..
245 * .. External Subroutines ..
246  EXTERNAL dcopy, zgemm, zggsvd3, zherk, zlacpy, zlaset
247 * ..
248 * .. Intrinsic Functions ..
249  INTRINSIC dble, max, min
250 * ..
251 * .. Executable Statements ..
252 *
253  ulp = dlamch( 'Precision' )
254  ulpinv = one / ulp
255  unfl = dlamch( 'Safe minimum' )
256 *
257 * Copy the matrix A to the array AF.
258 *
259  CALL zlacpy( 'Full', m, n, a, lda, af, lda )
260  CALL zlacpy( 'Full', p, n, b, ldb, bf, ldb )
261 *
262  anorm = max( zlange( '1', m, n, a, lda, rwork ), unfl )
263  bnorm = max( zlange( '1', p, n, b, ldb, rwork ), unfl )
264 *
265 * Factorize the matrices A and B in the arrays AF and BF.
266 *
267  CALL zggsvd3( '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 zgemm( 'No transpose', 'No transpose', m, n, n, cone, a, lda,
290  $ q, ldq, czero, work, lda )
291 *
292  CALL zgemm( '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 = zlange( '1', m, n, a, lda, rwork )
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 zgemm( 'No transpose', 'No transpose', p, n, n, cone, b, ldb,
320  $ q, ldq, czero, work, ldb )
321 *
322  CALL zgemm( '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 = zlange( '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 zlaset( 'Full', m, m, czero, cone, work, ldq )
344  CALL zherk( 'Upper', 'Conjugate transpose', m, m, -one, u, ldu,
345  $ one, work, ldu )
346 *
347 * Compute norm( I - U'*U ) / ( M * ULP ) .
348 *
349  resid = zlanhe( '1', 'Upper', m, work, ldu, rwork )
350  result( 3 ) = ( resid / dble( max( 1, m ) ) ) / ulp
351 *
352 * Compute I - V'*V
353 *
354  CALL zlaset( 'Full', p, p, czero, cone, work, ldv )
355  CALL zherk( 'Upper', 'Conjugate transpose', p, p, -one, v, ldv,
356  $ one, work, ldv )
357 *
358 * Compute norm( I - V'*V ) / ( P * ULP ) .
359 *
360  resid = zlanhe( '1', 'Upper', p, work, ldv, rwork )
361  result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp
362 *
363 * Compute I - Q'*Q
364 *
365  CALL zlaset( 'Full', n, n, czero, cone, work, ldq )
366  CALL zherk( 'Upper', 'Conjugate transpose', n, n, -one, q, ldq,
367  $ one, work, ldq )
368 *
369 * Compute norm( I - Q'*Q ) / ( N * ULP ) .
370 *
371  resid = zlanhe( '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, 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 ZGSVTS3
395 *
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
zggsvd3
subroutine zggsvd3(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)
ZGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
Definition: zggsvd3.f:355
zlanhe
double precision function zlanhe(NORM, UPLO, N, A, LDA, WORK)
ZLANHE 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: zlanhe.f:126
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
zlange
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
zherk
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:175

Here is the call graph for this function:

Here is the caller graph for this function: