228 SUBROUTINE ccsdts( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
229 $ ldv1t, v2t, ldv2t, theta, iwork, work, lwork,
238 INTEGER ldx, ldu1, ldu2, ldv1t, ldv2t, lwork, m, p, q
242 REAL result( 15 ), rwork( * ), theta( * )
243 COMPLEX u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
244 $ v2t( ldv2t, * ), work( lwork ),
x( ldx, * ),
251 REAL piover2, realone, realzero
252 parameter( piover2 = 1.57079632679489662e0,
253 $ realone = 1.0e0, realzero = 0.0e0 )
255 parameter( zero = (0.0e0,0.0e0), one = (1.0e0,0.0e0) )
259 REAL eps2, resid, ulp, ulpinv
269 INTRINSIC cmplx, cos, max, min,
REAL, sin
273 ulp =
slamch(
'Precision' )
274 ulpinv = realone / ulp
278 CALL
claset(
'Full', m, m, zero, one, work, ldx )
279 CALL
cherk(
'Upper',
'Conjugate transpose', m, m, -realone,
280 $
x, ldx, realone, work, ldx )
283 $
clange(
'1', m, m, work, ldx, rwork ) /
REAL( M ) )
287 r = min( p, m-p, q, m-q )
291 CALL
clacpy(
'Full', m, m,
x, ldx, xf, ldx )
295 CALL
cuncsd(
'Y',
'Y',
'Y',
'Y',
'N',
'D', m, p, q, xf(1,1), ldx,
296 $ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
297 $ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
298 $ work, lwork, rwork, 17*(r+2), iwork, info )
302 CALL
clacpy(
'Full', m, m,
x, ldx, xf, ldx )
304 CALL
cgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
305 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
307 CALL
cgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
308 $ u1, ldu1, work, ldx, zero, xf, ldx )
311 xf(i,i) = xf(i,i) - one
314 xf(min(p,q)-r+i,min(p,q)-r+i) =
315 $ xf(min(p,q)-r+i,min(p,q)-r+i) - cmplx( cos(theta(i)),
319 CALL
cgemm(
'No transpose',
'Conjugate transpose', p, m-q, m-q,
320 $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
322 CALL
cgemm(
'Conjugate transpose',
'No transpose', p, m-q, p,
323 $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
325 DO i = 1, min(p,m-q)-r
326 xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
329 xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
330 $ xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
331 $ cmplx( sin(theta(r-i+1)), 0.0e0 )
334 CALL
cgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
335 $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
337 CALL
cgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
338 $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
340 DO i = 1, min(m-p,q)-r
341 xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
344 xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
345 $ xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
346 $ cmplx( sin(theta(r-i+1)), 0.0e0 )
349 CALL
cgemm(
'No transpose',
'Conjugate transpose', m-p, m-q, m-q,
350 $ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
352 CALL
cgemm(
'Conjugate transpose',
'No transpose', m-p, m-q, m-p,
353 $ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
355 DO i = 1, min(m-p,m-q)-r
356 xf(p+i,q+i) = xf(p+i,q+i) - one
359 xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
360 $ xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
361 $ cmplx( cos(theta(i)), 0.0e0 )
366 resid =
clange(
'1', p, q, xf, ldx, rwork )
367 result( 1 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
371 resid =
clange(
'1', p, m-q, xf(1,q+1), ldx, rwork )
372 result( 2 ) = ( resid /
REAL(MAX(1,P,M-Q)) ) / eps2
376 resid =
clange(
'1', m-p, q, xf(p+1,1), ldx, rwork )
377 result( 3 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
381 resid =
clange(
'1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
382 result( 4 ) = ( resid /
REAL(MAX(1,M-P,M-Q)) ) / eps2
386 CALL
claset(
'Full', p, p, zero, one, work, ldu1 )
387 CALL
cherk(
'Upper',
'Conjugate transpose', p, p, -realone,
388 $ u1, ldu1, realone, work, ldu1 )
392 resid =
clanhe(
'1',
'Upper', p, work, ldu1, rwork )
393 result( 5 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
397 CALL
claset(
'Full', m-p, m-p, zero, one, work, ldu2 )
398 CALL
cherk(
'Upper',
'Conjugate transpose', m-p, m-p, -realone,
399 $ u2, ldu2, realone, work, ldu2 )
403 resid =
clanhe(
'1',
'Upper', m-p, work, ldu2, rwork )
404 result( 6 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
408 CALL
claset(
'Full', q, q, zero, one, work, ldv1t )
409 CALL
cherk(
'Upper',
'No transpose', q, q, -realone,
410 $ v1t, ldv1t, realone, work, ldv1t )
414 resid =
clanhe(
'1',
'Upper', q, work, ldv1t, rwork )
415 result( 7 ) = ( resid /
REAL(MAX(1,Q)) ) / ulp
419 CALL
claset(
'Full', m-q, m-q, zero, one, work, ldv2t )
420 CALL
cherk(
'Upper',
'No transpose', m-q, m-q, -realone,
421 $ v2t, ldv2t, realone, work, ldv2t )
425 resid =
clanhe(
'1',
'Upper', m-q, work, ldv2t, rwork )
426 result( 8 ) = ( resid /
REAL(MAX(1,M-Q)) ) / ulp
430 result( 9 ) = realzero
432 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 )
THEN
436 IF ( theta(i).LT.theta(i-1) )
THEN
444 CALL
claset(
'Full', q, q, zero, one, work, ldx )
445 CALL
cherk(
'Upper',
'Conjugate transpose', q, m, -realone,
446 $
x, ldx, realone, work, ldx )
449 $
clange(
'1', q, q, work, ldx, rwork ) /
REAL( M ) )
453 r = min( p, m-p, q, m-q )
457 CALL
clacpy(
'Full', m, q,
x, ldx, xf, ldx )
461 CALL
cuncsd2by1(
'Y',
'Y',
'Y', m, p, q, xf(1,1), ldx, xf(p+1,1),
462 $ ldx, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work,
463 $ lwork, rwork, 17*(r+2), iwork, info )
467 CALL
cgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
468 $
x, ldx, v1t, ldv1t, zero, work, ldx )
470 CALL
cgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
471 $ u1, ldu1, work, ldx, zero,
x, ldx )
474 x(i,i) =
x(i,i) - one
477 x(min(p,q)-r+i,min(p,q)-r+i) =
478 $
x(min(p,q)-r+i,min(p,q)-r+i) - cmplx( cos(theta(i)),
482 CALL
cgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
483 $
x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
485 CALL
cgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
486 $ one, u2, ldu2, work, ldx, zero,
x(p+1,1), ldx )
488 DO i = 1, min(m-p,q)-r
489 x(m-i+1,q-i+1) =
x(m-i+1,q-i+1) - one
492 x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
493 $
x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
494 $ cmplx( sin(theta(r-i+1)), 0.0e0 )
499 resid =
clange(
'1', p, q,
x, ldx, rwork )
500 result( 10 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
504 resid =
clange(
'1', m-p, q,
x(p+1,1), ldx, rwork )
505 result( 11 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
509 CALL
claset(
'Full', p, p, zero, one, work, ldu1 )
510 CALL
cherk(
'Upper',
'Conjugate transpose', p, p, -realone,
511 $ u1, ldu1, realone, work, ldu1 )
515 resid =
clanhe(
'1',
'Upper', p, work, ldu1, rwork )
516 result( 12 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
520 CALL
claset(
'Full', m-p, m-p, zero, one, work, ldu2 )
521 CALL
cherk(
'Upper',
'Conjugate transpose', m-p, m-p, -realone,
522 $ u2, ldu2, realone, work, ldu2 )
526 resid =
clanhe(
'1',
'Upper', m-p, work, ldu2, rwork )
527 result( 13 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
531 CALL
claset(
'Full', q, q, zero, one, work, ldv1t )
532 CALL
cherk(
'Upper',
'No transpose', q, q, -realone,
533 $ v1t, ldv1t, realone, work, ldv1t )
537 resid =
clanhe(
'1',
'Upper', q, work, ldv1t, rwork )
538 result( 14 ) = ( resid /
REAL(MAX(1,Q)) ) / ulp
542 result( 15 ) = realzero
544 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 )
THEN
545 result( 15 ) = ulpinv
548 IF ( theta(i).LT.theta(i-1) )
THEN
549 result( 15 ) = ulpinv
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...
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
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.
REAL function slamch(CMACH)
SLAMCH
subroutine ccsdts(M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, RWORK, RESULT)
CCSDTS
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
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 ...
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK