LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 231 of file scsdts.f.

231 *
232 * -- LAPACK test routine (version 3.4.0) --
233 * -- LAPACK is a software package provided by Univ. of Tennessee, --
234 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
235 * November 2011
236 *
237 * .. Scalar Arguments ..
238  INTEGER ldx, ldu1, ldu2, ldv1t, ldv2t, lwork, m, p, q
239 * ..
240 * .. Array Arguments ..
241  INTEGER iwork( * )
242  REAL result( 15 ), rwork( * ), theta( * )
243  REAL u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
244  $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
245  $ xf( ldx, * )
246 * ..
247 *
248 * =====================================================================
249 *
250 * .. Parameters ..
251  REAL piover2, realone, realzero
252  parameter ( piover2 = 1.57079632679489662e0,
253  $ realone = 1.0e0, realzero = 0.0e0 )
254  REAL zero, one
255  parameter ( zero = 0.0e0, one = 1.0e0 )
256 * ..
257 * .. Local Scalars ..
258  INTEGER i, info, r
259  REAL eps2, resid, ulp, ulpinv
260 * ..
261 * .. External Functions ..
262  REAL slamch, slange, slansy
263  EXTERNAL slamch, slange, slansy
264 * ..
265 * .. External Subroutines ..
266  EXTERNAL sgemm, slacpy, slaset, sorcsd, sorcsd2by1,
267  $ ssyrk
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC cos, max, min, REAL, sin
271 * ..
272 * .. Executable Statements ..
273 *
274  ulp = slamch( 'Precision' )
275  ulpinv = realone / ulp
276 *
277 * The first half of the routine checks the 2-by-2 CSD
278 *
279  CALL slaset( 'Full', m, m, zero, one, work, ldx )
280  CALL ssyrk( 'Upper', 'Conjugate transpose', m, m, -one, x, ldx,
281  $ one, work, ldx )
282  IF (m.GT.0) THEN
283  eps2 = max( ulp,
284  $ slange( '1', m, m, work, ldx, rwork ) / REAL( M ) )
285  ELSE
286  eps2 = ulp
287  END IF
288  r = min( p, m-p, q, m-q )
289 *
290 * Copy the matrix X to the array XF.
291 *
292  CALL slacpy( 'Full', m, m, x, ldx, xf, ldx )
293 *
294 * Compute the CSD
295 *
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 )
300 *
301 * Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
302 *
303  CALL slacpy( 'Full', m, m, x, ldx, xf, ldx )
304 *
305  CALL sgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
306  $ xf, ldx, v1t, ldv1t, zero, work, ldx )
307 *
308  CALL sgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
309  $ u1, ldu1, work, ldx, zero, xf, ldx )
310 *
311  DO i = 1, min(p,q)-r
312  xf(i,i) = xf(i,i) - one
313  END DO
314  DO i = 1, r
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))
317  END DO
318 *
319  CALL sgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
320  $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
321 *
322  CALL sgemm( 'Conjugate transpose', 'No transpose', p, m-q, p,
323  $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
324 *
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
327  END DO
328  DO i = 1, r
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) +
331  $ sin(theta(r-i+1))
332  END DO
333 *
334  CALL sgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
335  $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
336 *
337  CALL sgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
338  $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
339 *
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
342  END DO
343  DO i = 1, r
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) -
346  $ sin(theta(r-i+1))
347  END DO
348 *
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 )
351 *
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 )
354 *
355  DO i = 1, min(m-p,m-q)-r
356  xf(p+i,q+i) = xf(p+i,q+i) - one
357  END DO
358  DO i = 1, r
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) -
361  $ cos(theta(i))
362  END DO
363 *
364 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
365 *
366  resid = slange( '1', p, q, xf, ldx, rwork )
367  result( 1 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
368 *
369 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
370 *
371  resid = slange( '1', p, m-q, xf(1,q+1), ldx, rwork )
372  result( 2 ) = ( resid / REAL(MAX(1,P,M-Q)) ) / eps2
373 *
374 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
375 *
376  resid = slange( '1', m-p, q, xf(p+1,1), ldx, rwork )
377  result( 3 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
378 *
379 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
380 *
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
383 *
384 * Compute I - U1'*U1
385 *
386  CALL slaset( 'Full', p, p, zero, one, work, ldu1 )
387  CALL ssyrk( 'Upper', 'Conjugate transpose', p, p, -one, u1, ldu1,
388  $ one, work, ldu1 )
389 *
390 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
391 *
392  resid = slansy( '1', 'Upper', p, work, ldu1, rwork )
393  result( 5 ) = ( resid / REAL(MAX(1,P)) ) / ulp
394 *
395 * Compute I - U2'*U2
396 *
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 )
400 *
401 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
402 *
403  resid = slansy( '1', 'Upper', m-p, work, ldu2, rwork )
404  result( 6 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
405 *
406 * Compute I - V1T*V1T'
407 *
408  CALL slaset( 'Full', q, q, zero, one, work, ldv1t )
409  CALL ssyrk( 'Upper', 'No transpose', q, q, -one, v1t, ldv1t, one,
410  $ work, ldv1t )
411 *
412 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
413 *
414  resid = slansy( '1', 'Upper', q, work, ldv1t, rwork )
415  result( 7 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
416 *
417 * Compute I - V2T*V2T'
418 *
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,
421  $ one, work, ldv2t )
422 *
423 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
424 *
425  resid = slansy( '1', 'Upper', m-q, work, ldv2t, rwork )
426  result( 8 ) = ( resid / REAL(MAX(1,M-Q)) ) / ulp
427 *
428 * Check sorting
429 *
430  result( 9 ) = realzero
431  DO i = 1, r
432  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
433  result( 9 ) = ulpinv
434  END IF
435  IF( i.GT.1 ) THEN
436  IF ( theta(i).LT.theta(i-1) ) THEN
437  result( 9 ) = ulpinv
438  END IF
439  END IF
440  END DO
441 *
442 * The second half of the routine checks the 2-by-1 CSD
443 *
444  CALL slaset( 'Full', q, q, zero, one, work, ldx )
445  CALL ssyrk( 'Upper', 'Conjugate transpose', q, m, -one, x, ldx,
446  $ one, work, ldx )
447  IF (m.GT.0) THEN
448  eps2 = max( ulp,
449  $ slange( '1', q, q, work, ldx, rwork ) / REAL( M ) )
450  ELSE
451  eps2 = ulp
452  END IF
453  r = min( p, m-p, q, m-q )
454 *
455 * Copy the matrix [X11;X21] to the array XF.
456 *
457  CALL slacpy( 'Full', m, q, x, ldx, xf, ldx )
458 *
459 * Compute the CSD
460 *
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 )
464 *
465 * Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
466 *
467  CALL sgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
468  $ x, ldx, v1t, ldv1t, zero, work, ldx )
469 *
470  CALL sgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
471  $ u1, ldu1, work, ldx, zero, x, ldx )
472 *
473  DO i = 1, min(p,q)-r
474  x(i,i) = x(i,i) - one
475  END DO
476  DO i = 1, r
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))
479  END DO
480 *
481  CALL sgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
482  $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
483 *
484  CALL sgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
485  $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
486 *
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
489  END DO
490  DO i = 1, r
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) -
493  $ sin(theta(r-i+1))
494  END DO
495 *
496 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
497 *
498  resid = slange( '1', p, q, x, ldx, rwork )
499  result( 10 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
500 *
501 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
502 *
503  resid = slange( '1', m-p, q, x(p+1,1), ldx, rwork )
504  result( 11 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
505 *
506 * Compute I - U1'*U1
507 *
508  CALL slaset( 'Full', p, p, zero, one, work, ldu1 )
509  CALL ssyrk( 'Upper', 'Conjugate transpose', p, p, -one, u1, ldu1,
510  $ one, work, ldu1 )
511 *
512 * Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
513 *
514  resid = slansy( '1', 'Upper', p, work, ldu1, rwork )
515  result( 12 ) = ( resid / REAL(MAX(1,P)) ) / ulp
516 *
517 * Compute I - U2'*U2
518 *
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 )
522 *
523 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
524 *
525  resid = slansy( '1', 'Upper', m-p, work, ldu2, rwork )
526  result( 13 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
527 *
528 * Compute I - V1T*V1T'
529 *
530  CALL slaset( 'Full', q, q, zero, one, work, ldv1t )
531  CALL ssyrk( 'Upper', 'No transpose', q, q, -one, v1t, ldv1t, one,
532  $ work, ldv1t )
533 *
534 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
535 *
536  resid = slansy( '1', 'Upper', q, work, ldv1t, rwork )
537  result( 14 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
538 *
539 * Check sorting
540 *
541  result( 15 ) = realzero
542  DO i = 1, r
543  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
544  result( 15 ) = ulpinv
545  END IF
546  IF( i.GT.1 ) THEN
547  IF ( theta(i).LT.theta(i-1) ) THEN
548  result( 15 ) = ulpinv
549  END IF
550  END IF
551  END DO
552 *
553  RETURN
554 *
555 * End of SCSDTS
556 *
ssyrk
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:171
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
sorcsd
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:302
slacpy
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
slaset
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:112
slange
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:116
sorcsd2by1
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:234
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
slansy
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.
Definition: slansy.f:124

Here is the call graph for this function:

Here is the caller graph for this function: