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
subroutine ccsdts(m, p, q, x, xf, ldx, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t, theta, iwork, work, lwork, rwork, result)
CCSDTS
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
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
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 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
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