229
230
231
232
233
234
235 INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
236
237
238 INTEGER IWORK( * )
239 REAL RESULT( 15 ), RWORK( * ), THETA( * )
240 REAL U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
241 $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
242 $ XF( LDX, * )
243
244
245
246
247
248 REAL REALONE, REALZERO
249 parameter( realone = 1.0e0, realzero = 0.0e0 )
250 REAL ZERO, ONE
251 parameter( zero = 0.0e0, one = 1.0e0 )
252 REAL PIOVER2
253 parameter( piover2 = 1.57079632679489661923132169163975144210e0 )
254
255
256 INTEGER I, INFO, R
257 REAL EPS2, RESID, ULP, ULPINV
258
259
260 REAL SLAMCH, SLANGE, SLANSY
262
263
266
267
268 INTRINSIC cos, max, min, real, sin
269
270
271
272 ulp =
slamch(
'Precision' )
273 ulpinv = realone / ulp
274
275
276
277 CALL slaset(
'Full', m, m, zero, one, work, ldx )
278 CALL ssyrk(
'Upper',
'Conjugate transpose', m, m, -one, x, ldx,
279 $ one, work, ldx )
280 IF (m.GT.0) THEN
281 eps2 = max( ulp,
282 $
slange(
'1', m, m, work, ldx, rwork ) / real( m ) )
283 ELSE
284 eps2 = ulp
285 END IF
286 r = min( p, m-p, q, m-q )
287
288
289
290 CALL slacpy(
'Full', m, m, x, ldx, xf, ldx )
291
292
293
294 CALL sorcsd(
'Y',
'Y',
'Y',
'Y',
'N',
'D', m, p, q, xf(1,1), ldx,
295 $ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
296 $ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
297 $ work, lwork, iwork, info )
298
299
300
301 CALL slacpy(
'Full', m, m, x, ldx, xf, ldx )
302
303 CALL sgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
304 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
305
306 CALL sgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
307 $ u1, ldu1, work, ldx, zero, xf, ldx )
308
309 DO i = 1, min(p,q)-r
310 xf(i,i) = xf(i,i) - one
311 END DO
312 DO i = 1, r
313 xf(min(p,q)-r+i,min(p,q)-r+i) =
314 $ xf(min(p,q)-r+i,min(p,q)-r+i) - cos(theta(i))
315 END DO
316
317 CALL sgemm(
'No transpose',
'Conjugate transpose', p, m-q, m-q,
318 $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
319
320 CALL sgemm(
'Conjugate transpose',
'No transpose', p, m-q, p,
321 $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
322
323 DO i = 1, min(p,m-q)-r
324 xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
325 END DO
326 DO i = 1, r
327 xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
328 $ xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
329 $ sin(theta(r-i+1))
330 END DO
331
332 CALL sgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
333 $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
334
335 CALL sgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
336 $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
337
338 DO i = 1, min(m-p,q)-r
339 xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
340 END DO
341 DO i = 1, r
342 xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
343 $ xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
344 $ sin(theta(r-i+1))
345 END DO
346
347 CALL sgemm(
'No transpose',
'Conjugate transpose', m-p, m-q, m-q,
348 $ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
349
350 CALL sgemm(
'Conjugate transpose',
'No transpose', m-p, m-q, m-p,
351 $ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
352
353 DO i = 1, min(m-p,m-q)-r
354 xf(p+i,q+i) = xf(p+i,q+i) - one
355 END DO
356 DO i = 1, r
357 xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
358 $ xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
359 $ cos(theta(i))
360 END DO
361
362
363
364 resid =
slange(
'1', p, q, xf, ldx, rwork )
365 result( 1 ) = ( resid / real(max(1,p,q)) ) / eps2
366
367
368
369 resid =
slange(
'1', p, m-q, xf(1,q+1), ldx, rwork )
370 result( 2 ) = ( resid / real(max(1,p,m-q)) ) / eps2
371
372
373
374 resid =
slange(
'1', m-p, q, xf(p+1,1), ldx, rwork )
375 result( 3 ) = ( resid / real(max(1,m-p,q)) ) / eps2
376
377
378
379 resid =
slange(
'1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
380 result( 4 ) = ( resid / real(max(1,m-p,m-q)) ) / eps2
381
382
383
384 CALL slaset(
'Full', p, p, zero, one, work, ldu1 )
385 CALL ssyrk(
'Upper',
'Conjugate transpose', p, p, -one, u1, ldu1,
386 $ one, work, ldu1 )
387
388
389
390 resid =
slansy(
'1',
'Upper', p, work, ldu1, rwork )
391 result( 5 ) = ( resid / real(max(1,p)) ) / ulp
392
393
394
395 CALL slaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
396 CALL ssyrk(
'Upper',
'Conjugate transpose', m-p, m-p, -one, u2,
397 $ ldu2, one, work, ldu2 )
398
399
400
401 resid =
slansy(
'1',
'Upper', m-p, work, ldu2, rwork )
402 result( 6 ) = ( resid / real(max(1,m-p)) ) / ulp
403
404
405
406 CALL slaset(
'Full', q, q, zero, one, work, ldv1t )
407 CALL ssyrk(
'Upper',
'No transpose', q, q, -one, v1t, ldv1t, one,
408 $ work, ldv1t )
409
410
411
412 resid =
slansy(
'1',
'Upper', q, work, ldv1t, rwork )
413 result( 7 ) = ( resid / real(max(1,q)) ) / ulp
414
415
416
417 CALL slaset(
'Full', m-q, m-q, zero, one, work, ldv2t )
418 CALL ssyrk(
'Upper',
'No transpose', m-q, m-q, -one, v2t, ldv2t,
419 $ one, work, ldv2t )
420
421
422
423 resid =
slansy(
'1',
'Upper', m-q, work, ldv2t, rwork )
424 result( 8 ) = ( resid / real(max(1,m-q)) ) / ulp
425
426
427
428 result( 9 ) = realzero
429 DO i = 1, r
430 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
431 result( 9 ) = ulpinv
432 END IF
433 IF( i.GT.1 ) THEN
434 IF ( theta(i).LT.theta(i-1) ) THEN
435 result( 9 ) = ulpinv
436 END IF
437 END IF
438 END DO
439
440
441
442 CALL slaset(
'Full', q, q, zero, one, work, ldx )
443 CALL ssyrk(
'Upper',
'Conjugate transpose', q, m, -one, x, ldx,
444 $ one, work, ldx )
445 IF (m.GT.0) THEN
446 eps2 = max( ulp,
447 $
slange(
'1', q, q, work, ldx, rwork ) / real( m ) )
448 ELSE
449 eps2 = ulp
450 END IF
451 r = min( p, m-p, q, m-q )
452
453
454
455 CALL slacpy(
'Full', m, q, x, ldx, xf, ldx )
456
457
458
459 CALL sorcsd2by1(
'Y',
'Y',
'Y', m, p, q, xf(1,1), ldx, xf(p+1,1),
460 $ ldx, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work,
461 $ lwork, iwork, info )
462
463
464
465 CALL sgemm(
'No transpose',
'Conjugate transpose', p, q, q, one,
466 $ x, ldx, v1t, ldv1t, zero, work, ldx )
467
468 CALL sgemm(
'Conjugate transpose',
'No transpose', p, q, p, one,
469 $ u1, ldu1, work, ldx, zero, x, ldx )
470
471 DO i = 1, min(p,q)-r
472 x(i,i) = x(i,i) - one
473 END DO
474 DO i = 1, r
475 x(min(p,q)-r+i,min(p,q)-r+i) =
476 $ x(min(p,q)-r+i,min(p,q)-r+i) - cos(theta(i))
477 END DO
478
479 CALL sgemm(
'No transpose',
'Conjugate transpose', m-p, q, q, one,
480 $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
481
482 CALL sgemm(
'Conjugate transpose',
'No transpose', m-p, q, m-p,
483 $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
484
485 DO i = 1, min(m-p,q)-r
486 x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
487 END DO
488 DO i = 1, r
489 x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
490 $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
491 $ sin(theta(r-i+1))
492 END DO
493
494
495
496 resid =
slange(
'1', p, q, x, ldx, rwork )
497 result( 10 ) = ( resid / real(max(1,p,q)) ) / eps2
498
499
500
501 resid =
slange(
'1', m-p, q, x(p+1,1), ldx, rwork )
502 result( 11 ) = ( resid / real(max(1,m-p,q)) ) / eps2
503
504
505
506 CALL slaset(
'Full', p, p, zero, one, work, ldu1 )
507 CALL ssyrk(
'Upper',
'Conjugate transpose', p, p, -one, u1, ldu1,
508 $ one, work, ldu1 )
509
510
511
512 resid =
slansy(
'1',
'Upper', p, work, ldu1, rwork )
513 result( 12 ) = ( resid / real(max(1,p)) ) / ulp
514
515
516
517 CALL slaset(
'Full', m-p, m-p, zero, one, work, ldu2 )
518 CALL ssyrk(
'Upper',
'Conjugate transpose', m-p, m-p, -one, u2,
519 $ ldu2, one, work, ldu2 )
520
521
522
523 resid =
slansy(
'1',
'Upper', m-p, work, ldu2, rwork )
524 result( 13 ) = ( resid / real(max(1,m-p)) ) / ulp
525
526
527
528 CALL slaset(
'Full', q, q, zero, one, work, ldv1t )
529 CALL ssyrk(
'Upper',
'No transpose', q, q, -one, v1t, ldv1t, one,
530 $ work, ldv1t )
531
532
533
534 resid =
slansy(
'1',
'Upper', q, work, ldv1t, rwork )
535 result( 14 ) = ( resid / real(max(1,q)) ) / ulp
536
537
538
539 result( 15 ) = realzero
540 DO i = 1, r
541 IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
542 result( 15 ) = ulpinv
543 END IF
544 IF( i.GT.1 ) THEN
545 IF ( theta(i).LT.theta(i-1) ) THEN
546 result( 15 ) = ulpinv
547 END IF
548 END IF
549 END DO
550
551 RETURN
552
553
554
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
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 ...
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,...
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.
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