LAPACK 3.3.1
Linear Algebra PACKage

dorcsd.f

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