226 SUBROUTINE ccsdts( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
227 $ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
235 INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
239 REAL RESULT( 15 ), RWORK( * ), THETA( * )
240 COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
241 $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
248 REAL REALONE, REALZERO
249 PARAMETER ( REALONE = 1.0e0, realzero = 0.0e0 )
251 parameter( zero = (0.0e0,0.0e0), one = (1.0e0,0.0e0) )
253 parameter( piover2 = 1.57079632679489661923132169163975144210e0 )
257 REAL EPS2, RESID, ULP, ULPINV
260 REAL SLAMCH, CLANGE, CLANHE
261 EXTERNAL SLAMCH, CLANGE, CLANHE
268 INTRINSIC cmplx, cos, max, min, real, sin
272 ulp = slamch(
'Precision' )
273 ulpinv = realone / ulp
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 )
282 $ clange(
'1', m, m, work, ldx, rwork ) / real( m ) )
286 r = min( p, m-p, q, m-q )
290 CALL clacpy(
'Full', m, m, x, ldx, xf, ldx )
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 )
301 CALL clacpy(
'Full', m, m, x, ldx, xf, ldx )
303 CALL cgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
304 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
306 CALL cgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
307 $ u1, ldu1, work, ldx, zero, xf, ldx )
310 xf(i,i) = xf(i,i) - one
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)),
318 CALL cgemm(
'No transpose',
'Conjugate transpose', p, m-q, m-q,
319 $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
321 CALL cgemm(
'Conjugate transpose',
'No transpose', p, m-q, p,
322 $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
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
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 )
333 CALL cgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
334 $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
336 CALL cgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
337 $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
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
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 )
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 )
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 )
354 DO i = 1, min(m-p,m-q)-r
355 xf(p+i,q+i) = xf(p+i,q+i) - one
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 )
365 resid = clange(
'1', p, q, xf, ldx, rwork )
366 result( 1 ) = ( resid / real(max(1,p,q)) ) / eps2
370 resid = clange(
'1', p, m-q, xf(1,q+1), ldx, rwork )
371 result( 2 ) = ( resid / real(max(1,p,m-q)) ) / eps2
375 resid = clange(
'1', m-p, q, xf(p+1,1), ldx, rwork )
376 result( 3 ) = ( resid / real(max(1,m-p,q)) ) / eps2
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
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 )
391 resid = clanhe(
'1',
'Upper', p, work, ldu1, rwork )
392 result( 5 ) = ( resid / real(max(1,p)) ) / ulp
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 )
402 resid = clanhe(
'1',
'Upper', m-p, work, ldu2, rwork )
403 result( 6 ) = ( resid / real(max(1,m-p)) ) / ulp
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 )
413 resid = clanhe(
'1',
'Upper', q, work, ldv1t, rwork )
414 result( 7 ) = ( resid / real(max(1,q)) ) / ulp
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 )
424 resid = clanhe(
'1',
'Upper', m-q, work, ldv2t, rwork )
425 result( 8 ) = ( resid / real(max(1,m-q)) ) / ulp
429 result( 9 ) = realzero
431 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 )
THEN
435 IF ( theta(i).LT.theta(i-1) )
THEN
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 )
448 $ clange(
'1', q, q, work, ldx, rwork ) / real( m ) )
452 r = min( p, m-p, q, m-q )
456 CALL clacpy(
'Full', m, q, x, ldx, xf, ldx )
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 )
466 CALL cgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
467 $ x, ldx, v1t, ldv1t, zero, work, ldx )
469 CALL cgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
470 $ u1, ldu1, work, ldx, zero, x, ldx )
473 x(i,i) = x(i,i) - one
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)),
481 CALL cgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
482 $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
484 CALL cgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
485 $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
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
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 )
498 resid = clange(
'1', p, q, x, ldx, rwork )
499 result( 10 ) = ( resid / real(max(1,p,q)) ) / eps2
503 resid = clange(
'1', m-p, q, x(p+1,1), ldx, rwork )
504 result( 11 ) = ( resid / real(max(1,m-p,q)) ) / eps2
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 )
514 resid = clanhe(
'1',
'Upper', p, work, ldu1, rwork )
515 result( 12 ) = ( resid / real(max(1,p)) ) / ulp
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 )
525 resid = clanhe(
'1',
'Upper', m-p, work, ldu2, rwork )
526 result( 13 ) = ( resid / real(max(1,m-p)) ) / ulp
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 )
536 resid = clanhe(
'1',
'Upper', q, work, ldv1t, rwork )
537 result( 14 ) = ( resid / real(max(1,q)) ) / ulp
541 result( 15 ) = realzero
543 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 )
THEN
544 result( 15 ) = ulpinv
547 IF ( theta(i).LT.theta(i-1) )
THEN
548 result( 15 ) = ulpinv