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
270 INTRINSIC cos, dble, dcmplx, max, min,
REAL, sin
274 ulp =
dlamch(
'Precision' )
275 ulpinv = realone / ulp
279 CALL zlaset(
'Full', m, m, zero, one, work, ldx )
280 CALL zherk(
'Upper',
'Conjugate transpose', m, m, -realone,
281 $ x, ldx, realone, work, ldx )
284 $
zlange(
'1', m, m, work, ldx, rwork ) / dble( m ) )
288 r = min( p, m-p, q, m-q )
292 CALL zlacpy(
'Full', m, m, x, ldx, xf, ldx )
296 CALL zuncsd(
'Y',
'Y',
'Y',
'Y',
'N',
'D', m, p, q, xf(1,1), ldx,
297 $ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
298 $ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
299 $ work, lwork, rwork, 17*(r+2), iwork, info )
303 CALL zlacpy(
'Full', m, m, x, ldx, xf, ldx )
305 CALL zgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
306 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
308 CALL zgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
309 $ u1, ldu1, work, ldx, zero, xf, ldx )
312 xf(i,i) = xf(i,i) - one
315 xf(min(p,q)-r+i,min(p,q)-r+i) =
316 $ xf(min(p,q)-r+i,min(p,q)-r+i) - dcmplx( cos(theta(i)),
320 CALL zgemm(
'No transpose',
'Conjugate transpose', p, m-q, m-q,
321 $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
323 CALL zgemm(
'Conjugate transpose',
'No transpose', p, m-q, p,
324 $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
326 DO i = 1, min(p,m-q)-r
327 xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
330 xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
331 $ xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
332 $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
335 CALL zgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
336 $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
338 CALL zgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
339 $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
341 DO i = 1, min(m-p,q)-r
342 xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
345 xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
346 $ xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
347 $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
350 CALL zgemm(
'No transpose',
'Conjugate transpose', m-p, m-q, m-q,
351 $ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
353 CALL zgemm(
'Conjugate transpose',
'No transpose', m-p, m-q, m-p,
354 $ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
356 DO i = 1, min(m-p,m-q)-r
357 xf(p+i,q+i) = xf(p+i,q+i) - one
360 xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
361 $ xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
362 $ dcmplx( cos(theta(i)), 0.0d0 )
367 resid =
zlange(
'1', p, q, xf, ldx, rwork )
368 result( 1 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
372 resid =
zlange(
'1', p, m-q, xf(1,q+1), ldx, rwork )
373 result( 2 ) = ( resid /
REAL(MAX(1,P,M-Q)) ) / eps2
377 resid =
zlange(
'1', m-p, q, xf(p+1,1), ldx, rwork )
378 result( 3 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
382 resid =
zlange(
'1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
383 result( 4 ) = ( resid /
REAL(MAX(1,M-P,M-Q)) ) / eps2
387 CALL zlaset(
'Full', p, p, zero, one, work, ldu1 )
388 CALL zherk(
'Upper',
'Conjugate transpose', p, p, -realone,
389 $ u1, ldu1, realone, work, ldu1 )
393 resid =
zlanhe(
'1',
'Upper', p, work, ldu1, rwork )
394 result( 5 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
398 CALL zlaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
399 CALL zherk(
'Upper',
'Conjugate transpose', m-p, m-p, -realone,
400 $ u2, ldu2, realone, work, ldu2 )
404 resid =
zlanhe(
'1',
'Upper', m-p, work, ldu2, rwork )
405 result( 6 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
409 CALL zlaset(
'Full', q, q, zero, one, work, ldv1t )
410 CALL zherk(
'Upper',
'No transpose', q, q, -realone,
411 $ v1t, ldv1t, realone, work, ldv1t )
415 resid =
zlanhe(
'1',
'Upper', q, work, ldv1t, rwork )
416 result( 7 ) = ( resid /
REAL(MAX(1,Q)) ) / ulp
420 CALL zlaset(
'Full', m-q, m-q, zero, one, work, ldv2t )
421 CALL zherk(
'Upper',
'No transpose', m-q, m-q, -realone,
422 $ v2t, ldv2t, realone, work, ldv2t )
426 resid =
zlanhe(
'1',
'Upper', m-q, work, ldv2t, rwork )
427 result( 8 ) = ( resid /
REAL(MAX(1,M-Q)) ) / ulp
431 result( 9 ) = realzero
433 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 )
THEN
437 IF ( theta(i).LT.theta(i-1) )
THEN
445 CALL zlaset(
'Full', q, q, zero, one, work, ldx )
446 CALL zherk(
'Upper',
'Conjugate transpose', q, m, -realone,
447 $ x, ldx, realone, work, ldx )
450 $
zlange(
'1', q, q, work, ldx, rwork ) / dble( m ) )
454 r = min( p, m-p, q, m-q )
458 CALL zlacpy(
'Full', m, m, x, ldx, xf, ldx )
462 CALL zuncsd2by1(
'Y',
'Y',
'Y', m, p, q, xf(1,1), ldx, xf(p+1,1),
463 $ ldx, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work,
464 $ lwork, rwork, 17*(r+2), iwork, info )
468 CALL zgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
469 $ x, ldx, v1t, ldv1t, zero, work, ldx )
471 CALL zgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
472 $ u1, ldu1, work, ldx, zero, x, ldx )
475 x(i,i) = x(i,i) - one
478 x(min(p,q)-r+i,min(p,q)-r+i) =
479 $ x(min(p,q)-r+i,min(p,q)-r+i) - dcmplx( cos(theta(i)),
483 CALL zgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
484 $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
486 CALL zgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
487 $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
489 DO i = 1, min(m-p,q)-r
490 x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
493 x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
494 $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
495 $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
500 resid =
zlange(
'1', p, q, x, ldx, rwork )
501 result( 10 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
505 resid =
zlange(
'1', m-p, q, x(p+1,1), ldx, rwork )
506 result( 11 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
510 CALL zlaset(
'Full', p, p, zero, one, work, ldu1 )
511 CALL zherk(
'Upper',
'Conjugate transpose', p, p, -realone,
512 $ u1, ldu1, realone, work, ldu1 )
516 resid =
zlanhe(
'1',
'Upper', p, work, ldu1, rwork )
517 result( 12 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
521 CALL zlaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
522 CALL zherk(
'Upper',
'Conjugate transpose', m-p, m-p, -realone,
523 $ u2, ldu2, realone, work, ldu2 )
527 resid =
zlanhe(
'1',
'Upper', m-p, work, ldu2, rwork )
528 result( 13 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
532 CALL zlaset(
'Full', q, q, zero, one, work, ldv1t )
533 CALL zherk(
'Upper',
'No transpose', q, q, -realone,
534 $ v1t, ldv1t, realone, work, ldv1t )
538 resid =
zlanhe(
'1',
'Upper', q, work, ldv1t, rwork )
539 result( 14 ) = ( resid /
REAL(MAX(1,Q)) ) / ulp
543 result( 15 ) = realzero
545 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 )
THEN
546 result( 15 ) = ulpinv
549 IF ( theta(i).LT.theta(i-1) )
THEN
550 result( 15 ) = ulpinv
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(CMACH)
DLAMCH
double precision function zlanhe(NORM, UPLO, N, A, LDA, WORK)
ZLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix.
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...
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
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
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
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK