LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ scsdts()

subroutine scsdts ( integer  M,
integer  P,
integer  Q,
real, dimension( ldx, * )  X,
real, dimension( ldx, * )  XF,
integer  LDX,
real, dimension( ldu1, * )  U1,
integer  LDU1,
real, dimension( ldu2, * )  U2,
integer  LDU2,
real, dimension( ldv1t, * )  V1T,
integer  LDV1T,
real, dimension( ldv2t, * )  V2T,
integer  LDV2T,
real, dimension( * )  THETA,
integer, dimension( * )  IWORK,
real, dimension( lwork )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
real, dimension( 15 )  RESULT 
)

SCSDTS

Purpose:
 SCSDTS tests SORCSD, which, given an M-by-M partitioned orthogonal
 matrix X,
              Q  M-Q
       X = [ X11 X12 ] P   ,
           [ X21 X22 ] M-P

 computes the CSD

       [ U1    ]**T * [ X11 X12 ] * [ V1    ]
       [    U2 ]      [ X21 X22 ]   [    V2 ]

                             [  I  0  0 |  0  0  0 ]
                             [  0  C  0 |  0 -S  0 ]
                             [  0  0  0 |  0  0 -I ]
                           = [---------------------] = [ D11 D12 ] .
                             [  0  0  0 |  I  0  0 ]   [ D21 D22 ]
                             [  0  S  0 |  0  C  0 ]
                             [  0  0  I |  0  0  0 ]

 and also SORCSD2BY1, which, given
          Q
       [ X11 ] P   ,
       [ X21 ] M-P

 computes the 2-by-1 CSD

                                     [  I  0  0 ]
                                     [  0  C  0 ]
                                     [  0  0  0 ]
       [ U1    ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
       [    U2 ]      [ X21 ]        [  0  0  0 ]   [ D21 ]
                                     [  0  S  0 ]
                                     [  0  0  I ]
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix X.  M >= 0.
[in]P
          P is INTEGER
          The number of rows of the matrix X11.  P >= 0.
[in]Q
          Q is INTEGER
          The number of columns of the matrix X11.  Q >= 0.
[in]X
          X is REAL array, dimension (LDX,M)
          The M-by-M matrix X.
[out]XF
          XF is REAL array, dimension (LDX,M)
          Details of the CSD of X, as returned by SORCSD;
          see SORCSD for further details.
[in]LDX
          LDX is INTEGER
          The leading dimension of the arrays X and XF.
          LDX >= max( 1,M ).
[out]U1
          U1 is REAL array, dimension(LDU1,P)
          The P-by-P orthogonal matrix U1.
[in]LDU1
          LDU1 is INTEGER
          The leading dimension of the array U1. LDU >= max(1,P).
[out]U2
          U2 is REAL array, dimension(LDU2,M-P)
          The (M-P)-by-(M-P) orthogonal matrix U2.
[in]LDU2
          LDU2 is INTEGER
          The leading dimension of the array U2. LDU >= max(1,M-P).
[out]V1T
          V1T is REAL array, dimension(LDV1T,Q)
          The Q-by-Q orthogonal matrix V1T.
[in]LDV1T
          LDV1T is INTEGER
          The leading dimension of the array V1T. LDV1T >=
          max(1,Q).
[out]V2T
          V2T is REAL array, dimension(LDV2T,M-Q)
          The (M-Q)-by-(M-Q) orthogonal matrix V2T.
[in]LDV2T
          LDV2T is INTEGER
          The leading dimension of the array V2T. LDV2T >=
          max(1,M-Q).
[out]THETA
          THETA is REAL array, dimension MIN(P,M-P,Q,M-Q)
          The CS values of X; the essentially diagonal matrices C and
          S are constructed from THETA; see subroutine SORCSD for
          details.
[out]IWORK
          IWORK is INTEGER array, dimension (M)
[out]WORK
          WORK is REAL array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK
[out]RWORK
          RWORK is REAL array
[out]RESULT
          RESULT is REAL array, dimension (15)
          The test ratios:
          First, the 2-by-2 CSD:
          RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
          RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
          RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
          RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
          RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
          RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
          RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
          RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
          RESULT(9) = 0        if THETA is in increasing order and
                               all angles are in [0,pi/2];
                    = ULPINV   otherwise.
          Then, the 2-by-1 CSD:
          RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
          RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
          RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
          RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
          RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
          RESULT(15) = 0        if THETA is in increasing order and
                                all angles are in [0,pi/2];
                     = ULPINV   otherwise.
          ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 226 of file scsdts.f.

229 *
230 * -- LAPACK test routine --
231 * -- LAPACK is a software package provided by Univ. of Tennessee, --
232 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233 *
234 * .. Scalar Arguments ..
235  INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
236 * ..
237 * .. Array Arguments ..
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 * .. Parameters ..
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 * .. Local Scalars ..
256  INTEGER I, INFO, R
257  REAL EPS2, RESID, ULP, ULPINV
258 * ..
259 * .. External Functions ..
260  REAL SLAMCH, SLANGE, SLANSY
261  EXTERNAL slamch, slange, slansy
262 * ..
263 * .. External Subroutines ..
264  EXTERNAL sgemm, slacpy, slaset, sorcsd, sorcsd2by1,
265  $ ssyrk
266 * ..
267 * .. Intrinsic Functions ..
268  INTRINSIC cos, max, min, real, sin
269 * ..
270 * .. Executable Statements ..
271 *
272  ulp = slamch( 'Precision' )
273  ulpinv = realone / ulp
274 *
275 * The first half of the routine checks the 2-by-2 CSD
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 * Copy the matrix X to the array XF.
289 *
290  CALL slacpy( 'Full', m, m, x, ldx, xf, ldx )
291 *
292 * Compute the CSD
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 * Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
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 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
363 *
364  resid = slange( '1', p, q, xf, ldx, rwork )
365  result( 1 ) = ( resid / real(max(1,p,q)) ) / eps2
366 *
367 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
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 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
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 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
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 * Compute I - U1'*U1
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 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
389 *
390  resid = slansy( '1', 'Upper', p, work, ldu1, rwork )
391  result( 5 ) = ( resid / real(max(1,p)) ) / ulp
392 *
393 * Compute I - U2'*U2
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 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
400 *
401  resid = slansy( '1', 'Upper', m-p, work, ldu2, rwork )
402  result( 6 ) = ( resid / real(max(1,m-p)) ) / ulp
403 *
404 * Compute I - V1T*V1T'
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 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
411 *
412  resid = slansy( '1', 'Upper', q, work, ldv1t, rwork )
413  result( 7 ) = ( resid / real(max(1,q)) ) / ulp
414 *
415 * Compute I - V2T*V2T'
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 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
422 *
423  resid = slansy( '1', 'Upper', m-q, work, ldv2t, rwork )
424  result( 8 ) = ( resid / real(max(1,m-q)) ) / ulp
425 *
426 * Check sorting
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 * The second half of the routine checks the 2-by-1 CSD
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 * Copy the matrix [X11;X21] to the array XF.
454 *
455  CALL slacpy( 'Full', m, q, x, ldx, xf, ldx )
456 *
457 * Compute the CSD
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 * Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
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 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
495 *
496  resid = slange( '1', p, q, x, ldx, rwork )
497  result( 10 ) = ( resid / real(max(1,p,q)) ) / eps2
498 *
499 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
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 * Compute I - U1'*U1
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 * Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
511 *
512  resid = slansy( '1', 'Upper', p, work, ldu1, rwork )
513  result( 12 ) = ( resid / real(max(1,p)) ) / ulp
514 *
515 * Compute I - U2'*U2
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 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
522 *
523  resid = slansy( '1', 'Upper', m-p, work, ldu2, rwork )
524  result( 13 ) = ( resid / real(max(1,m-p)) ) / ulp
525 *
526 * Compute I - V1T*V1T'
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 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
533 *
534  resid = slansy( '1', 'Upper', q, work, ldv1t, rwork )
535  result( 14 ) = ( resid / real(max(1,q)) ) / ulp
536 *
537 * Check sorting
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 * End of SCSDTS
554 *
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.
Definition: slaset.f:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
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 ...
Definition: slange.f:114
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
Definition: sorcsd.f:300
subroutine sorcsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, IWORK, INFO)
SORCSD2BY1
Definition: sorcsd2by1.f:233
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,...
Definition: slansy.f:122
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:169
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: