LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ ccsdts()

subroutine ccsdts ( integer  M,
integer  P,
integer  Q,
complex, dimension( ldx, * )  X,
complex, dimension( ldx, * )  XF,
integer  LDX,
complex, dimension( ldu1, * )  U1,
integer  LDU1,
complex, dimension( ldu2, * )  U2,
integer  LDU2,
complex, dimension( ldv1t, * )  V1T,
integer  LDV1T,
complex, dimension( ldv2t, * )  V2T,
integer  LDV2T,
real, dimension( * )  THETA,
integer, dimension( * )  IWORK,
complex, dimension( lwork )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
real, dimension( 15 )  RESULT 
)

CCSDTS

Purpose:
 CCSDTS tests CUNCSD, 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 array, dimension (LDX,M)
          The M-by-M matrix X.
[out]XF
          XF is COMPLEX array, dimension (LDX,M)
          Details of the CSD of X, as returned by CUNCSD;
          see CUNCSD 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 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 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 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 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 REAL 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 CUNCSD for
          details.
[out]IWORK
          IWORK is INTEGER array, dimension (M)
[out]WORK
          WORK is COMPLEX array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK
[out]RWORK
          RWORK is REAL array
[out]RESULT
          RESULT is REAL 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 ). )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 226 of file ccsdts.f.

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