LAPACK: SRC/cuncsd.f Source File

Go to the documentation of this file.
00001       RECURSIVE SUBROUTINE CUNCSD( 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, RWORK, LRWORK,
00006      $                             IWORK, INFO )
00007       IMPLICIT NONE
00008 *
00009 *  -- LAPACK routine (version 3.3.0) --
00010 *
00011 *  -- Contributed by Brian Sutton of the Randolph-Macon College --
00012 *  -- November 2010
00013 *
00014 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00015 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--     
00016 *
00017 *     .. Scalar Arguments ..
00018       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
00019       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
00020      $                   LDX21, LDX22, LRWORK, LWORK, M, P, Q
00021 *     ..
00022 *     .. Array Arguments ..
00023       INTEGER            IWORK( * )
00024       REAL               THETA( * )
00025       REAL               RWORK( * )
00026       COMPLEX            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 *  CUNCSD computes the CS decomposition of an M-by-M partitioned
00036 *  unitary 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 |    ]**H
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 unitary 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) COMPLEX array, dimension (LDX,M)
00092 *          On entry, the unitary matrix whose CSD is desired.
00093 *
00094 *  LDX     (input) INTEGER
00095 *          The leading dimension of X. LDX >= MAX(1,M).
00096 *
00097 *  THETA   (output) REAL 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) COMPLEX array, dimension (P)
00103 *          If JOBU1 = 'Y', U1 contains the P-by-P unitary 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) COMPLEX array, dimension (M-P)
00110 *          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
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) COMPLEX array, dimension (Q)
00118 *          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
00119 *          matrix V1**H.
00120 *
00121 *  LDV1T   (input) INTEGER
00122 *          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
00123 *          MAX(1,Q).
00124 *
00125 *  V2T     (output) COMPLEX array, dimension (M-Q)
00126 *          If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) unitary
00127 *          matrix V2**H.
00128 *
00129 *  LDV2T   (input) INTEGER
00130 *          The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
00131 *          MAX(1,M-Q).
00132 *
00133 *  WORK    (workspace) COMPLEX array, dimension (MAX(1,LWORK))
00134 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00135 *
00136 *  LWORK   (input) INTEGER
00137 *          The dimension of the array WORK.
00138 *
00139 *          If LWORK = -1, then a workspace query is assumed; the routine
00140 *          only calculates the optimal size of the WORK array, returns
00141 *          this value as the first entry of the work array, and no error
00142 *          message related to LWORK is issued by XERBLA.
00143 *
00144 *  RWORK   (workspace) REAL array, dimension MAX(1,LRWORK)
00145 *          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
00146 *          If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
00147 *          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
00148 *          define the matrix in intermediate bidiagonal-block form
00149 *          remaining after nonconvergence. INFO specifies the number
00150 *          of nonzero PHI's.
00151 *
00152 *  LRWORK  (input) INTEGER
00153 *          The dimension of the array RWORK.
00154 *
00155 *          If LRWORK = -1, then a workspace query is assumed; the routine
00156 *          only calculates the optimal size of the RWORK array, returns
00157 *          this value as the first entry of the work array, and no error
00158 *          message related to LRWORK is issued by XERBLA.
00159 *
00160 *  IWORK   (workspace) INTEGER array, dimension (M-Q)
00161 *
00162 *  INFO    (output) INTEGER
00163 *          = 0:  successful exit.
00164 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00165 *          > 0:  CBBCSD did not converge. See the description of RWORK
00166 *                above for details.
00167 *
00168 *  Reference
00169 *  =========
00170 *
00171 *  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
00172 *      Algorithms, 50(1):33-65, 2009.
00173 *
00174 *  ===================================================================
00175 *
00176 *     .. Parameters ..
00177       REAL               REALONE
00178       PARAMETER          ( REALONE = 1.0E0 )
00179       COMPLEX            NEGONE, ONE, PIOVER2, ZERO
00180       PARAMETER          ( NEGONE = (-1.0E0,0.0E0), ONE = (1.0E0,0.0E0),
00181      $                     PIOVER2 = 1.57079632679489662E0,
00182      $                     ZERO = (0.0E0,0.0E0) )
00183 *     ..
00184 *     .. Local Scalars ..
00185       CHARACTER          TRANST, SIGNST
00186       INTEGER            CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
00187      $                   IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
00188      $                   IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
00189      $                   ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN,
00190      $                   LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN,
00191      $                   LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN,
00192      $                   LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN,
00193      $                   LORGQRWORKOPT, LWORKMIN, LWORKOPT
00194       LOGICAL            COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
00195      $                   WANTV1T, WANTV2T
00196       INTEGER            LRWORKMIN, LRWORKOPT
00197       LOGICAL            LRQUERY
00198 *     ..
00199 *     .. External Subroutines ..
00200       EXTERNAL           CBBCSD, CLACPY, CLAPMR, CLAPMT, CLASCL, CLASET,
00201      $                   CUNBDB, CUNGLQ, CUNGQR, XERBLA
00202 *     ..
00203 *     .. External Functions ..
00204       LOGICAL            LSAME
00205       EXTERNAL           LSAME
00206 *     ..
00207 *     .. Intrinsic Functions
00208       INTRINSIC          COS, INT, MAX, MIN, SIN
00209 *     ..
00210 *     .. Executable Statements ..
00211 *
00212 *     Test input arguments
00213 *
00214       INFO = 0
00215       WANTU1 = LSAME( JOBU1, 'Y' )
00216       WANTU2 = LSAME( JOBU2, 'Y' )
00217       WANTV1T = LSAME( JOBV1T, 'Y' )
00218       WANTV2T = LSAME( JOBV2T, 'Y' )
00219       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
00220       DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' )
00221       LQUERY = LWORK .EQ. -1
00222       LRQUERY = LRWORK .EQ. -1
00223       IF( M .LT. 0 ) THEN
00224          INFO = -7
00225       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
00226          INFO = -8
00227       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
00228          INFO = -9
00229       ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR.
00230      $         ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN
00231          INFO = -11
00232       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
00233          INFO = -14
00234       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
00235          INFO = -16
00236       ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
00237          INFO = -18
00238       ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
00239          INFO = -20
00240       END IF
00241 *
00242 *     Work with transpose if convenient
00243 *
00244       IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN
00245          IF( COLMAJOR ) THEN
00246             TRANST = 'T'
00247          ELSE
00248             TRANST = 'N'
00249          END IF
00250          IF( DEFAULTSIGNS ) THEN
00251             SIGNST = 'O'
00252          ELSE
00253             SIGNST = 'D'
00254          END IF
00255          CALL CUNCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M,
00256      $                Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22,
00257      $                LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1,
00258      $                U2, LDU2, WORK, LWORK, RWORK, LRWORK, IWORK,
00259      $                INFO )
00260          RETURN
00261       END IF
00262 *
00263 *     Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
00264 *     convenient
00265 *
00266       IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN
00267          IF( DEFAULTSIGNS ) THEN
00268             SIGNST = 'O'
00269          ELSE
00270             SIGNST = 'D'
00271          END IF
00272          CALL CUNCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M,
00273      $                M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11,
00274      $                LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T,
00275      $                LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO )
00276          RETURN
00277       END IF
00278 *
00279 *     Compute workspace
00280 *
00281       IF( INFO .EQ. 0 ) THEN
00282 *
00283 *        Real workspace
00284 *
00285          IPHI = 2
00286          IB11D = IPHI + MAX( 1, Q - 1 )
00287          IB11E = IB11D + MAX( 1, Q )
00288          IB12D = IB11E + MAX( 1, Q - 1 )
00289          IB12E = IB12D + MAX( 1, Q )
00290          IB21D = IB12E + MAX( 1, Q - 1 )
00291          IB21E = IB21D + MAX( 1, Q )
00292          IB22D = IB21E + MAX( 1, Q - 1 )
00293          IB22E = IB22D + MAX( 1, Q )
00294          IBBCSD = IB22E + MAX( 1, Q - 1 )
00295          CALL CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0,
00296      $                0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0,
00297      $                0, 0, 0, 0, 0, 0, 0, RWORK, -1, CHILDINFO )
00298          LBBCSDWORKOPT = INT( RWORK(1) )
00299          LBBCSDWORKMIN = LBBCSDWORKOPT
00300          LRWORKOPT = IBBCSD + LBBCSDWORKOPT - 1
00301          LRWORKMIN = IBBCSD + LBBCSDWORKMIN - 1
00302          RWORK(1) = LRWORKOPT
00303 *
00304 *        Complex workspace
00305 *
00306          ITAUP1 = 2
00307          ITAUP2 = ITAUP1 + MAX( 1, P )
00308          ITAUQ1 = ITAUP2 + MAX( 1, M - P )
00309          ITAUQ2 = ITAUQ1 + MAX( 1, Q )
00310          IORGQR = ITAUQ2 + MAX( 1, M - Q )
00311          CALL CUNGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
00312      $                CHILDINFO )
00313          LORGQRWORKOPT = INT( WORK(1) )
00314          LORGQRWORKMIN = MAX( 1, M - Q )
00315          IORGLQ = ITAUQ2 + MAX( 1, M - Q )
00316          CALL CUNGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
00317      $                CHILDINFO )
00318          LORGLQWORKOPT = INT( WORK(1) )
00319          LORGLQWORKMIN = MAX( 1, M - Q )
00320          IORBDB = ITAUQ2 + MAX( 1, M - Q )
00321          CALL CUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
00322      $                X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK,
00323      $                -1, CHILDINFO )
00324          LORBDBWORKOPT = INT( WORK(1) )
00325          LORBDBWORKMIN = LORBDBWORKOPT
00326          LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT,
00327      $              IORBDB + LORBDBWORKOPT ) - 1
00328          LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN,
00329      $              IORBDB + LORBDBWORKMIN ) - 1
00330          WORK(1) = LWORKOPT
00331 *
00332          IF( LWORK .LT. LWORKMIN
00333      $       .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN
00334             INFO = -22
00335          ELSE IF( LRWORK .LT. LRWORKMIN
00336      $            .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN
00337             INFO = -24
00338          ELSE
00339             LORGQRWORK = LWORK - IORGQR + 1
00340             LORGLQWORK = LWORK - IORGLQ + 1
00341             LORBDBWORK = LWORK - IORBDB + 1
00342             LBBCSDWORK = LRWORK - IBBCSD + 1
00343          END IF
00344       END IF
00345 *
00346 *     Abort if any illegal arguments
00347 *
00348       IF( INFO .NE. 0 ) THEN
00349          CALL XERBLA( 'CUNCSD', -INFO )
00350          RETURN
00351       ELSE IF( LQUERY .OR. LRQUERY ) THEN
00352          RETURN
00353       END IF
00354 *
00355 *     Transform to bidiagonal block form
00356 *
00357       CALL CUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21,
00358      $             LDX21, X22, LDX22, THETA, RWORK(IPHI), WORK(ITAUP1),
00359      $             WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2),
00360      $             WORK(IORBDB), LORBDBWORK, CHILDINFO )
00361 *
00362 *     Accumulate Householder reflectors
00363 *
00364       IF( COLMAJOR ) THEN
00365          IF( WANTU1 .AND. P .GT. 0 ) THEN
00366             CALL CLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
00367             CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
00368      $                   LORGQRWORK, INFO)
00369          END IF
00370          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
00371             CALL CLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
00372             CALL CUNGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
00373      $                   WORK(IORGQR), LORGQRWORK, INFO )
00374          END IF
00375          IF( WANTV1T .AND. Q .GT. 0 ) THEN
00376             CALL CLACPY( 'U', Q-1, Q-1, X11(1,2), 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 CUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
00384      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00385          END IF
00386          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
00387             CALL CLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T )
00388             CALL CLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
00389      $                   V2T(P+1,P+1), LDV2T )
00390             CALL CUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
00391      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00392          END IF
00393       ELSE
00394          IF( WANTU1 .AND. P .GT. 0 ) THEN
00395             CALL CLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 )
00396             CALL CUNGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ),
00397      $                   LORGLQWORK, INFO)
00398          END IF
00399          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
00400             CALL CLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 )
00401             CALL CUNGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
00402      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00403          END IF
00404          IF( WANTV1T .AND. Q .GT. 0 ) THEN
00405             CALL CLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2),
00406      $                   LDV1T )
00407             V1T(1, 1) = ONE
00408             DO J = 2, Q
00409                V1T(1,J) = ZERO
00410                V1T(J,1) = ZERO
00411             END DO
00412             CALL CUNGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
00413      $                   WORK(IORGQR), LORGQRWORK, INFO )
00414          END IF
00415          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
00416             CALL CLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T )
00417             CALL CLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22,
00418      $                   V2T(P+1,P+1), LDV2T )
00419             CALL CUNGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
00420      $                   WORK(IORGQR), LORGQRWORK, INFO )
00421          END IF
00422       END IF
00423 *
00424 *     Compute the CSD of the matrix in bidiagonal-block form
00425 *
00426       CALL CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA,
00427      $             RWORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
00428      $             LDV2T, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D),
00429      $             RWORK(IB12E), RWORK(IB21D), RWORK(IB21E),
00430      $             RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD),
00431      $             LBBCSDWORK, INFO )
00432 *
00433 *     Permute rows and columns to place identity submatrices in top-
00434 *     left corner of (1,1)-block and/or bottom-right corner of (1,2)-
00435 *     block and/or bottom-right corner of (2,1)-block and/or top-left
00436 *     corner of (2,2)-block 
00437 *
00438       IF( Q .GT. 0 .AND. WANTU2 ) THEN
00439          DO I = 1, Q
00440             IWORK(I) = M - P - Q + I
00441          END DO
00442          DO I = Q + 1, M - P
00443             IWORK(I) = I - Q
00444          END DO
00445          IF( COLMAJOR ) THEN
00446             CALL CLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
00447          ELSE
00448             CALL CLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK )
00449          END IF
00450       END IF
00451       IF( M .GT. 0 .AND. WANTV2T ) THEN
00452          DO I = 1, P
00453             IWORK(I) = M - P - Q + I
00454          END DO
00455          DO I = P + 1, M - Q
00456             IWORK(I) = I - P
00457          END DO
00458          IF( .NOT. COLMAJOR ) THEN
00459             CALL CLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
00460          ELSE
00461             CALL CLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
00462          END IF
00463       END IF
00464 *
00465       RETURN
00466 *
00467 *     End CUNCSD
00468 *
00469       END
00470