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

◆ dorcsd()

recursive subroutine dorcsd ( character  jobu1,
character  jobu2,
character  jobv1t,
character  jobv2t,
character  trans,
character  signs,
integer  m,
integer  p,
integer  q,
double precision, dimension( ldx11, * )  x11,
integer  ldx11,
double precision, dimension( ldx12, * )  x12,
integer  ldx12,
double precision, dimension( ldx21, * )  x21,
integer  ldx21,
double precision, dimension( ldx22, * )  x22,
integer  ldx22,
double precision, dimension( * )  theta,
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( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
integer  info 
)

DORCSD

Download DORCSD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DORCSD computes the CS decomposition of an M-by-M partitioned
 orthogonal matrix X:

                                 [  I  0  0 |  0  0  0 ]
                                 [  0  C  0 |  0 -S  0 ]
     [ X11 | X12 ]   [ U1 |    ] [  0  0  0 |  0  0 -I ] [ V1 |    ]**T
 X = [-----------] = [---------] [---------------------] [---------]   .
     [ X21 | X22 ]   [    | U2 ] [  0  0  0 |  I  0  0 ] [    | V2 ]
                                 [  0  S  0 |  0  C  0 ]
                                 [  0  0  I |  0  0  0 ]

 X11 is P-by-Q. The orthogonal matrices U1, U2, V1, and V2 are P-by-P,
 (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
 R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
 which R = MIN(P,M-P,Q,M-Q).
Parameters
[in]JOBU1
          JOBU1 is CHARACTER
          = 'Y':      U1 is computed;
          otherwise:  U1 is not computed.
[in]JOBU2
          JOBU2 is CHARACTER
          = 'Y':      U2 is computed;
          otherwise:  U2 is not computed.
[in]JOBV1T
          JOBV1T is CHARACTER
          = 'Y':      V1T is computed;
          otherwise:  V1T is not computed.
[in]JOBV2T
          JOBV2T is CHARACTER
          = 'Y':      V2T is computed;
          otherwise:  V2T is not computed.
[in]TRANS
          TRANS is CHARACTER
          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
                      order;
          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
                      major order.
[in]SIGNS
          SIGNS is CHARACTER
          = 'O':      The lower-left block is made nonpositive (the
                      "other" convention);
          otherwise:  The upper-right block is made nonpositive (the
                      "default" convention).
[in]M
          M is INTEGER
          The number of rows and columns in X.
[in]P
          P is INTEGER
          The number of rows in X11 and X12. 0 <= P <= M.
[in]Q
          Q is INTEGER
          The number of columns in X11 and X21. 0 <= Q <= M.
[in,out]X11
          X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
          On entry, part of the orthogonal matrix whose CSD is desired.
[in]LDX11
          LDX11 is INTEGER
          The leading dimension of X11. LDX11 >= MAX(1,P).
[in,out]X12
          X12 is DOUBLE PRECISION array, dimension (LDX12,M-Q)
          On entry, part of the orthogonal matrix whose CSD is desired.
[in]LDX12
          LDX12 is INTEGER
          The leading dimension of X12. LDX12 >= MAX(1,P).
[in,out]X21
          X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
          On entry, part of the orthogonal matrix whose CSD is desired.
[in]LDX21
          LDX21 is INTEGER
          The leading dimension of X11. LDX21 >= MAX(1,M-P).
[in,out]X22
          X22 is DOUBLE PRECISION array, dimension (LDX22,M-Q)
          On entry, part of the orthogonal matrix whose CSD is desired.
[in]LDX22
          LDX22 is INTEGER
          The leading dimension of X11. LDX22 >= MAX(1,M-P).
[out]THETA
          THETA is DOUBLE PRECISION array, dimension (R), in which R =
          MIN(P,M-P,Q,M-Q).
          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
[out]U1
          U1 is DOUBLE PRECISION array, dimension (LDU1,P)
          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
[in]LDU1
          LDU1 is INTEGER
          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
          MAX(1,P).
[out]U2
          U2 is DOUBLE PRECISION array, dimension (LDU2,M-P)
          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
          matrix U2.
[in]LDU2
          LDU2 is INTEGER
          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
          MAX(1,M-P).
[out]V1T
          V1T is DOUBLE PRECISION array, dimension (LDV1T,Q)
          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
          matrix V1**T.
[in]LDV1T
          LDV1T is INTEGER
          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
          MAX(1,Q).
[out]V2T
          V2T is DOUBLE PRECISION array, dimension (LDV2T,M-Q)
          If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) orthogonal
          matrix V2**T.
[in]LDV2T
          LDV2T is INTEGER
          The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
          MAX(1,M-Q).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
          If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
          define the matrix in intermediate bidiagonal-block form
          remaining after nonconvergence. INFO specifies the number
          of nonzero PHI's.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the work array, and no error
          message related to LWORK is issued by XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (M-MIN(P, M-P, Q, M-Q))
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  DBBCSD did not converge. See the description of WORK
                above for details.
References:
[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 295 of file dorcsd.f.

300*
301* -- LAPACK computational routine --
302* -- LAPACK is a software package provided by Univ. of Tennessee, --
303* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
304*
305* .. Scalar Arguments ..
306 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
307 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
308 $ LDX21, LDX22, LWORK, M, P, Q
309* ..
310* .. Array Arguments ..
311 INTEGER IWORK( * )
312 DOUBLE PRECISION THETA( * )
313 DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
314 $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
315 $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
316 $ * )
317* ..
318*
319* ===================================================================
320*
321* .. Parameters ..
322 DOUBLE PRECISION ONE, ZERO
323 parameter( one = 1.0d0,
324 $ zero = 0.0d0 )
325* ..
326* .. Local Scalars ..
327 CHARACTER TRANST, SIGNST
328 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
329 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
330 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
331 $ ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN,
332 $ LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN,
333 $ LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN,
334 $ LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN,
335 $ LORGQRWORKOPT, LWORKMIN, LWORKOPT
336 LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
337 $ WANTV1T, WANTV2T
338* ..
339* .. External Subroutines ..
340 EXTERNAL dbbcsd, dlacpy, dlapmr, dlapmt,
342* ..
343* .. External Functions ..
344 LOGICAL LSAME
345 EXTERNAL lsame
346* ..
347* .. Intrinsic Functions
348 INTRINSIC int, max, min
349* ..
350* .. Executable Statements ..
351*
352* Test input arguments
353*
354 info = 0
355 wantu1 = lsame( jobu1, 'Y' )
356 wantu2 = lsame( jobu2, 'Y' )
357 wantv1t = lsame( jobv1t, 'Y' )
358 wantv2t = lsame( jobv2t, 'Y' )
359 colmajor = .NOT. lsame( trans, 'T' )
360 defaultsigns = .NOT. lsame( signs, 'O' )
361 lquery = lwork .EQ. -1
362 IF( m .LT. 0 ) THEN
363 info = -7
364 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
365 info = -8
366 ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
367 info = -9
368 ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
369 info = -11
370 ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
371 info = -11
372 ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
373 info = -13
374 ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
375 info = -13
376 ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
377 info = -15
378 ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
379 info = -15
380 ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
381 info = -17
382 ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
383 info = -17
384 ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
385 info = -20
386 ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
387 info = -22
388 ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
389 info = -24
390 ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
391 info = -26
392 END IF
393*
394* Work with transpose if convenient
395*
396 IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) ) THEN
397 IF( colmajor ) THEN
398 transt = 'T'
399 ELSE
400 transt = 'N'
401 END IF
402 IF( defaultsigns ) THEN
403 signst = 'O'
404 ELSE
405 signst = 'D'
406 END IF
407 CALL dorcsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
408 $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
409 $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
410 $ u2, ldu2, work, lwork, iwork, info )
411 RETURN
412 END IF
413*
414* Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
415* convenient
416*
417 IF( info .EQ. 0 .AND. m-q .LT. q ) THEN
418 IF( defaultsigns ) THEN
419 signst = 'O'
420 ELSE
421 signst = 'D'
422 END IF
423 CALL dorcsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
424 $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
425 $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
426 $ ldv1t, work, lwork, iwork, info )
427 RETURN
428 END IF
429*
430* Compute workspace
431*
432 IF( info .EQ. 0 ) THEN
433*
434 iphi = 2
435 itaup1 = iphi + max( 1, q - 1 )
436 itaup2 = itaup1 + max( 1, p )
437 itauq1 = itaup2 + max( 1, m - p )
438 itauq2 = itauq1 + max( 1, q )
439 iorgqr = itauq2 + max( 1, m - q )
440 CALL dorgqr( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
441 $ childinfo )
442 lorgqrworkopt = int( work(1) )
443 lorgqrworkmin = max( 1, m - q )
444 iorglq = itauq2 + max( 1, m - q )
445 CALL dorglq( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
446 $ childinfo )
447 lorglqworkopt = int( work(1) )
448 lorglqworkmin = max( 1, m - q )
449 iorbdb = itauq2 + max( 1, m - q )
450 CALL dorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
451 $ x21, ldx21, x22, ldx22, theta, v1t, u1, u2, v1t,
452 $ v2t, work, -1, childinfo )
453 lorbdbworkopt = int( work(1) )
454 lorbdbworkmin = lorbdbworkopt
455 ib11d = itauq2 + max( 1, m - q )
456 ib11e = ib11d + max( 1, q )
457 ib12d = ib11e + max( 1, q - 1 )
458 ib12e = ib12d + max( 1, q )
459 ib21d = ib12e + max( 1, q - 1 )
460 ib21e = ib21d + max( 1, q )
461 ib22d = ib21e + max( 1, q - 1 )
462 ib22e = ib22d + max( 1, q )
463 ibbcsd = ib22e + max( 1, q - 1 )
464 CALL dbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
465 $ theta, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
466 $ ldv2t, u1, u1, u1, u1, u1, u1, u1, u1, work, -1,
467 $ childinfo )
468 lbbcsdworkopt = int( work(1) )
469 lbbcsdworkmin = lbbcsdworkopt
470 lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
471 $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkopt ) - 1
472 lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
473 $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkmin ) - 1
474 work(1) = max(lworkopt,lworkmin)
475*
476 IF( lwork .LT. lworkmin .AND. .NOT. lquery ) THEN
477 info = -22
478 ELSE
479 lorgqrwork = lwork - iorgqr + 1
480 lorglqwork = lwork - iorglq + 1
481 lorbdbwork = lwork - iorbdb + 1
482 lbbcsdwork = lwork - ibbcsd + 1
483 END IF
484 END IF
485*
486* Abort if any illegal arguments
487*
488 IF( info .NE. 0 ) THEN
489 CALL xerbla( 'DORCSD', -info )
490 RETURN
491 ELSE IF( lquery ) THEN
492 RETURN
493 END IF
494*
495* Transform to bidiagonal block form
496*
497 CALL dorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
498 $ ldx21, x22, ldx22, theta, work(iphi), work(itaup1),
499 $ work(itaup2), work(itauq1), work(itauq2),
500 $ work(iorbdb), lorbdbwork, childinfo )
501*
502* Accumulate Householder reflectors
503*
504 IF( colmajor ) THEN
505 IF( wantu1 .AND. p .GT. 0 ) THEN
506 CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
507 CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
508 $ lorgqrwork, info)
509 END IF
510 IF( wantu2 .AND. m-p .GT. 0 ) THEN
511 CALL dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
512 CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
513 $ work(iorgqr), lorgqrwork, info )
514 END IF
515 IF( wantv1t .AND. q .GT. 0 ) THEN
516 CALL dlacpy( 'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
517 $ ldv1t )
518 v1t(1, 1) = one
519 DO j = 2, q
520 v1t(1,j) = zero
521 v1t(j,1) = zero
522 END DO
523 CALL dorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
524 $ work(iorglq), lorglqwork, info )
525 END IF
526 IF( wantv2t .AND. m-q .GT. 0 ) THEN
527 CALL dlacpy( 'U', p, m-q, x12, ldx12, v2t, ldv2t )
528 IF (m-p .GT. q) Then
529 CALL dlacpy( 'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
530 $ v2t(p+1,p+1), ldv2t )
531 END IF
532 IF (m .GT. q) THEN
533 CALL dorglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
534 $ work(iorglq), lorglqwork, info )
535 END IF
536 END IF
537 ELSE
538 IF( wantu1 .AND. p .GT. 0 ) THEN
539 CALL dlacpy( 'U', q, p, x11, ldx11, u1, ldu1 )
540 CALL dorglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
541 $ lorglqwork, info)
542 END IF
543 IF( wantu2 .AND. m-p .GT. 0 ) THEN
544 CALL dlacpy( 'U', q, m-p, x21, ldx21, u2, ldu2 )
545 CALL dorglq( m-p, m-p, q, u2, ldu2, work(itaup2),
546 $ work(iorglq), lorglqwork, info )
547 END IF
548 IF( wantv1t .AND. q .GT. 0 ) THEN
549 CALL dlacpy( 'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
550 $ ldv1t )
551 v1t(1, 1) = one
552 DO j = 2, q
553 v1t(1,j) = zero
554 v1t(j,1) = zero
555 END DO
556 CALL dorgqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
557 $ work(iorgqr), lorgqrwork, info )
558 END IF
559 IF( wantv2t .AND. m-q .GT. 0 ) THEN
560 CALL dlacpy( 'L', m-q, p, x12, ldx12, v2t, ldv2t )
561 CALL dlacpy( 'L', m-p-q, m-p-q, x22(p+1,q+1), ldx22,
562 $ v2t(p+1,p+1), ldv2t )
563 CALL dorgqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
564 $ work(iorgqr), lorgqrwork, info )
565 END IF
566 END IF
567*
568* Compute the CSD of the matrix in bidiagonal-block form
569*
570 CALL dbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
571 $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
572 $ ldv2t, work(ib11d), work(ib11e), work(ib12d),
573 $ work(ib12e), work(ib21d), work(ib21e), work(ib22d),
574 $ work(ib22e), work(ibbcsd), lbbcsdwork, info )
575*
576* Permute rows and columns to place identity submatrices in top-
577* left corner of (1,1)-block and/or bottom-right corner of (1,2)-
578* block and/or bottom-right corner of (2,1)-block and/or top-left
579* corner of (2,2)-block
580*
581 IF( q .GT. 0 .AND. wantu2 ) THEN
582 DO i = 1, q
583 iwork(i) = m - p - q + i
584 END DO
585 DO i = q + 1, m - p
586 iwork(i) = i - q
587 END DO
588 IF( colmajor ) THEN
589 CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
590 ELSE
591 CALL dlapmr( .false., m-p, m-p, u2, ldu2, iwork )
592 END IF
593 END IF
594 IF( m .GT. 0 .AND. wantv2t ) THEN
595 DO i = 1, p
596 iwork(i) = m - p - q + i
597 END DO
598 DO i = p + 1, m - q
599 iwork(i) = i - p
600 END DO
601 IF( .NOT. colmajor ) THEN
602 CALL dlapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
603 ELSE
604 CALL dlapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
605 END IF
606 END IF
607*
608 RETURN
609*
610* End DORCSD
611*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dbbcsd(jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta, phi, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t, b11d, b11e, b12d, b12e, b21d, b21e, b22d, b22e, work, lwork, info)
DBBCSD
Definition dbbcsd.f:332
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
subroutine dlapmr(forwrd, m, n, x, ldx, k)
DLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition dlapmr.f:104
subroutine dlapmt(forwrd, m, n, x, ldx, k)
DLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition dlapmt.f:104
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dorbdb(trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21, ldx21, x22, ldx22, theta, phi, taup1, taup2, tauq1, tauq2, work, lwork, info)
DORBDB
Definition dorbdb.f:287
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
subroutine dorglq(m, n, k, a, lda, tau, work, lwork, info)
DORGLQ
Definition dorglq.f:127
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:128
Here is the call graph for this function:
Here is the caller graph for this function: