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

◆ 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*
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
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:103
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:115
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:124
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:106
cuncsd2by1
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
cuncsd
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
Here is the call graph for this function:
Here is the caller graph for this function: