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
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
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
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