228 SUBROUTINE zcsdts( 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 DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
243 COMPLEX*16 U1( ldu1, * ), U2( ldu2, * ), V1T( ldv1t, * ),
244 $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
251 DOUBLE PRECISION PIOVER2, REALONE, REALZERO
252 parameter( piover2 = 1.57079632679489662d0,
253 $ realone = 1.0d0, realzero = 0.0d0 )
255 parameter( zero = (0.0d0,0.0d0), one = (1.0d0,0.0d0) )
259 DOUBLE PRECISION EPS2, RESID, ULP, ULPINV
262 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
263 EXTERNAL dlamch, zlange, zlanhe
269 INTRINSIC cos, dble, dcmplx, max, min,
REAL, SIN
273 ulp = dlamch(
'Precision' )
274 ulpinv = realone / ulp
278 CALL
zlaset(
'Full', m, m, zero, one, work, ldx )
279 CALL
zherk(
'Upper',
'Conjugate transpose', m, m, -realone,
280 $ x, ldx, realone, work, ldx )
283 $ zlange(
'1', m, m, work, ldx, rwork ) / dble( m ) )
287 r = min( p, m-p, q, m-q )
291 CALL
zlacpy(
'Full', m, m, x, ldx, xf, ldx )
295 CALL
zuncsd(
'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
zlacpy(
'Full', m, m, x, ldx, xf, ldx )
304 CALL
zgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
305 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
307 CALL
zgemm(
'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) - dcmplx( cos(theta(i)),
319 CALL
zgemm(
'No transpose',
'Conjugate transpose', p, m-q, m-q,
320 $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
322 CALL
zgemm(
'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 $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
334 CALL
zgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
335 $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
337 CALL
zgemm(
'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 $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
349 CALL
zgemm(
'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
zgemm(
'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 $ dcmplx( cos(theta(i)), 0.0d0 )
366 resid = zlange(
'1', p, q, xf, ldx, rwork )
367 result( 1 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
371 resid = zlange(
'1', p, m-q, xf(1,q+1), ldx, rwork )
372 result( 2 ) = ( resid /
REAL(MAX(1,P,M-Q)) ) / eps2
376 resid = zlange(
'1', m-p, q, xf(p+1,1), ldx, rwork )
377 result( 3 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
381 resid = zlange(
'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
zlaset(
'Full', p, p, zero, one, work, ldu1 )
387 CALL
zherk(
'Upper',
'Conjugate transpose', p, p, -realone,
388 $ u1, ldu1, realone, work, ldu1 )
392 resid = zlanhe(
'1',
'Upper', p, work, ldu1, rwork )
393 result( 5 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
397 CALL
zlaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
398 CALL
zherk(
'Upper',
'Conjugate transpose', m-p, m-p, -realone,
399 $ u2, ldu2, realone, work, ldu2 )
403 resid = zlanhe(
'1',
'Upper', m-p, work, ldu2, rwork )
404 result( 6 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
408 CALL
zlaset(
'Full', q, q, zero, one, work, ldv1t )
409 CALL
zherk(
'Upper',
'No transpose', q, q, -realone,
410 $ v1t, ldv1t, realone, work, ldv1t )
414 resid = zlanhe(
'1',
'Upper', q, work, ldv1t, rwork )
415 result( 7 ) = ( resid /
REAL(MAX(1,Q)) ) / ulp
419 CALL
zlaset(
'Full', m-q, m-q, zero, one, work, ldv2t )
420 CALL
zherk(
'Upper',
'No transpose', m-q, m-q, -realone,
421 $ v2t, ldv2t, realone, work, ldv2t )
425 resid = zlanhe(
'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
zlaset(
'Full', q, q, zero, one, work, ldx )
445 CALL
zherk(
'Upper',
'Conjugate transpose', q, m, -realone,
446 $ x, ldx, realone, work, ldx )
449 $ zlange(
'1', q, q, work, ldx, rwork ) / dble( m ) )
453 r = min( p, m-p, q, m-q )
457 CALL
zlacpy(
'Full', m, m, x, ldx, xf, ldx )
461 CALL
zuncsd2by1(
'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
zgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
468 $ x, ldx, v1t, ldv1t, zero, work, ldx )
470 CALL
zgemm(
'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) - dcmplx( cos(theta(i)),
482 CALL
zgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
483 $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
485 CALL
zgemm(
'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 $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
499 resid = zlange(
'1', p, q, x, ldx, rwork )
500 result( 10 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
504 resid = zlange(
'1', m-p, q, x(p+1,1), ldx, rwork )
505 result( 11 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
509 CALL
zlaset(
'Full', p, p, zero, one, work, ldu1 )
510 CALL
zherk(
'Upper',
'Conjugate transpose', p, p, -realone,
511 $ u1, ldu1, realone, work, ldu1 )
515 resid = zlanhe(
'1',
'Upper', p, work, ldu1, rwork )
516 result( 12 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
520 CALL
zlaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
521 CALL
zherk(
'Upper',
'Conjugate transpose', m-p, m-p, -realone,
522 $ u2, ldu2, realone, work, ldu2 )
526 resid = zlanhe(
'1',
'Upper', m-p, work, ldu2, rwork )
527 result( 13 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
531 CALL
zlaset(
'Full', q, q, zero, one, work, ldv1t )
532 CALL
zherk(
'Upper',
'No transpose', q, q, -realone,
533 $ v1t, ldv1t, realone, work, ldv1t )
537 resid = zlanhe(
'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 zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zcsdts(M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, RWORK, RESULT)
ZCSDTS
recursive subroutine zuncsd(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)
ZUNCSD
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine zuncsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
ZUNCSD2BY1