LAPACK: SRC/dorcsd.f Source File

Go to the documentation of this file.
00001       RECURSIVE SUBROUTINE DORCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
00002      $                             SIGNS, M, P, Q, X11, LDX11, X12,
00003      $                             LDX12, X21, LDX21, X22, LDX22, THETA,
00004      $                             U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
00005      $                             LDV2T, WORK, LWORK, IWORK, INFO )
00006       IMPLICIT NONE
00007 *
00008 *  -- LAPACK routine (version 3.3.0) --
00009 *
00010 *  -- Contributed by Brian Sutton of the Randolph-Macon College --
00011 *  -- November 2010
00012 *
00013 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00014 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--     
00015 *
00016 *     .. Scalar Arguments ..
00017       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
00018       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
00019      $                   LDX21, LDX22, LWORK, M, P, Q
00020 *     ..
00021 *     .. Array Arguments ..
00022       INTEGER            IWORK( * )
00023       DOUBLE PRECISION   THETA( * )
00024       DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
00025      $                   V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
00026      $                   X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
00027      $                   * )
00028 *     ..
00029 *
00030 *  Purpose
00031 *  =======
00032 *
00033 *  DORCSD computes the CS decomposition of an M-by-M partitioned
00034 *  orthogonal matrix X:
00035 *
00036 *                                  [  I  0  0 |  0  0  0 ]
00037 *                                  [  0  C  0 |  0 -S  0 ]
00038 *      [ X11 | X12 ]   [ U1 |    ] [  0  0  0 |  0  0 -I ] [ V1 |    ]**T
00039 *  X = [-----------] = [---------] [---------------------] [---------]   .
00040 *      [ X21 | X22 ]   [    | U2 ] [  0  0  0 |  I  0  0 ] [    | V2 ]
00041 *                                  [  0  S  0 |  0  C  0 ]
00042 *                                  [  0  0  I |  0  0  0 ]
00043 *
00044 *  X11 is P-by-Q. The orthogonal matrices U1, U2, V1, and V2 are P-by-P,
00045 *  (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
00046 *  R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
00047 *  which R = MIN(P,M-P,Q,M-Q).
00048 *
00049 *  Arguments
00050 *  =========
00051 *
00052 *  JOBU1   (input) CHARACTER
00053 *          = 'Y':      U1 is computed;
00054 *          otherwise:  U1 is not computed.
00055 *
00056 *  JOBU2   (input) CHARACTER
00057 *          = 'Y':      U2 is computed;
00058 *          otherwise:  U2 is not computed.
00059 *
00060 *  JOBV1T  (input) CHARACTER
00061 *          = 'Y':      V1T is computed;
00062 *          otherwise:  V1T is not computed.
00063 *
00064 *  JOBV2T  (input) CHARACTER
00065 *          = 'Y':      V2T is computed;
00066 *          otherwise:  V2T is not computed.
00067 *
00068 *  TRANS   (input) CHARACTER
00069 *          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
00070 *                      order;
00071 *          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
00072 *                      major order.
00073 *
00074 *  SIGNS   (input) CHARACTER
00075 *          = 'O':      The lower-left block is made nonpositive (the
00076 *                      "other" convention);
00077 *          otherwise:  The upper-right block is made nonpositive (the
00078 *                      "default" convention).
00079 *
00080 *  M       (input) INTEGER
00081 *          The number of rows and columns in X.
00082 *
00083 *  P       (input) INTEGER
00084 *          The number of rows in X11 and X12. 0 <= P <= M.
00085 *
00086 *  Q       (input) INTEGER
00087 *          The number of columns in X11 and X21. 0 <= Q <= M.
00088 *
00089 *  X       (input/workspace) DOUBLE PRECISION array, dimension (LDX,M)
00090 *          On entry, the orthogonal matrix whose CSD is desired.
00091 *
00092 *  LDX     (input) INTEGER
00093 *          The leading dimension of X. LDX >= MAX(1,M).
00094 *
00095 *  THETA   (output) DOUBLE PRECISION array, dimension (R), in which R =
00096 *          MIN(P,M-P,Q,M-Q).
00097 *          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
00098 *          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
00099 *
00100 *  U1      (output) DOUBLE PRECISION array, dimension (P)
00101 *          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
00102 *
00103 *  LDU1    (input) INTEGER
00104 *          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
00105 *          MAX(1,P).
00106 *
00107 *  U2      (output) DOUBLE PRECISION array, dimension (M-P)
00108 *          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
00109 *          matrix U2.
00110 *
00111 *  LDU2    (input) INTEGER
00112 *          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
00113 *          MAX(1,M-P).
00114 *
00115 *  V1T     (output) DOUBLE PRECISION array, dimension (Q)
00116 *          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
00117 *          matrix V1**T.
00118 *
00119 *  LDV1T   (input) INTEGER
00120 *          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
00121 *          MAX(1,Q).
00122 *
00123 *  V2T     (output) DOUBLE PRECISION array, dimension (M-Q)
00124 *          If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) orthogonal
00125 *          matrix V2**T.
00126 *
00127 *  LDV2T   (input) INTEGER
00128 *          The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
00129 *          MAX(1,M-Q).
00130 *
00131 *  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00132 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00133 *          If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
00134 *          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
00135 *          define the matrix in intermediate bidiagonal-block form
00136 *          remaining after nonconvergence. INFO specifies the number
00137 *          of nonzero PHI's.
00138 *
00139 *  LWORK   (input) INTEGER
00140 *          The dimension of the array WORK.
00141 *
00142 *          If LWORK = -1, then a workspace query is assumed; the routine
00143 *          only calculates the optimal size of the WORK array, returns
00144 *          this value as the first entry of the work array, and no error
00145 *          message related to LWORK is issued by XERBLA.
00146 *
00147 *  IWORK   (workspace) INTEGER array, dimension (M-Q)
00148 *
00149 *  INFO    (output) INTEGER
00150 *          = 0:  successful exit.
00151 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00152 *          > 0:  DBBCSD did not converge. See the description of WORK
00153 *                above for details.
00154 *
00155 *  Reference
00156 *  =========
00157 *
00158 *  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
00159 *      Algorithms, 50(1):33-65, 2009.
00160 *
00161 *  ===================================================================
00162 *
00163 *     .. Parameters ..
00164       DOUBLE PRECISION   REALONE
00165       PARAMETER          ( REALONE = 1.0D0 )
00166       DOUBLE PRECISION   NEGONE, ONE, PIOVER2, ZERO
00167       PARAMETER          ( NEGONE = -1.0D0, ONE = 1.0D0,
00168      $                     PIOVER2 = 1.57079632679489662D0,
00169      $                     ZERO = 0.0D0 )
00170 *     ..
00171 *     .. Local Scalars ..
00172       CHARACTER          TRANST, SIGNST
00173       INTEGER            CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
00174      $                   IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
00175      $                   IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
00176      $                   ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN,
00177      $                   LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN,
00178      $                   LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN,
00179      $                   LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN,
00180      $                   LORGQRWORKOPT, LWORKMIN, LWORKOPT
00181       LOGICAL            COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
00182      $                   WANTV1T, WANTV2T
00183 *     ..
00184 *     .. External Subroutines ..
00185       EXTERNAL           DBBCSD, DLACPY, DLAPMR, DLAPMT, DLASCL, DLASET,
00186      $                   DORBDB, DORGLQ, DORGQR, XERBLA
00187 *     ..
00188 *     .. External Functions ..
00189       LOGICAL            LSAME
00190       EXTERNAL           LSAME
00191 *     ..
00192 *     .. Intrinsic Functions
00193       INTRINSIC          COS, INT, MAX, MIN, SIN
00194 *     ..
00195 *     .. Executable Statements ..
00196 *
00197 *     Test input arguments
00198 *
00199       INFO = 0
00200       WANTU1 = LSAME( JOBU1, 'Y' )
00201       WANTU2 = LSAME( JOBU2, 'Y' )
00202       WANTV1T = LSAME( JOBV1T, 'Y' )
00203       WANTV2T = LSAME( JOBV2T, 'Y' )
00204       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
00205       DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' )
00206       LQUERY = LWORK .EQ. -1
00207       IF( M .LT. 0 ) THEN
00208          INFO = -7
00209       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
00210          INFO = -8
00211       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
00212          INFO = -9
00213       ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR.
00214      $         ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN
00215          INFO = -11
00216       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
00217          INFO = -14
00218       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
00219          INFO = -16
00220       ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
00221          INFO = -18
00222       ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
00223          INFO = -20
00224       END IF
00225 *
00226 *     Work with transpose if convenient
00227 *
00228       IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN
00229          IF( COLMAJOR ) THEN
00230             TRANST = 'T'
00231          ELSE
00232             TRANST = 'N'
00233          END IF
00234          IF( DEFAULTSIGNS ) THEN
00235             SIGNST = 'O'
00236          ELSE
00237             SIGNST = 'D'
00238          END IF
00239          CALL DORCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M,
00240      $                Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22,
00241      $                LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1,
00242      $                U2, LDU2, WORK, LWORK, IWORK, INFO )
00243          RETURN
00244       END IF
00245 *
00246 *     Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
00247 *     convenient
00248 *
00249       IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN
00250          IF( DEFAULTSIGNS ) THEN
00251             SIGNST = 'O'
00252          ELSE
00253             SIGNST = 'D'
00254          END IF
00255          CALL DORCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M,
00256      $                M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11,
00257      $                LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T,
00258      $                LDV1T, WORK, LWORK, IWORK, INFO )
00259          RETURN
00260       END IF
00261 *
00262 *     Compute workspace
00263 *
00264       IF( INFO .EQ. 0 ) THEN
00265 *
00266          IPHI = 2
00267          ITAUP1 = IPHI + MAX( 1, Q - 1 )
00268          ITAUP2 = ITAUP1 + MAX( 1, P )
00269          ITAUQ1 = ITAUP2 + MAX( 1, M - P )
00270          ITAUQ2 = ITAUQ1 + MAX( 1, Q )
00271          IORGQR = ITAUQ2 + MAX( 1, M - Q )
00272          CALL DORGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
00273      $                CHILDINFO )
00274          LORGQRWORKOPT = INT( WORK(1) )
00275          LORGQRWORKMIN = MAX( 1, M - Q )
00276          IORGLQ = ITAUQ2 + MAX( 1, M - Q )
00277          CALL DORGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
00278      $                CHILDINFO )
00279          LORGLQWORKOPT = INT( WORK(1) )
00280          LORGLQWORKMIN = MAX( 1, M - Q )
00281          IORBDB = ITAUQ2 + MAX( 1, M - Q )
00282          CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
00283      $                X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK,
00284      $                -1, CHILDINFO )
00285          LORBDBWORKOPT = INT( WORK(1) )
00286          LORBDBWORKMIN = LORBDBWORKOPT
00287          IB11D = ITAUQ2 + MAX( 1, M - Q )
00288          IB11E = IB11D + MAX( 1, Q )
00289          IB12D = IB11E + MAX( 1, Q - 1 )
00290          IB12E = IB12D + MAX( 1, Q )
00291          IB21D = IB12E + MAX( 1, Q - 1 )
00292          IB21E = IB21D + MAX( 1, Q )
00293          IB22D = IB21E + MAX( 1, Q - 1 )
00294          IB22E = IB22D + MAX( 1, Q )
00295          IBBCSD = IB22E + MAX( 1, Q - 1 )
00296          CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0,
00297      $                0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0,
00298      $                0, 0, 0, 0, 0, 0, 0, WORK, -1, CHILDINFO )
00299          LBBCSDWORKOPT = INT( WORK(1) )
00300          LBBCSDWORKMIN = LBBCSDWORKOPT
00301          LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT,
00302      $              IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKOPT ) - 1
00303          LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN,
00304      $              IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKMIN ) - 1
00305          WORK(1) = LWORKOPT
00306 *
00307          IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
00308             INFO = -22
00309          ELSE
00310             LORGQRWORK = LWORK - IORGQR + 1
00311             LORGLQWORK = LWORK - IORGLQ + 1
00312             LORBDBWORK = LWORK - IORBDB + 1
00313             LBBCSDWORK = LWORK - IBBCSD + 1
00314          END IF
00315       END IF
00316 *
00317 *     Abort if any illegal arguments
00318 *
00319       IF( INFO .NE. 0 ) THEN
00320          CALL XERBLA( 'DORCSD', -INFO )
00321          RETURN
00322       ELSE IF( LQUERY ) THEN
00323          RETURN
00324       END IF
00325 *
00326 *     Transform to bidiagonal block form
00327 *
00328       CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21,
00329      $             LDX21, X22, LDX22, THETA, WORK(IPHI), WORK(ITAUP1),
00330      $             WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2),
00331      $             WORK(IORBDB), LORBDBWORK, CHILDINFO )
00332 *
00333 *     Accumulate Householder reflectors
00334 *
00335       IF( COLMAJOR ) THEN
00336          IF( WANTU1 .AND. P .GT. 0 ) THEN
00337             CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
00338             CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
00339      $                   LORGQRWORK, INFO)
00340          END IF
00341          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
00342             CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
00343             CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
00344      $                   WORK(IORGQR), LORGQRWORK, INFO )
00345          END IF
00346          IF( WANTV1T .AND. Q .GT. 0 ) THEN
00347             CALL DLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2),
00348      $                   LDV1T )
00349             V1T(1, 1) = ONE
00350             DO J = 2, Q
00351                V1T(1,J) = ZERO
00352                V1T(J,1) = ZERO
00353             END DO
00354             CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
00355      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00356          END IF
00357          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
00358             CALL DLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T )
00359             CALL DLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
00360      $                   V2T(P+1,P+1), LDV2T )
00361             CALL DORGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
00362      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00363          END IF
00364       ELSE
00365          IF( WANTU1 .AND. P .GT. 0 ) THEN
00366             CALL DLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 )
00367             CALL DORGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ),
00368      $                   LORGLQWORK, INFO)
00369          END IF
00370          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
00371             CALL DLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 )
00372             CALL DORGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
00373      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00374          END IF
00375          IF( WANTV1T .AND. Q .GT. 0 ) THEN
00376             CALL DLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2),
00377      $                   LDV1T )
00378             V1T(1, 1) = ONE
00379             DO J = 2, Q
00380                V1T(1,J) = ZERO
00381                V1T(J,1) = ZERO
00382             END DO
00383             CALL DORGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
00384      $                   WORK(IORGQR), LORGQRWORK, INFO )
00385          END IF
00386          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
00387             CALL DLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T )
00388             CALL DLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22,
00389      $                   V2T(P+1,P+1), LDV2T )
00390             CALL DORGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
00391      $                   WORK(IORGQR), LORGQRWORK, INFO )
00392          END IF
00393       END IF
00394 *
00395 *     Compute the CSD of the matrix in bidiagonal-block form
00396 *
00397       CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA,
00398      $             WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
00399      $             LDV2T, WORK(IB11D), WORK(IB11E), WORK(IB12D),
00400      $             WORK(IB12E), WORK(IB21D), WORK(IB21E), WORK(IB22D),
00401      $             WORK(IB22E), WORK(IBBCSD), LBBCSDWORK, INFO )
00402 *
00403 *     Permute rows and columns to place identity submatrices in top-
00404 *     left corner of (1,1)-block and/or bottom-right corner of (1,2)-
00405 *     block and/or bottom-right corner of (2,1)-block and/or top-left
00406 *     corner of (2,2)-block 
00407 *
00408       IF( Q .GT. 0 .AND. WANTU2 ) THEN
00409          DO I = 1, Q
00410             IWORK(I) = M - P - Q + I
00411          END DO
00412          DO I = Q + 1, M - P
00413             IWORK(I) = I - Q
00414          END DO
00415          IF( COLMAJOR ) THEN
00416             CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
00417          ELSE
00418             CALL DLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK )
00419          END IF
00420       END IF
00421       IF( M .GT. 0 .AND. WANTV2T ) THEN
00422          DO I = 1, P
00423             IWORK(I) = M - P - Q + I
00424          END DO
00425          DO I = P + 1, M - Q
00426             IWORK(I) = I - P
00427          END DO
00428          IF( .NOT. COLMAJOR ) THEN
00429             CALL DLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
00430          ELSE
00431             CALL DLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
00432          END IF
00433       END IF
00434 *
00435       RETURN
00436 *
00437 *     End DORCSD
00438 *
00439       END
00440