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

◆ sorcsd()

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

SORCSD

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

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