238 INTEGER ldx, ldu1, ldu2, ldv1t, ldv2t, lwork, m, p, q
242 REAL result( 15 ), rwork( * ), theta( * )
243 REAL u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
244 $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
251 REAL piover2, realone, realzero
252 parameter ( piover2 = 1.57079632679489662e0,
253 $ realone = 1.0e0, realzero = 0.0e0 )
255 parameter ( zero = 0.0e0, one = 1.0e0 )
259 REAL eps2, resid, ulp, ulpinv
270 INTRINSIC cos, max, min,
REAL, sin
274 ulp =
slamch(
'Precision' )
275 ulpinv = realone / ulp
279 CALL slaset(
'Full', m, m, zero, one, work, ldx )
280 CALL ssyrk(
'Upper',
'Conjugate transpose', m, m, -one, x, ldx,
284 $
slange(
'1', m, m, work, ldx, rwork ) /
REAL( M ) )
288 r = min( p, m-p, q, m-q )
292 CALL slacpy(
'Full', m, m, x, ldx, xf, ldx )
296 CALL sorcsd(
'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 slacpy(
'Full', m, m, x, ldx, xf, ldx )
305 CALL sgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
306 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
308 CALL sgemm(
'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 sgemm(
'No transpose',
'Conjugate transpose', p, m-q, m-q,
320 $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
322 CALL sgemm(
'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 sgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
335 $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
337 CALL sgemm(
'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 sgemm(
'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 sgemm(
'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 =
slange(
'1', p, q, xf, ldx, rwork )
367 result( 1 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
371 resid =
slange(
'1', p, m-q, xf(1,q+1), ldx, rwork )
372 result( 2 ) = ( resid /
REAL(MAX(1,P,M-Q)) ) / eps2
376 resid =
slange(
'1', m-p, q, xf(p+1,1), ldx, rwork )
377 result( 3 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
381 resid =
slange(
'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 slaset(
'Full', p, p, zero, one, work, ldu1 )
387 CALL ssyrk(
'Upper',
'Conjugate transpose', p, p, -one, u1, ldu1,
392 resid =
slansy(
'1',
'Upper', p, work, ldu1, rwork )
393 result( 5 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
397 CALL slaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
398 CALL ssyrk(
'Upper',
'Conjugate transpose', m-p, m-p, -one, u2,
399 $ ldu2, one, work, ldu2 )
403 resid =
slansy(
'1',
'Upper', m-p, work, ldu2, rwork )
404 result( 6 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
408 CALL slaset(
'Full', q, q, zero, one, work, ldv1t )
409 CALL ssyrk(
'Upper',
'No transpose', q, q, -one, v1t, ldv1t, one,
414 resid =
slansy(
'1',
'Upper', q, work, ldv1t, rwork )
415 result( 7 ) = ( resid /
REAL(MAX(1,Q)) ) / ulp
419 CALL slaset(
'Full', m-q, m-q, zero, one, work, ldv2t )
420 CALL ssyrk(
'Upper',
'No transpose', m-q, m-q, -one, v2t, ldv2t,
425 resid =
slansy(
'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 slaset(
'Full', q, q, zero, one, work, ldx )
445 CALL ssyrk(
'Upper',
'Conjugate transpose', q, m, -one, x, ldx,
449 $
slange(
'1', q, q, work, ldx, rwork ) /
REAL( M ) )
453 r = min( p, m-p, q, m-q )
457 CALL slacpy(
'Full', m, q, x, ldx, xf, ldx )
461 CALL sorcsd2by1(
'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 sgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
468 $ x, ldx, v1t, ldv1t, zero, work, ldx )
470 CALL sgemm(
'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 sgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
482 $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
484 CALL sgemm(
'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 =
slange(
'1', p, q, x, ldx, rwork )
499 result( 10 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
503 resid =
slange(
'1', m-p, q, x(p+1,1), ldx, rwork )
504 result( 11 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
508 CALL slaset(
'Full', p, p, zero, one, work, ldu1 )
509 CALL ssyrk(
'Upper',
'Conjugate transpose', p, p, -one, u1, ldu1,
514 resid =
slansy(
'1',
'Upper', p, work, ldu1, rwork )
515 result( 12 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
519 CALL slaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
520 CALL ssyrk(
'Upper',
'Conjugate transpose', m-p, m-p, -one, u2,
521 $ ldu2, one, work, ldu2 )
525 resid =
slansy(
'1',
'Upper', m-p, work, ldu2, rwork )
526 result( 13 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
530 CALL slaset(
'Full', q, q, zero, one, work, ldv1t )
531 CALL ssyrk(
'Upper',
'No transpose', q, q, -one, v1t, ldv1t, one,
536 resid =
slansy(
'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 ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
real function slamch(CMACH)
SLAMCH
subroutine sorcsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, IWORK, INFO)
SORCSD2BY1
recursive subroutine sorcsd(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)
SORCSD
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
real function slansy(NORM, UPLO, N, A, LDA, WORK)
SLANSY 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.