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
262 REAL SLAMCH, CLANGE, CLANHE
263 EXTERNAL slamch, clange, clanhe
270 INTRINSIC cmplx, cos, max, min,
REAL, SIN
274 ulp = slamch(
'Precision' )
275 ulpinv = realone / ulp
279 CALL claset(
'Full', m, m, zero, one, work, ldx )
280 CALL cherk(
'Upper',
'Conjugate transpose', m, m, -realone,
281 $ x, ldx, realone, work, ldx )
284 $ clange(
'1', m, m, work, ldx, rwork ) /
REAL( M ) )
288 r = min( p, m-p, q, m-q )
292 CALL clacpy(
'Full', m, m, x, ldx, xf, ldx )
296 CALL cuncsd(
'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 )
303 CALL clacpy(
'Full', m, m, x, ldx, xf, ldx )
305 CALL cgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
306 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
308 CALL cgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
309 $ u1, ldu1, work, ldx, zero, xf, ldx )
312 xf(i,i) = xf(i,i) - one
315 xf(min(p,q)-r+i,min(p,q)-r+i) =
316 $ xf(min(p,q)-r+i,min(p,q)-r+i) - cmplx( cos(theta(i)),
320 CALL cgemm(
'No transpose',
'Conjugate transpose', p, m-q, m-q,
321 $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
323 CALL cgemm(
'Conjugate transpose',
'No transpose', p, m-q, p,
324 $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
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
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 $ cmplx( sin(theta(r-i+1)), 0.0e0 )
335 CALL cgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
336 $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
338 CALL cgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
339 $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
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
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 $ cmplx( sin(theta(r-i+1)), 0.0e0 )
350 CALL cgemm(
'No transpose',
'Conjugate transpose', m-p, m-q, m-q,
351 $ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
353 CALL cgemm(
'Conjugate transpose',
'No transpose', m-p, m-q, m-p,
354 $ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
356 DO i = 1, min(m-p,m-q)-r
357 xf(p+i,q+i) = xf(p+i,q+i) - one
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 $ cmplx( cos(theta(i)), 0.0e0 )
367 resid = clange(
'1', p, q, xf, ldx, rwork )
368 result( 1 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
372 resid = clange(
'1', p, m-q, xf(1,q+1), ldx, rwork )
373 result( 2 ) = ( resid /
REAL(MAX(1,P,M-Q)) ) / eps2
377 resid = clange(
'1', m-p, q, xf(p+1,1), ldx, rwork )
378 result( 3 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
382 resid = clange(
'1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
383 result( 4 ) = ( resid /
REAL(MAX(1,M-P,M-Q)) ) / eps2
387 CALL claset(
'Full', p, p, zero, one, work, ldu1 )
388 CALL cherk(
'Upper',
'Conjugate transpose', p, p, -realone,
389 $ u1, ldu1, realone, work, ldu1 )
393 resid = clanhe(
'1',
'Upper', p, work, ldu1, rwork )
394 result( 5 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
398 CALL claset(
'Full', m-p, m-p, zero, one, work, ldu2 )
399 CALL cherk(
'Upper',
'Conjugate transpose', m-p, m-p, -realone,
400 $ u2, ldu2, realone, work, ldu2 )
404 resid = clanhe(
'1',
'Upper', m-p, work, ldu2, rwork )
405 result( 6 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
409 CALL claset(
'Full', q, q, zero, one, work, ldv1t )
410 CALL cherk(
'Upper',
'No transpose', q, q, -realone,
411 $ v1t, ldv1t, realone, work, ldv1t )
415 resid = clanhe(
'1',
'Upper', q, work, ldv1t, rwork )
416 result( 7 ) = ( resid /
REAL(MAX(1,Q)) ) / ulp
420 CALL claset(
'Full', m-q, m-q, zero, one, work, ldv2t )
421 CALL cherk(
'Upper',
'No transpose', m-q, m-q, -realone,
422 $ v2t, ldv2t, realone, work, ldv2t )
426 resid = clanhe(
'1',
'Upper', m-q, work, ldv2t, rwork )
427 result( 8 ) = ( resid /
REAL(MAX(1,M-Q)) ) / ulp
431 result( 9 ) = realzero
433 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 )
THEN
437 IF ( theta(i).LT.theta(i-1) )
THEN
445 CALL claset(
'Full', q, q, zero, one, work, ldx )
446 CALL cherk(
'Upper',
'Conjugate transpose', q, m, -realone,
447 $ x, ldx, realone, work, ldx )
450 $ clange(
'1', q, q, work, ldx, rwork ) /
REAL( M ) )
454 r = min( p, m-p, q, m-q )
458 CALL clacpy(
'Full', m, q, x, ldx, xf, ldx )
462 CALL cuncsd2by1(
'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 )
468 CALL cgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
469 $ x, ldx, v1t, ldv1t, zero, work, ldx )
471 CALL cgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
472 $ u1, ldu1, work, ldx, zero, x, ldx )
475 x(i,i) = x(i,i) - one
478 x(min(p,q)-r+i,min(p,q)-r+i) =
479 $ x(min(p,q)-r+i,min(p,q)-r+i) - cmplx( cos(theta(i)),
483 CALL cgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
484 $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
486 CALL cgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
487 $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
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
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 $ cmplx( sin(theta(r-i+1)), 0.0e0 )
500 resid = clange(
'1', p, q, x, ldx, rwork )
501 result( 10 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
505 resid = clange(
'1', m-p, q, x(p+1,1), ldx, rwork )
506 result( 11 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
510 CALL claset(
'Full', p, p, zero, one, work, ldu1 )
511 CALL cherk(
'Upper',
'Conjugate transpose', p, p, -realone,
512 $ u1, ldu1, realone, work, ldu1 )
516 resid = clanhe(
'1',
'Upper', p, work, ldu1, rwork )
517 result( 12 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
521 CALL claset(
'Full', m-p, m-p, zero, one, work, ldu2 )
522 CALL cherk(
'Upper',
'Conjugate transpose', m-p, m-p, -realone,
523 $ u2, ldu2, realone, work, ldu2 )
527 resid = clanhe(
'1',
'Upper', m-p, work, ldu2, rwork )
528 result( 13 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
532 CALL claset(
'Full', q, q, zero, one, work, ldv1t )
533 CALL cherk(
'Upper',
'No transpose', q, q, -realone,
534 $ v1t, ldv1t, realone, work, ldv1t )
538 resid = clanhe(
'1',
'Upper', q, work, ldv1t, rwork )
539 result( 14 ) = ( resid /
REAL(MAX(1,Q)) ) / ulp
543 result( 15 ) = realzero
545 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 )
THEN
546 result( 15 ) = ulpinv
549 IF ( theta(i).LT.theta(i-1) )
THEN
550 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...
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
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
subroutine ccsdts(M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, RWORK, RESULT)
CCSDTS
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
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
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM