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