LAPACK 3.3.1
Linear Algebra PACKage

cuncsd.f

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