LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dorcsd.f
Go to the documentation of this file.
00001 *> \brief \b DORCSD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DORCSD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       RECURSIVE SUBROUTINE DORCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
00022 *                                    SIGNS, M, P, Q, X11, LDX11, X12,
00023 *                                    LDX12, X21, LDX21, X22, LDX22, THETA,
00024 *                                    U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
00025 *                                    LDV2T, WORK, LWORK, IWORK, INFO )
00026 * 
00027 *       .. Scalar Arguments ..
00028 *       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
00029 *       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
00030 *      $                   LDX21, LDX22, LWORK, M, P, Q
00031 *       ..
00032 *       .. Array Arguments ..
00033 *       INTEGER            IWORK( * )
00034 *       DOUBLE PRECISION   THETA( * )
00035 *       DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
00036 *      $                   V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
00037 *      $                   X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
00038 *      $                   * )
00039 *       ..
00040 *  
00041 *
00042 *> \par Purpose:
00043 *  =============
00044 *>
00045 *> \verbatim
00046 *>
00047 *> DORCSD computes the CS decomposition of an M-by-M partitioned
00048 *> orthogonal matrix X:
00049 *>
00050 *>                                 [  I  0  0 |  0  0  0 ]
00051 *>                                 [  0  C  0 |  0 -S  0 ]
00052 *>     [ X11 | X12 ]   [ U1 |    ] [  0  0  0 |  0  0 -I ] [ V1 |    ]**T
00053 *> X = [-----------] = [---------] [---------------------] [---------]   .
00054 *>     [ X21 | X22 ]   [    | U2 ] [  0  0  0 |  I  0  0 ] [    | V2 ]
00055 *>                                 [  0  S  0 |  0  C  0 ]
00056 *>                                 [  0  0  I |  0  0  0 ]
00057 *>
00058 *> X11 is P-by-Q. The orthogonal matrices U1, U2, V1, and V2 are P-by-P,
00059 *> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
00060 *> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
00061 *> which R = MIN(P,M-P,Q,M-Q).
00062 *> \endverbatim
00063 *
00064 *  Arguments:
00065 *  ==========
00066 *
00067 *> \param[in] JOBU1
00068 *> \verbatim
00069 *>          JOBU1 is CHARACTER
00070 *>          = 'Y':      U1 is computed;
00071 *>          otherwise:  U1 is not computed.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] JOBU2
00075 *> \verbatim
00076 *>          JOBU2 is CHARACTER
00077 *>          = 'Y':      U2 is computed;
00078 *>          otherwise:  U2 is not computed.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] JOBV1T
00082 *> \verbatim
00083 *>          JOBV1T is CHARACTER
00084 *>          = 'Y':      V1T is computed;
00085 *>          otherwise:  V1T is not computed.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] JOBV2T
00089 *> \verbatim
00090 *>          JOBV2T is CHARACTER
00091 *>          = 'Y':      V2T is computed;
00092 *>          otherwise:  V2T is not computed.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] TRANS
00096 *> \verbatim
00097 *>          TRANS is CHARACTER
00098 *>          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
00099 *>                      order;
00100 *>          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
00101 *>                      major order.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] SIGNS
00105 *> \verbatim
00106 *>          SIGNS is CHARACTER
00107 *>          = 'O':      The lower-left block is made nonpositive (the
00108 *>                      "other" convention);
00109 *>          otherwise:  The upper-right block is made nonpositive (the
00110 *>                      "default" convention).
00111 *> \endverbatim
00112 *>
00113 *> \param[in] M
00114 *> \verbatim
00115 *>          M is INTEGER
00116 *>          The number of rows and columns in X.
00117 *> \endverbatim
00118 *>
00119 *> \param[in] P
00120 *> \verbatim
00121 *>          P is INTEGER
00122 *>          The number of rows in X11 and X12. 0 <= P <= M.
00123 *> \endverbatim
00124 *>
00125 *> \param[in] Q
00126 *> \verbatim
00127 *>          Q is INTEGER
00128 *>          The number of columns in X11 and X21. 0 <= Q <= M.
00129 *> \endverbatim
00130 *>
00131 *> \param[in,out] X11
00132 *> \verbatim
00133 *>          X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
00134 *>          On entry, part of the orthogonal matrix whose CSD is desired.
00135 *> \endverbatim
00136 *>
00137 *> \param[in] LDX11
00138 *> \verbatim
00139 *>          LDX11 is INTEGER
00140 *>          The leading dimension of X11. LDX11 >= MAX(1,P).
00141 *> \endverbatim
00142 *>
00143 *> \param[in,out] X12
00144 *> \verbatim
00145 *>          X12 is DOUBLE PRECISION array, dimension (LDX12,M-Q)
00146 *>          On entry, part of the orthogonal matrix whose CSD is desired.
00147 *> \endverbatim
00148 *>
00149 *> \param[in] LDX12
00150 *> \verbatim
00151 *>          LDX12 is INTEGER
00152 *>          The leading dimension of X12. LDX12 >= MAX(1,P).
00153 *> \endverbatim
00154 *>
00155 *> \param[in,out] X21
00156 *> \verbatim
00157 *>          X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
00158 *>          On entry, part of the orthogonal matrix whose CSD is desired.
00159 *> \endverbatim
00160 *>
00161 *> \param[in] LDX21
00162 *> \verbatim
00163 *>          LDX21 is INTEGER
00164 *>          The leading dimension of X11. LDX21 >= MAX(1,M-P).
00165 *> \endverbatim
00166 *>
00167 *> \param[in,out] X22
00168 *> \verbatim
00169 *>          X22 is DOUBLE PRECISION array, dimension (LDX22,M-Q)
00170 *>          On entry, part of the orthogonal matrix whose CSD is desired.
00171 *> \endverbatim
00172 *>
00173 *> \param[in] LDX22
00174 *> \verbatim
00175 *>          LDX22 is INTEGER
00176 *>          The leading dimension of X11. LDX22 >= MAX(1,M-P).
00177 *> \endverbatim
00178 *>
00179 *> \param[out] THETA
00180 *> \verbatim
00181 *>          THETA is DOUBLE PRECISION array, dimension (R), in which R =
00182 *>          MIN(P,M-P,Q,M-Q).
00183 *>          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
00184 *>          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
00185 *> \endverbatim
00186 *>
00187 *> \param[out] U1
00188 *> \verbatim
00189 *>          U1 is DOUBLE PRECISION array, dimension (P)
00190 *>          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
00191 *> \endverbatim
00192 *>
00193 *> \param[in] LDU1
00194 *> \verbatim
00195 *>          LDU1 is INTEGER
00196 *>          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
00197 *>          MAX(1,P).
00198 *> \endverbatim
00199 *>
00200 *> \param[out] U2
00201 *> \verbatim
00202 *>          U2 is DOUBLE PRECISION array, dimension (M-P)
00203 *>          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
00204 *>          matrix U2.
00205 *> \endverbatim
00206 *>
00207 *> \param[in] LDU2
00208 *> \verbatim
00209 *>          LDU2 is INTEGER
00210 *>          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
00211 *>          MAX(1,M-P).
00212 *> \endverbatim
00213 *>
00214 *> \param[out] V1T
00215 *> \verbatim
00216 *>          V1T is DOUBLE PRECISION array, dimension (Q)
00217 *>          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
00218 *>          matrix V1**T.
00219 *> \endverbatim
00220 *>
00221 *> \param[in] LDV1T
00222 *> \verbatim
00223 *>          LDV1T is INTEGER
00224 *>          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
00225 *>          MAX(1,Q).
00226 *> \endverbatim
00227 *>
00228 *> \param[out] V2T
00229 *> \verbatim
00230 *>          V2T is DOUBLE PRECISION array, dimension (M-Q)
00231 *>          If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) orthogonal
00232 *>          matrix V2**T.
00233 *> \endverbatim
00234 *>
00235 *> \param[in] LDV2T
00236 *> \verbatim
00237 *>          LDV2T is INTEGER
00238 *>          The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
00239 *>          MAX(1,M-Q).
00240 *> \endverbatim
00241 *>
00242 *> \param[out] WORK
00243 *> \verbatim
00244 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00245 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00246 *>          If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
00247 *>          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
00248 *>          define the matrix in intermediate bidiagonal-block form
00249 *>          remaining after nonconvergence. INFO specifies the number
00250 *>          of nonzero PHI's.
00251 *> \endverbatim
00252 *>
00253 *> \param[in] LWORK
00254 *> \verbatim
00255 *>          LWORK is INTEGER
00256 *>          The dimension of the array WORK.
00257 *>
00258 *>          If LWORK = -1, then a workspace query is assumed; the routine
00259 *>          only calculates the optimal size of the WORK array, returns
00260 *>          this value as the first entry of the work array, and no error
00261 *>          message related to LWORK is issued by XERBLA.
00262 *> \endverbatim
00263 *>
00264 *> \param[out] IWORK
00265 *> \verbatim
00266 *>          IWORK is INTEGER array, dimension (M-MIN(P, M-P, Q, M-Q))
00267 *> \endverbatim
00268 *>
00269 *> \param[out] INFO
00270 *> \verbatim
00271 *>          INFO is INTEGER
00272 *>          = 0:  successful exit.
00273 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00274 *>          > 0:  DBBCSD did not converge. See the description of WORK
00275 *>                above for details.
00276 *> \endverbatim
00277 *
00278 *> \par References:
00279 *  ================
00280 *>
00281 *>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
00282 *>      Algorithms, 50(1):33-65, 2009.
00283 *
00284 *  Authors:
00285 *  ========
00286 *
00287 *> \author Univ. of Tennessee 
00288 *> \author Univ. of California Berkeley 
00289 *> \author Univ. of Colorado Denver 
00290 *> \author NAG Ltd. 
00291 *
00292 *> \date November 2011
00293 *
00294 *> \ingroup doubleOTHERcomputational
00295 *
00296 *  =====================================================================
00297       RECURSIVE SUBROUTINE DORCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
00298      $                             SIGNS, M, P, Q, X11, LDX11, X12,
00299      $                             LDX12, X21, LDX21, X22, LDX22, THETA,
00300      $                             U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
00301      $                             LDV2T, WORK, LWORK, IWORK, INFO )
00302 *
00303 *  -- LAPACK computational routine (version 3.4.0) --
00304 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00305 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00306 *     November 2011
00307 *
00308 *     .. Scalar Arguments ..
00309       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
00310       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
00311      $                   LDX21, LDX22, LWORK, M, P, Q
00312 *     ..
00313 *     .. Array Arguments ..
00314       INTEGER            IWORK( * )
00315       DOUBLE PRECISION   THETA( * )
00316       DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
00317      $                   V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
00318      $                   X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
00319      $                   * )
00320 *     ..
00321 *
00322 *  ===================================================================
00323 *
00324 *     .. Parameters ..
00325       DOUBLE PRECISION   ONE, ZERO
00326       PARAMETER          ( ONE = 1.0D0,
00327      $                     ZERO = 0.0D0 )
00328 *     ..
00329 *     .. Local Scalars ..
00330       CHARACTER          TRANST, SIGNST
00331       INTEGER            CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
00332      $                   IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
00333      $                   IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
00334      $                   ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN,
00335      $                   LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN,
00336      $                   LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN,
00337      $                   LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN,
00338      $                   LORGQRWORKOPT, LWORKMIN, LWORKOPT
00339       LOGICAL            COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
00340      $                   WANTV1T, WANTV2T
00341 *     ..
00342 *     .. External Subroutines ..
00343       EXTERNAL           DBBCSD, DLACPY, DLAPMR, DLAPMT, DLASCL, DLASET,
00344      $                   DORBDB, DORGLQ, DORGQR, XERBLA
00345 *     ..
00346 *     .. External Functions ..
00347       LOGICAL            LSAME
00348       EXTERNAL           LSAME
00349 *     ..
00350 *     .. Intrinsic Functions
00351       INTRINSIC          INT, MAX, MIN
00352 *     ..
00353 *     .. Executable Statements ..
00354 *
00355 *     Test input arguments
00356 *
00357       INFO = 0
00358       WANTU1 = LSAME( JOBU1, 'Y' )
00359       WANTU2 = LSAME( JOBU2, 'Y' )
00360       WANTV1T = LSAME( JOBV1T, 'Y' )
00361       WANTV2T = LSAME( JOBV2T, 'Y' )
00362       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
00363       DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' )
00364       LQUERY = LWORK .EQ. -1
00365       IF( M .LT. 0 ) THEN
00366          INFO = -7
00367       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
00368          INFO = -8
00369       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
00370          INFO = -9
00371       ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR.
00372      $         ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN
00373          INFO = -11
00374       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
00375          INFO = -20
00376       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
00377          INFO = -22
00378       ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
00379          INFO = -24
00380       ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
00381          INFO = -26
00382       END IF
00383 *
00384 *     Work with transpose if convenient
00385 *
00386       IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN
00387          IF( COLMAJOR ) THEN
00388             TRANST = 'T'
00389          ELSE
00390             TRANST = 'N'
00391          END IF
00392          IF( DEFAULTSIGNS ) THEN
00393             SIGNST = 'O'
00394          ELSE
00395             SIGNST = 'D'
00396          END IF
00397          CALL DORCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M,
00398      $                Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22,
00399      $                LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1,
00400      $                U2, LDU2, WORK, LWORK, IWORK, INFO )
00401          RETURN
00402       END IF
00403 *
00404 *     Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
00405 *     convenient
00406 *
00407       IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN
00408          IF( DEFAULTSIGNS ) THEN
00409             SIGNST = 'O'
00410          ELSE
00411             SIGNST = 'D'
00412          END IF
00413          CALL DORCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M,
00414      $                M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11,
00415      $                LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T,
00416      $                LDV1T, WORK, LWORK, IWORK, INFO )
00417          RETURN
00418       END IF
00419 *
00420 *     Compute workspace
00421 *
00422       IF( INFO .EQ. 0 ) THEN
00423 *
00424          IPHI = 2
00425          ITAUP1 = IPHI + MAX( 1, Q - 1 )
00426          ITAUP2 = ITAUP1 + MAX( 1, P )
00427          ITAUQ1 = ITAUP2 + MAX( 1, M - P )
00428          ITAUQ2 = ITAUQ1 + MAX( 1, Q )
00429          IORGQR = ITAUQ2 + MAX( 1, M - Q )
00430          CALL DORGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
00431      $                CHILDINFO )
00432          LORGQRWORKOPT = INT( WORK(1) )
00433          LORGQRWORKMIN = MAX( 1, M - Q )
00434          IORGLQ = ITAUQ2 + MAX( 1, M - Q )
00435          CALL DORGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
00436      $                CHILDINFO )
00437          LORGLQWORKOPT = INT( WORK(1) )
00438          LORGLQWORKMIN = MAX( 1, M - Q )
00439          IORBDB = ITAUQ2 + MAX( 1, M - Q )
00440          CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
00441      $                X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK,
00442      $                -1, CHILDINFO )
00443          LORBDBWORKOPT = INT( WORK(1) )
00444          LORBDBWORKMIN = LORBDBWORKOPT
00445          IB11D = ITAUQ2 + MAX( 1, M - Q )
00446          IB11E = IB11D + MAX( 1, Q )
00447          IB12D = IB11E + MAX( 1, Q - 1 )
00448          IB12E = IB12D + MAX( 1, Q )
00449          IB21D = IB12E + MAX( 1, Q - 1 )
00450          IB21E = IB21D + MAX( 1, Q )
00451          IB22D = IB21E + MAX( 1, Q - 1 )
00452          IB22E = IB22D + MAX( 1, Q )
00453          IBBCSD = IB22E + MAX( 1, Q - 1 )
00454          CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0,
00455      $                0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0,
00456      $                0, 0, 0, 0, 0, 0, 0, WORK, -1, CHILDINFO )
00457          LBBCSDWORKOPT = INT( WORK(1) )
00458          LBBCSDWORKMIN = LBBCSDWORKOPT
00459          LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT,
00460      $              IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKOPT ) - 1
00461          LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN,
00462      $              IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKMIN ) - 1
00463          WORK(1) = MAX(LWORKOPT,LWORKMIN)
00464 *
00465          IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
00466             INFO = -22
00467          ELSE
00468             LORGQRWORK = LWORK - IORGQR + 1
00469             LORGLQWORK = LWORK - IORGLQ + 1
00470             LORBDBWORK = LWORK - IORBDB + 1
00471             LBBCSDWORK = LWORK - IBBCSD + 1
00472          END IF
00473       END IF
00474 *
00475 *     Abort if any illegal arguments
00476 *
00477       IF( INFO .NE. 0 ) THEN
00478          CALL XERBLA( 'DORCSD', -INFO )
00479          RETURN
00480       ELSE IF( LQUERY ) THEN
00481          RETURN
00482       END IF
00483 *
00484 *     Transform to bidiagonal block form
00485 *
00486       CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21,
00487      $             LDX21, X22, LDX22, THETA, WORK(IPHI), WORK(ITAUP1),
00488      $             WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2),
00489      $             WORK(IORBDB), LORBDBWORK, CHILDINFO )
00490 *
00491 *     Accumulate Householder reflectors
00492 *
00493       IF( COLMAJOR ) THEN
00494          IF( WANTU1 .AND. P .GT. 0 ) THEN
00495             CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
00496             CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
00497      $                   LORGQRWORK, INFO)
00498          END IF
00499          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
00500             CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
00501             CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
00502      $                   WORK(IORGQR), LORGQRWORK, INFO )
00503          END IF
00504          IF( WANTV1T .AND. Q .GT. 0 ) THEN
00505             CALL DLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2),
00506      $                   LDV1T )
00507             V1T(1, 1) = ONE
00508             DO J = 2, Q
00509                V1T(1,J) = ZERO
00510                V1T(J,1) = ZERO
00511             END DO
00512             CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
00513      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00514          END IF
00515          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
00516             CALL DLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T )
00517             CALL DLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
00518      $                   V2T(P+1,P+1), LDV2T )
00519             CALL DORGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
00520      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00521          END IF
00522       ELSE
00523          IF( WANTU1 .AND. P .GT. 0 ) THEN
00524             CALL DLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 )
00525             CALL DORGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ),
00526      $                   LORGLQWORK, INFO)
00527          END IF
00528          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
00529             CALL DLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 )
00530             CALL DORGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
00531      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00532          END IF
00533          IF( WANTV1T .AND. Q .GT. 0 ) THEN
00534             CALL DLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2),
00535      $                   LDV1T )
00536             V1T(1, 1) = ONE
00537             DO J = 2, Q
00538                V1T(1,J) = ZERO
00539                V1T(J,1) = ZERO
00540             END DO
00541             CALL DORGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
00542      $                   WORK(IORGQR), LORGQRWORK, INFO )
00543          END IF
00544          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
00545             CALL DLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T )
00546             CALL DLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22,
00547      $                   V2T(P+1,P+1), LDV2T )
00548             CALL DORGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
00549      $                   WORK(IORGQR), LORGQRWORK, INFO )
00550          END IF
00551       END IF
00552 *
00553 *     Compute the CSD of the matrix in bidiagonal-block form
00554 *
00555       CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA,
00556      $             WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
00557      $             LDV2T, WORK(IB11D), WORK(IB11E), WORK(IB12D),
00558      $             WORK(IB12E), WORK(IB21D), WORK(IB21E), WORK(IB22D),
00559      $             WORK(IB22E), WORK(IBBCSD), LBBCSDWORK, INFO )
00560 *
00561 *     Permute rows and columns to place identity submatrices in top-
00562 *     left corner of (1,1)-block and/or bottom-right corner of (1,2)-
00563 *     block and/or bottom-right corner of (2,1)-block and/or top-left
00564 *     corner of (2,2)-block 
00565 *
00566       IF( Q .GT. 0 .AND. WANTU2 ) THEN
00567          DO I = 1, Q
00568             IWORK(I) = M - P - Q + I
00569          END DO
00570          DO I = Q + 1, M - P
00571             IWORK(I) = I - Q
00572          END DO
00573          IF( COLMAJOR ) THEN
00574             CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
00575          ELSE
00576             CALL DLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK )
00577          END IF
00578       END IF
00579       IF( M .GT. 0 .AND. WANTV2T ) THEN
00580          DO I = 1, P
00581             IWORK(I) = M - P - Q + I
00582          END DO
00583          DO I = P + 1, M - Q
00584             IWORK(I) = I - P
00585          END DO
00586          IF( .NOT. COLMAJOR ) THEN
00587             CALL DLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
00588          ELSE
00589             CALL DLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
00590          END IF
00591       END IF
00592 *
00593       RETURN
00594 *
00595 *     End DORCSD
00596 *
00597       END
00598 
 All Files Functions