228 SUBROUTINE dcsdts( 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 DOUBLE PRECISION 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 )
254 DOUBLE PRECISION zero, one
255 parameter( zero = 0.0d0, one = 1.0d0 )
259 DOUBLE PRECISION eps2, resid, ulp, ulpinv
270 INTRINSIC cos, dble, max, min, sin
274 ulp =
dlamch(
'Precision' )
275 ulpinv = realone / ulp
279 CALL
dlaset(
'Full', m, m, zero, one, work, ldx )
280 CALL
dsyrk(
'Upper',
'Conjugate transpose', m, m, -one,
x, ldx,
284 $
dlange(
'1', m, m, work, ldx, rwork ) / dble( m ) )
288 r = min( p, m-p, q, m-q )
292 CALL
dlacpy(
'Full', m, m,
x, ldx, xf, ldx )
296 CALL
dorcsd(
'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, iwork, info )
303 CALL
dlacpy(
'Full', m, m,
x, ldx, xf, ldx )
305 CALL
dgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
306 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
308 CALL
dgemm(
'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) - cos(theta(i))
319 CALL
dgemm(
'No transpose',
'Conjugate transpose', p, m-q, m-q,
320 $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
322 CALL
dgemm(
'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) +
334 CALL
dgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
335 $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
337 CALL
dgemm(
'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) -
349 CALL
dgemm(
'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
dgemm(
'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) -
366 resid =
dlange(
'1', p, q, xf, ldx, rwork )
367 result( 1 ) = ( resid / dble(max(1,p,q)) ) / eps2
371 resid =
dlange(
'1', p, m-q, xf(1,q+1), ldx, rwork )
372 result( 2 ) = ( resid / dble(max(1,p,m-q)) ) / eps2
376 resid =
dlange(
'1', m-p, q, xf(p+1,1), ldx, rwork )
377 result( 3 ) = ( resid / dble(max(1,m-p,q)) ) / eps2
381 resid =
dlange(
'1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
382 result( 4 ) = ( resid / dble(max(1,m-p,m-q)) ) / eps2
386 CALL
dlaset(
'Full', p, p, zero, one, work, ldu1 )
387 CALL
dsyrk(
'Upper',
'Conjugate transpose', p, p, -one, u1, ldu1,
392 resid =
dlansy(
'1',
'Upper', p, work, ldu1, rwork )
393 result( 5 ) = ( resid / dble(max(1,p)) ) / ulp
397 CALL
dlaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
398 CALL
dsyrk(
'Upper',
'Conjugate transpose', m-p, m-p, -one, u2,
399 $ ldu2, one, work, ldu2 )
403 resid =
dlansy(
'1',
'Upper', m-p, work, ldu2, rwork )
404 result( 6 ) = ( resid / dble(max(1,m-p)) ) / ulp
408 CALL
dlaset(
'Full', q, q, zero, one, work, ldv1t )
409 CALL
dsyrk(
'Upper',
'No transpose', q, q, -one, v1t, ldv1t, one,
414 resid =
dlansy(
'1',
'Upper', q, work, ldv1t, rwork )
415 result( 7 ) = ( resid / dble(max(1,q)) ) / ulp
419 CALL
dlaset(
'Full', m-q, m-q, zero, one, work, ldv2t )
420 CALL
dsyrk(
'Upper',
'No transpose', m-q, m-q, -one, v2t, ldv2t,
425 resid =
dlansy(
'1',
'Upper', m-q, work, ldv2t, rwork )
426 result( 8 ) = ( resid / dble(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
dlaset(
'Full', q, q, zero, one, work, ldx )
445 CALL
dsyrk(
'Upper',
'Conjugate transpose', q, m, -one,
x, ldx,
449 $
dlange(
'1', q, q, work, ldx, rwork ) / dble( m ) )
453 r = min( p, m-p, q, m-q )
457 CALL
dlacpy(
'Full', m, q,
x, ldx, xf, ldx )
461 CALL
dorcsd2by1(
'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, iwork, info )
467 CALL
dgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
468 $
x, ldx, v1t, ldv1t, zero, work, ldx )
470 CALL
dgemm(
'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) - cos(theta(i))
481 CALL
dgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
482 $
x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
484 CALL
dgemm(
'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) -
498 resid =
dlange(
'1', p, q,
x, ldx, rwork )
499 result( 10 ) = ( resid / dble(max(1,p,q)) ) / eps2
503 resid =
dlange(
'1', m-p, q,
x(p+1,1), ldx, rwork )
504 result( 11 ) = ( resid / dble(max(1,m-p,q)) ) / eps2
508 CALL
dlaset(
'Full', p, p, zero, one, work, ldu1 )
509 CALL
dsyrk(
'Upper',
'Conjugate transpose', p, p, -one, u1, ldu1,
514 resid =
dlansy(
'1',
'Upper', p, work, ldu1, rwork )
515 result( 12 ) = ( resid / dble(max(1,p)) ) / ulp
519 CALL
dlaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
520 CALL
dsyrk(
'Upper',
'Conjugate transpose', m-p, m-p, -one, u2,
521 $ ldu2, one, work, ldu2 )
525 resid =
dlansy(
'1',
'Upper', m-p, work, ldu2, rwork )
526 result( 13 ) = ( resid / dble(max(1,m-p)) ) / ulp
530 CALL
dlaset(
'Full', q, q, zero, one, work, ldv1t )
531 CALL
dsyrk(
'Upper',
'No transpose', q, q, -one, v1t, ldv1t, one,
536 resid =
dlansy(
'1',
'Upper', q, work, ldv1t, rwork )
537 result( 14 ) = ( resid / dble(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
DOUBLE PRECISION function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
recursive subroutine dorcsd(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, IWORK, INFO)
DORCSD
subroutine dorcsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, IWORK, INFO)
DORCSD2BY1
DOUBLE PRECISION function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dcsdts(M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, RWORK, RESULT)
DCSDTS
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...