LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dcsdts()

subroutine dcsdts ( integer  m,
integer  p,
integer  q,
double precision, dimension( ldx, * )  x,
double precision, dimension( ldx, * )  xf,
integer  ldx,
double precision, dimension( ldu1, * )  u1,
integer  ldu1,
double precision, dimension( ldu2, * )  u2,
integer  ldu2,
double precision, dimension( ldv1t, * )  v1t,
integer  ldv1t,
double precision, dimension( ldv2t, * )  v2t,
integer  ldv2t,
double precision, dimension( * )  theta,
integer, dimension( * )  iwork,
double precision, dimension( lwork )  work,
integer  lwork,
double precision, dimension( * )  rwork,
double precision, dimension( 15 )  result 
)

DCSDTS

Purpose:
 DCSDTS tests DORCSD, 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 DORCSD2BY1, 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 DOUBLE PRECISION array, dimension (LDX,M)
          The M-by-M matrix X.
[out]XF
          XF is DOUBLE PRECISION array, dimension (LDX,M)
          Details of the CSD of X, as returned by DORCSD;
          see DORCSD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DORCSD for
          details.
[out]IWORK
          IWORK is INTEGER array, dimension (M)
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK
[out]RWORK
          RWORK is DOUBLE PRECISION array
[out]RESULT
          RESULT is DOUBLE PRECISION 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 dcsdts.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 DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
240 DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
241 $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
242 $ XF( LDX, * )
243* ..
244*
245* =====================================================================
246*
247* .. Parameters ..
248 DOUBLE PRECISION REALONE, REALZERO
249 parameter( realone = 1.0d0, realzero = 0.0d0 )
250 DOUBLE PRECISION ZERO, ONE
251 parameter( zero = 0.0d0, one = 1.0d0 )
252 DOUBLE PRECISION PIOVER2
253 parameter( piover2 = 1.57079632679489661923132169163975144210d0 )
254* ..
255* .. Local Scalars ..
256 INTEGER I, INFO, R
257 DOUBLE PRECISION EPS2, RESID, ULP, ULPINV
258* ..
259* .. External Functions ..
260 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
261 EXTERNAL dlamch, dlange, dlansy
262* ..
263* .. External Subroutines ..
264 EXTERNAL dgemm, dlacpy, dlaset, dorcsd, dorcsd2by1,
265 $ dsyrk
266* ..
267* .. Intrinsic Functions ..
268 INTRINSIC cos, dble, max, min, sin
269* ..
270* .. Executable Statements ..
271*
272 ulp = dlamch( 'Precision' )
273 ulpinv = realone / ulp
274*
275* The first half of the routine checks the 2-by-2 CSD
276*
277 CALL dlaset( 'Full', m, m, zero, one, work, ldx )
278 CALL dsyrk( 'Upper', 'Conjugate transpose', m, m, -one, x, ldx,
279 $ one, work, ldx )
280 IF (m.GT.0) THEN
281 eps2 = max( ulp,
282 $ dlange( '1', m, m, work, ldx, rwork ) / dble( 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 dlacpy( 'Full', m, m, x, ldx, xf, ldx )
291*
292* Compute the CSD
293*
294 CALL dorcsd( '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 dlacpy( 'Full', m, m, x, ldx, xf, ldx )
302*
303 CALL dgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
304 $ xf, ldx, v1t, ldv1t, zero, work, ldx )
305*
306 CALL dgemm( '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 dgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
318 $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
319*
320 CALL dgemm( '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 dgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
333 $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
334*
335 CALL dgemm( '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 dgemm( '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 dgemm( '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 = dlange( '1', p, q, xf, ldx, rwork )
365 result( 1 ) = ( resid / dble(max(1,p,q)) ) / eps2
366*
367* Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
368*
369 resid = dlange( '1', p, m-q, xf(1,q+1), ldx, rwork )
370 result( 2 ) = ( resid / dble(max(1,p,m-q)) ) / eps2
371*
372* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
373*
374 resid = dlange( '1', m-p, q, xf(p+1,1), ldx, rwork )
375 result( 3 ) = ( resid / dble(max(1,m-p,q)) ) / eps2
376*
377* Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
378*
379 resid = dlange( '1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
380 result( 4 ) = ( resid / dble(max(1,m-p,m-q)) ) / eps2
381*
382* Compute I - U1'*U1
383*
384 CALL dlaset( 'Full', p, p, zero, one, work, ldu1 )
385 CALL dsyrk( '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 = dlansy( '1', 'Upper', p, work, ldu1, rwork )
391 result( 5 ) = ( resid / dble(max(1,p)) ) / ulp
392*
393* Compute I - U2'*U2
394*
395 CALL dlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
396 CALL dsyrk( '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 = dlansy( '1', 'Upper', m-p, work, ldu2, rwork )
402 result( 6 ) = ( resid / dble(max(1,m-p)) ) / ulp
403*
404* Compute I - V1T*V1T'
405*
406 CALL dlaset( 'Full', q, q, zero, one, work, ldv1t )
407 CALL dsyrk( '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 = dlansy( '1', 'Upper', q, work, ldv1t, rwork )
413 result( 7 ) = ( resid / dble(max(1,q)) ) / ulp
414*
415* Compute I - V2T*V2T'
416*
417 CALL dlaset( 'Full', m-q, m-q, zero, one, work, ldv2t )
418 CALL dsyrk( '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 = dlansy( '1', 'Upper', m-q, work, ldv2t, rwork )
424 result( 8 ) = ( resid / dble(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 dlaset( 'Full', q, q, zero, one, work, ldx )
443 CALL dsyrk( 'Upper', 'Conjugate transpose', q, m, -one, x, ldx,
444 $ one, work, ldx )
445 IF( m.GT.0 ) THEN
446 eps2 = max( ulp,
447 $ dlange( '1', q, q, work, ldx, rwork ) / dble( 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 dlacpy( 'Full', m, q, x, ldx, xf, ldx )
456*
457* Compute the CSD
458*
459 CALL dorcsd2by1( '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 dgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
466 $ x, ldx, v1t, ldv1t, zero, work, ldx )
467*
468 CALL dgemm( '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 dgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
480 $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
481*
482 CALL dgemm( '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 = dlange( '1', p, q, x, ldx, rwork )
497 result( 10 ) = ( resid / dble(max(1,p,q)) ) / eps2
498*
499* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
500*
501 resid = dlange( '1', m-p, q, x(p+1,1), ldx, rwork )
502 result( 11 ) = ( resid / dble(max(1,m-p,q)) ) / eps2
503*
504* Compute I - U1'*U1
505*
506 CALL dlaset( 'Full', p, p, zero, one, work, ldu1 )
507 CALL dsyrk( '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 = dlansy( '1', 'Upper', p, work, ldu1, rwork )
513 result( 12 ) = ( resid / dble(max(1,p)) ) / ulp
514*
515* Compute I - U2'*U2
516*
517 CALL dlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
518 CALL dsyrk( '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 = dlansy( '1', 'Upper', m-p, work, ldu2, rwork )
524 result( 13 ) = ( resid / dble(max(1,m-p)) ) / ulp
525*
526* Compute I - V1T*V1T'
527*
528 CALL dlaset( 'Full', q, q, zero, one, work, ldv1t )
529 CALL dsyrk( '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 = dlansy( '1', 'Upper', q, work, ldv1t, rwork )
535 result( 14 ) = ( resid / dble(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 DCSDTS
554*
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dsyrk
subroutine dsyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
DSYRK
Definition dsyrk.f:169
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlange
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlange.f:114
dlansy
double precision function dlansy(norm, uplo, n, a, lda, work)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansy.f:122
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
dorcsd2by1
subroutine dorcsd2by1(jobu1, jobu2, jobv1t, m, p, q, x11, ldx11, x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work, lwork, iwork, info)
DORCSD2BY1
Definition dorcsd2by1.f:233
dorcsd
recursive subroutine dorcsd(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)
DORCSD
Definition dorcsd.f:300
Here is the call graph for this function:
Here is the caller graph for this function: