203 SUBROUTINE scsdts( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
204 $ ldv1t, v2t, ldv2t, theta, iwork, work, lwork,
213 INTEGER ldx, ldu1, ldu2, ldv1t, ldv2t, lwork, m, p, q
217 REAL result( 9 ), rwork( * ), theta( * )
218 REAL u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
219 $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
226 REAL piover2, realone, realzero
227 parameter( piover2 = 1.57079632679489662e0,
228 $ realone = 1.0e0, realzero = 0.0e0 )
230 parameter( zero = 0.0e0, one = 1.0e0 )
234 REAL eps2, resid, ulp, ulpinv
244 INTRINSIC REAL, max, min
248 ulp =
slamch(
'Precision' )
249 ulpinv = realone / ulp
250 CALL
slaset(
'Full', m, m, zero, one, work, ldx )
251 CALL
ssyrk(
'Upper',
'Conjugate transpose', m, m, -one, x, ldx,
255 $
slange(
'1', m, m, work, ldx, rwork ) /
REAL( M ) )
259 r = min( p, m-p, q, m-q )
263 CALL
slacpy(
'Full', m, m, x, ldx, xf, ldx )
267 CALL
sorcsd(
'Y',
'Y',
'Y',
'Y',
'N',
'D', m, p, q, xf(1,1), ldx,
268 $ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
269 $ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
270 $ work, lwork, iwork, info )
274 CALL
sgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
275 $ x, ldx, v1t, ldv1t, zero, work, ldx )
277 CALL
sgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
278 $ u1, ldu1, work, ldx, zero, x, ldx )
281 x(i,i) = x(i,i) - one
284 x(min(p,q)-r+i,min(p,q)-r+i) =
285 $ x(min(p,q)-r+i,min(p,q)-r+i) - cos(theta(i))
288 CALL
sgemm(
'No transpose',
'Conjugate transpose', p, m-q, m-q,
289 $ one, x(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
291 CALL
sgemm(
'Conjugate transpose',
'No transpose', p, m-q, p,
292 $ one, u1, ldu1, work, ldx, zero, x(1,q+1), ldx )
294 DO i = 1, min(p,m-q)-r
295 x(p-i+1,m-i+1) = x(p-i+1,m-i+1) + one
298 x(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
299 $ x(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
303 CALL
sgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
304 $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
306 CALL
sgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
307 $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
309 DO i = 1, min(m-p,q)-r
310 x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
313 x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
314 $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
318 CALL
sgemm(
'No transpose',
'Conjugate transpose', m-p, m-q, m-q,
319 $ one, x(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
321 CALL
sgemm(
'Conjugate transpose',
'No transpose', m-p, m-q, m-p,
322 $ one, u2, ldu2, work, ldx, zero, x(p+1,q+1), ldx )
324 DO i = 1, min(m-p,m-q)-r
325 x(p+i,q+i) = x(p+i,q+i) - one
328 x(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
329 $ x(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
335 resid =
slange(
'1', p, q, x, ldx, rwork )
336 result( 1 ) = ( resid /
REAL(MAX(1,P,Q)) ) / eps2
340 resid =
slange(
'1', p, m-q, x(1,q+1), ldx, rwork )
341 result( 2 ) = ( resid /
REAL(MAX(1,P,M-Q)) ) / eps2
345 resid =
slange(
'1', m-p, q, x(p+1,1), ldx, rwork )
346 result( 3 ) = ( resid /
REAL(MAX(1,M-P,Q)) ) / eps2
350 resid =
slange(
'1', m-p, m-q, x(p+1,q+1), ldx, rwork )
351 result( 4 ) = ( resid /
REAL(MAX(1,M-P,M-Q)) ) / eps2
355 CALL
slaset(
'Full', p, p, zero, one, work, ldu1 )
356 CALL
ssyrk(
'Upper',
'Conjugate transpose', p, p, -one, u1, ldu1,
361 resid =
slansy(
'1',
'Upper', p, work, ldu1, rwork )
362 result( 5 ) = ( resid /
REAL(MAX(1,P)) ) / ulp
366 CALL
slaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
367 CALL
ssyrk(
'Upper',
'Conjugate transpose', m-p, m-p, -one, u2,
368 $ ldu2, one, work, ldu2 )
372 resid =
slansy(
'1',
'Upper', m-p, work, ldu2, rwork )
373 result( 6 ) = ( resid /
REAL(MAX(1,M-P)) ) / ulp
377 CALL
slaset(
'Full', q, q, zero, one, work, ldv1t )
378 CALL
ssyrk(
'Upper',
'No transpose', q, q, -one, v1t, ldv1t, one,
383 resid =
slansy(
'1',
'Upper', q, work, ldv1t, rwork )
384 result( 7 ) = ( resid /
REAL(MAX(1,Q)) ) / ulp
388 CALL
slaset(
'Full', m-q, m-q, zero, one, work, ldv2t )
389 CALL
ssyrk(
'Upper',
'No transpose', m-q, m-q, -one, v2t, ldv2t,
394 resid =
slansy(
'1',
'Upper', m-q, work, ldv2t, rwork )
395 result( 8 ) = ( resid /
REAL(MAX(1,M-Q)) ) / ulp
401 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 )
THEN
405 IF ( theta(i).LT.theta(i-1) )
THEN