LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine zcsdts ( integer M, integer P, integer Q, complex*16, dimension( ldx, * ) X, complex*16, dimension( ldx, * ) XF, integer LDX, complex*16, dimension( ldu1, * ) U1, integer LDU1, complex*16, dimension( ldu2, * ) U2, integer LDU2, complex*16, dimension( ldv1t, * ) V1T, integer LDV1T, complex*16, dimension( ldv2t, * ) V2T, integer LDV2T, double precision, dimension( * ) THETA, integer, dimension( * ) IWORK, complex*16, dimension( lwork ) WORK, integer LWORK, double precision, dimension( * ) RWORK, double precision, dimension( 15 ) RESULT )

ZCSDTS

Purpose:
``` ZCSDTS tests ZUNCSD, which, given an M-by-M partitioned unitary
matrix X,
Q  M-Q
X = [ X11 X12 ] P   ,
[ X21 X22 ] M-P

computes the CSD

[ U1    ]**T * [ X11 X12 ] * [ V1    ]
[    U2 ]      [ X21 X22 ]   [    V2 ]

[  I  0  0 |  0  0  0 ]
[  0  C  0 |  0 -S  0 ]
[  0  0  0 |  0  0 -I ]
= [---------------------] = [ D11 D12 ] .
[  0  0  0 |  I  0  0 ]   [ D21 D22 ]
[  0  S  0 |  0  C  0 ]
[  0  0  I |  0  0  0 ]

and also SORCSD2BY1, which, given
Q
[ X11 ] P   ,
[ X21 ] M-P

computes the 2-by-1 CSD

[  I  0  0 ]
[  0  C  0 ]
[  0  0  0 ]
[ U1    ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
[    U2 ]      [ X21 ]        [  0  0  0 ]   [ D21 ]
[  0  S  0 ]
[  0  0  I ]```
Parameters
 [in] M ``` M is INTEGER The number of rows of the matrix X. M >= 0.``` [in] P ``` P is INTEGER The number of rows of the matrix X11. P >= 0.``` [in] Q ``` Q is INTEGER The number of columns of the matrix X11. Q >= 0.``` [in] X ``` X is COMPLEX*16 array, dimension (LDX,M) The M-by-M matrix X.``` [out] XF ``` XF is COMPLEX*16 array, dimension (LDX,M) Details of the CSD of X, as returned by ZUNCSD; see ZUNCSD for further details.``` [in] LDX ``` LDX is INTEGER The leading dimension of the arrays X and XF. LDX >= max( 1,M ).``` [out] U1 ``` U1 is COMPLEX*16 array, dimension(LDU1,P) The P-by-P unitary matrix U1.``` [in] LDU1 ``` LDU1 is INTEGER The leading dimension of the array U1. LDU >= max(1,P).``` [out] U2 ``` U2 is COMPLEX*16 array, dimension(LDU2,M-P) The (M-P)-by-(M-P) unitary matrix U2.``` [in] LDU2 ``` LDU2 is INTEGER The leading dimension of the array U2. LDU >= max(1,M-P).``` [out] V1T ``` V1T is COMPLEX*16 array, dimension(LDV1T,Q) The Q-by-Q unitary matrix V1T.``` [in] LDV1T ``` LDV1T is INTEGER The leading dimension of the array V1T. LDV1T >= max(1,Q).``` [out] V2T ``` V2T is COMPLEX*16 array, dimension(LDV2T,M-Q) The (M-Q)-by-(M-Q) unitary matrix V2T.``` [in] LDV2T ``` LDV2T is INTEGER The leading dimension of the array V2T. LDV2T >= max(1,M-Q).``` [out] THETA ``` THETA is DOUBLE PRECISION array, dimension MIN(P,M-P,Q,M-Q) The CS values of X; the essentially diagonal matrices C and S are constructed from THETA; see subroutine ZUNCSD for details.``` [out] IWORK ` IWORK is INTEGER array, dimension (M)` [out] WORK ` WORK is COMPLEX*16 array, dimension (LWORK)` [in] LWORK ``` LWORK is INTEGER The dimension of the array WORK``` [out] RWORK ` RWORK is DOUBLE PRECISION array` [out] RESULT ``` RESULT is DOUBLE PRECISION array, dimension (15) The test ratios: First, the 2-by-2 CSD: RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP ) RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP ) RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP ) RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP ) RESULT(9) = 0 if THETA is in increasing order and all angles are in [0,pi/2]; = ULPINV otherwise. Then, the 2-by-1 CSD: RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP ) RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP ) RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP ) RESULT(15) = 0 if THETA is in increasing order and all angles are in [0,pi/2]; = ULPINV otherwise. ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )```
Date
November 2015

Definition at line 231 of file zcsdts.f.

231 *
232 * -- LAPACK test routine (version 3.6.0) --
233 * -- LAPACK is a software package provided by Univ. of Tennessee, --
234 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
235 * November 2015
236 *
237 * .. Scalar Arguments ..
238  INTEGER ldx, ldu1, ldu2, ldv1t, ldv2t, lwork, m, p, q
239 * ..
240 * .. Array Arguments ..
241  INTEGER iwork( * )
242  DOUBLE PRECISION result( 15 ), rwork( * ), theta( * )
243  COMPLEX*16 u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
244  \$ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
245  \$ xf( ldx, * )
246 * ..
247 *
248 * =====================================================================
249 *
250 * .. Parameters ..
251  DOUBLE PRECISION piover2, realone, realzero
252  parameter ( piover2 = 1.57079632679489662d0,
253  \$ realone = 1.0d0, realzero = 0.0d0 )
254  COMPLEX*16 zero, one
255  parameter ( zero = (0.0d0,0.0d0), one = (1.0d0,0.0d0) )
256 * ..
257 * .. Local Scalars ..
258  INTEGER i, info, r
259  DOUBLE PRECISION eps2, resid, ulp, ulpinv
260 * ..
261 * .. External Functions ..
262  DOUBLE PRECISION dlamch, zlange, zlanhe
263  EXTERNAL dlamch, zlange, zlanhe
264 * ..
265 * .. External Subroutines ..
266  EXTERNAL zgemm, zherk, zlacpy, zlaset, zuncsd,
267  \$ zuncsd2by1
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC cos, dble, dcmplx, max, min, REAL, sin
271 * ..
272 * .. Executable Statements ..
273 *
274  ulp = dlamch( 'Precision' )
275  ulpinv = realone / ulp
276 *
277 * The first half of the routine checks the 2-by-2 CSD
278 *
279  CALL zlaset( 'Full', m, m, zero, one, work, ldx )
280  CALL zherk( 'Upper', 'Conjugate transpose', m, m, -realone,
281  \$ x, ldx, realone, work, ldx )
282  IF (m.GT.0) THEN
283  eps2 = max( ulp,
284  \$ zlange( '1', m, m, work, ldx, rwork ) / dble( m ) )
285  ELSE
286  eps2 = ulp
287  END IF
288  r = min( p, m-p, q, m-q )
289 *
290 * Copy the matrix X to the array XF.
291 *
292  CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
293 *
294 * Compute the CSD
295 *
296  CALL zuncsd( 'Y', 'Y', 'Y', 'Y', 'N', 'D', m, p, q, xf(1,1), ldx,
297  \$ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
298  \$ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
299  \$ work, lwork, rwork, 17*(r+2), iwork, info )
300 *
301 * Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
302 *
303  CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
304 *
305  CALL zgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
306  \$ xf, ldx, v1t, ldv1t, zero, work, ldx )
307 *
308  CALL zgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
309  \$ u1, ldu1, work, ldx, zero, xf, ldx )
310 *
311  DO i = 1, min(p,q)-r
312  xf(i,i) = xf(i,i) - one
313  END DO
314  DO i = 1, r
315  xf(min(p,q)-r+i,min(p,q)-r+i) =
316  \$ xf(min(p,q)-r+i,min(p,q)-r+i) - dcmplx( cos(theta(i)),
317  \$ 0.0d0 )
318  END DO
319 *
320  CALL zgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
321  \$ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
322 *
323  CALL zgemm( 'Conjugate transpose', 'No transpose', p, m-q, p,
324  \$ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
325 *
326  DO i = 1, min(p,m-q)-r
327  xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
328  END DO
329  DO i = 1, r
330  xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
331  \$ xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
332  \$ dcmplx( sin(theta(r-i+1)), 0.0d0 )
333  END DO
334 *
335  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
336  \$ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
337 *
338  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
339  \$ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
340 *
341  DO i = 1, min(m-p,q)-r
342  xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
343  END DO
344  DO i = 1, r
345  xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
346  \$ xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
347  \$ dcmplx( sin(theta(r-i+1)), 0.0d0 )
348  END DO
349 *
350  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, m-q, m-q,
351  \$ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
352 *
353  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, m-q, m-p,
354  \$ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
355 *
356  DO i = 1, min(m-p,m-q)-r
357  xf(p+i,q+i) = xf(p+i,q+i) - one
358  END DO
359  DO i = 1, r
360  xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
361  \$ xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
362  \$ dcmplx( cos(theta(i)), 0.0d0 )
363  END DO
364 *
365 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
366 *
367  resid = zlange( '1', p, q, xf, ldx, rwork )
368  result( 1 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
369 *
370 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
371 *
372  resid = zlange( '1', p, m-q, xf(1,q+1), ldx, rwork )
373  result( 2 ) = ( resid / REAL(MAX(1,P,M-Q)) ) / eps2
374 *
375 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
376 *
377  resid = zlange( '1', m-p, q, xf(p+1,1), ldx, rwork )
378  result( 3 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
379 *
380 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
381 *
382  resid = zlange( '1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
383  result( 4 ) = ( resid / REAL(MAX(1,M-P,M-Q)) ) / eps2
384 *
385 * Compute I - U1'*U1
386 *
387  CALL zlaset( 'Full', p, p, zero, one, work, ldu1 )
388  CALL zherk( 'Upper', 'Conjugate transpose', p, p, -realone,
389  \$ u1, ldu1, realone, work, ldu1 )
390 *
391 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
392 *
393  resid = zlanhe( '1', 'Upper', p, work, ldu1, rwork )
394  result( 5 ) = ( resid / REAL(MAX(1,P)) ) / ulp
395 *
396 * Compute I - U2'*U2
397 *
398  CALL zlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
399  CALL zherk( 'Upper', 'Conjugate transpose', m-p, m-p, -realone,
400  \$ u2, ldu2, realone, work, ldu2 )
401 *
402 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
403 *
404  resid = zlanhe( '1', 'Upper', m-p, work, ldu2, rwork )
405  result( 6 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
406 *
407 * Compute I - V1T*V1T'
408 *
409  CALL zlaset( 'Full', q, q, zero, one, work, ldv1t )
410  CALL zherk( 'Upper', 'No transpose', q, q, -realone,
411  \$ v1t, ldv1t, realone, work, ldv1t )
412 *
413 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
414 *
415  resid = zlanhe( '1', 'Upper', q, work, ldv1t, rwork )
416  result( 7 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
417 *
418 * Compute I - V2T*V2T'
419 *
420  CALL zlaset( 'Full', m-q, m-q, zero, one, work, ldv2t )
421  CALL zherk( 'Upper', 'No transpose', m-q, m-q, -realone,
422  \$ v2t, ldv2t, realone, work, ldv2t )
423 *
424 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
425 *
426  resid = zlanhe( '1', 'Upper', m-q, work, ldv2t, rwork )
427  result( 8 ) = ( resid / REAL(MAX(1,M-Q)) ) / ulp
428 *
429 * Check sorting
430 *
431  result( 9 ) = realzero
432  DO i = 1, r
433  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
434  result( 9 ) = ulpinv
435  END IF
436  IF( i.GT.1) THEN
437  IF ( theta(i).LT.theta(i-1) ) THEN
438  result( 9 ) = ulpinv
439  END IF
440  END IF
441  END DO
442 *
443 * The second half of the routine checks the 2-by-1 CSD
444 *
445  CALL zlaset( 'Full', q, q, zero, one, work, ldx )
446  CALL zherk( 'Upper', 'Conjugate transpose', q, m, -realone,
447  \$ x, ldx, realone, work, ldx )
448  IF (m.GT.0) THEN
449  eps2 = max( ulp,
450  \$ zlange( '1', q, q, work, ldx, rwork ) / dble( m ) )
451  ELSE
452  eps2 = ulp
453  END IF
454  r = min( p, m-p, q, m-q )
455 *
456 * Copy the matrix X to the array XF.
457 *
458  CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
459 *
460 * Compute the CSD
461 *
462  CALL zuncsd2by1( 'Y', 'Y', 'Y', m, p, q, xf(1,1), ldx, xf(p+1,1),
463  \$ ldx, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work,
464  \$ lwork, rwork, 17*(r+2), iwork, info )
465 *
466 * Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
467 *
468  CALL zgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
469  \$ x, ldx, v1t, ldv1t, zero, work, ldx )
470 *
471  CALL zgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
472  \$ u1, ldu1, work, ldx, zero, x, ldx )
473 *
474  DO i = 1, min(p,q)-r
475  x(i,i) = x(i,i) - one
476  END DO
477  DO i = 1, r
478  x(min(p,q)-r+i,min(p,q)-r+i) =
479  \$ x(min(p,q)-r+i,min(p,q)-r+i) - dcmplx( cos(theta(i)),
480  \$ 0.0d0 )
481  END DO
482 *
483  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
484  \$ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
485 *
486  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
487  \$ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
488 *
489  DO i = 1, min(m-p,q)-r
490  x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
491  END DO
492  DO i = 1, r
493  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
494  \$ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
495  \$ dcmplx( sin(theta(r-i+1)), 0.0d0 )
496  END DO
497 *
498 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
499 *
500  resid = zlange( '1', p, q, x, ldx, rwork )
501  result( 10 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
502 *
503 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
504 *
505  resid = zlange( '1', m-p, q, x(p+1,1), ldx, rwork )
506  result( 11 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
507 *
508 * Compute I - U1'*U1
509 *
510  CALL zlaset( 'Full', p, p, zero, one, work, ldu1 )
511  CALL zherk( 'Upper', 'Conjugate transpose', p, p, -realone,
512  \$ u1, ldu1, realone, work, ldu1 )
513 *
514 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
515 *
516  resid = zlanhe( '1', 'Upper', p, work, ldu1, rwork )
517  result( 12 ) = ( resid / REAL(MAX(1,P)) ) / ulp
518 *
519 * Compute I - U2'*U2
520 *
521  CALL zlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
522  CALL zherk( 'Upper', 'Conjugate transpose', m-p, m-p, -realone,
523  \$ u2, ldu2, realone, work, ldu2 )
524 *
525 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
526 *
527  resid = zlanhe( '1', 'Upper', m-p, work, ldu2, rwork )
528  result( 13 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
529 *
530 * Compute I - V1T*V1T'
531 *
532  CALL zlaset( 'Full', q, q, zero, one, work, ldv1t )
533  CALL zherk( 'Upper', 'No transpose', q, q, -realone,
534  \$ v1t, ldv1t, realone, work, ldv1t )
535 *
536 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
537 *
538  resid = zlanhe( '1', 'Upper', q, work, ldv1t, rwork )
539  result( 14 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
540 *
541 * Check sorting
542 *
543  result( 15 ) = realzero
544  DO i = 1, r
545  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
546  result( 15 ) = ulpinv
547  END IF
548  IF( i.GT.1) THEN
549  IF ( theta(i).LT.theta(i-1) ) THEN
550  result( 15 ) = ulpinv
551  END IF
552  END IF
553  END DO
554 *
555  RETURN
556 *
557 * End of ZCSDTS
558 *
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
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
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
recursive subroutine zuncsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
ZUNCSD
Definition: zuncsd.f:322
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
subroutine zuncsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
ZUNCSD2BY1
Definition: zuncsd2by1.f:255
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: