 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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

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:
 Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.

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 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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
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
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 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
Here is the call graph for this function:
Here is the caller graph for this function: