LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sorcsd.f
Go to the documentation of this file.
00001 *> \brief \b SORCSD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SORCSD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorcsd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorcsd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorcsd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       RECURSIVE SUBROUTINE SORCSD( 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 *       REAL               THETA( * )
00035 *       REAL               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 *> SORCSD 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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:  SBBCSD 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 realOTHERcomputational
00295 *
00296 *  =====================================================================
00297       RECURSIVE SUBROUTINE SORCSD( 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       REAL               THETA( * )
00316       REAL               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       REAL               ONE, ZERO
00326       PARAMETER          ( ONE = 1.0E+0,
00327      $                     ZERO = 0.0E+0 )
00328 *     ..
00329 *     .. Local Arrays ..
00330       REAL               DUMMY(1)
00331 *     ..
00332 *     .. Local Scalars ..
00333       CHARACTER          TRANST, SIGNST
00334       INTEGER            CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
00335      $                   IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
00336      $                   IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
00337      $                   ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN,
00338      $                   LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN,
00339      $                   LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN,
00340      $                   LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN,
00341      $                   LORGQRWORKOPT, LWORKMIN, LWORKOPT
00342       LOGICAL            COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
00343      $                   WANTV1T, WANTV2T
00344 *     ..
00345 *     .. External Subroutines ..
00346       EXTERNAL           SBBCSD, SLACPY, SLAPMR, SLAPMT, SLASCL, SLASET,
00347      $                   SORBDB, SORGLQ, SORGQR, XERBLA
00348 *     ..
00349 *     .. External Functions ..
00350       LOGICAL            LSAME
00351       EXTERNAL           LSAME
00352 *     ..
00353 *     .. Intrinsic Functions
00354       INTRINSIC          INT, MAX, MIN
00355 *     ..
00356 *     .. Executable Statements ..
00357 *
00358 *     Test input arguments
00359 *
00360       INFO = 0
00361       WANTU1 = LSAME( JOBU1, 'Y' )
00362       WANTU2 = LSAME( JOBU2, 'Y' )
00363       WANTV1T = LSAME( JOBV1T, 'Y' )
00364       WANTV2T = LSAME( JOBV2T, 'Y' )
00365       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
00366       DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' )
00367       LQUERY = LWORK .EQ. -1
00368       IF( M .LT. 0 ) THEN
00369          INFO = -7
00370       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
00371          INFO = -8
00372       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
00373          INFO = -9
00374       ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR.
00375      $         ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN
00376          INFO = -11
00377       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
00378          INFO = -20
00379       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
00380          INFO = -22
00381       ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
00382          INFO = -24
00383       ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
00384          INFO = -26
00385       END IF
00386 *
00387 *     Work with transpose if convenient
00388 *
00389       IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN
00390          IF( COLMAJOR ) THEN
00391             TRANST = 'T'
00392          ELSE
00393             TRANST = 'N'
00394          END IF
00395          IF( DEFAULTSIGNS ) THEN
00396             SIGNST = 'O'
00397          ELSE
00398             SIGNST = 'D'
00399          END IF
00400          CALL SORCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M,
00401      $                Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22,
00402      $                LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1,
00403      $                U2, LDU2, WORK, LWORK, IWORK, INFO )
00404          RETURN
00405       END IF
00406 *
00407 *     Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
00408 *     convenient
00409 *
00410       IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN
00411          IF( DEFAULTSIGNS ) THEN
00412             SIGNST = 'O'
00413          ELSE
00414             SIGNST = 'D'
00415          END IF
00416          CALL SORCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M,
00417      $                M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11,
00418      $                LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T,
00419      $                LDV1T, WORK, LWORK, IWORK, INFO )
00420          RETURN
00421       END IF
00422 *
00423 *     Compute workspace
00424 *
00425       IF( INFO .EQ. 0 ) THEN
00426 *
00427          IPHI = 2
00428          ITAUP1 = IPHI + MAX( 1, Q - 1 )
00429          ITAUP2 = ITAUP1 + MAX( 1, P )
00430          ITAUQ1 = ITAUP2 + MAX( 1, M - P )
00431          ITAUQ2 = ITAUQ1 + MAX( 1, Q )
00432          IORGQR = ITAUQ2 + MAX( 1, M - Q )
00433          CALL SORGQR( M-Q, M-Q, M-Q, DUMMY, MAX(1,M-Q), DUMMY, WORK, -1,
00434      $                CHILDINFO )
00435          LORGQRWORKOPT = INT( WORK(1) )
00436          LORGQRWORKMIN = MAX( 1, M - Q )
00437          IORGLQ = ITAUQ2 + MAX( 1, M - Q )
00438          CALL SORGLQ( M-Q, M-Q, M-Q, DUMMY, MAX(1,M-Q), DUMMY, WORK, -1,
00439      $                CHILDINFO )
00440          LORGLQWORKOPT = INT( WORK(1) )
00441          LORGLQWORKMIN = MAX( 1, M - Q )
00442          IORBDB = ITAUQ2 + MAX( 1, M - Q )
00443          CALL SORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
00444      $        X21, LDX21, X22, LDX22, DUMMY, DUMMY, DUMMY, DUMMY, DUMMY,
00445      $        DUMMY,WORK,-1,CHILDINFO )
00446          LORBDBWORKOPT = INT( WORK(1) )
00447          LORBDBWORKMIN = LORBDBWORKOPT
00448          IB11D = ITAUQ2 + MAX( 1, M - Q )
00449          IB11E = IB11D + MAX( 1, Q )
00450          IB12D = IB11E + MAX( 1, Q - 1 )
00451          IB12E = IB12D + MAX( 1, Q )
00452          IB21D = IB12E + MAX( 1, Q - 1 )
00453          IB21E = IB21D + MAX( 1, Q )
00454          IB22D = IB21E + MAX( 1, Q - 1 )
00455          IB22E = IB22D + MAX( 1, Q )
00456          IBBCSD = IB22E + MAX( 1, Q - 1 )
00457          CALL SBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
00458      $                DUMMY, DUMMY, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
00459      $                LDV2T, DUMMY, DUMMY, DUMMY, DUMMY, DUMMY, DUMMY,
00460      $                DUMMY, DUMMY, WORK, -1, CHILDINFO )
00461          LBBCSDWORKOPT = INT( WORK(1) )
00462          LBBCSDWORKMIN = LBBCSDWORKOPT
00463          LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT,
00464      $              IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKOPT ) - 1
00465          LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN,
00466      $              IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKMIN ) - 1
00467          WORK(1) = MAX(LWORKOPT,LWORKMIN)
00468 *
00469          IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
00470             INFO = -22
00471          ELSE
00472             LORGQRWORK = LWORK - IORGQR + 1
00473             LORGLQWORK = LWORK - IORGLQ + 1
00474             LORBDBWORK = LWORK - IORBDB + 1
00475             LBBCSDWORK = LWORK - IBBCSD + 1
00476          END IF
00477       END IF
00478 *
00479 *     Abort if any illegal arguments
00480 *
00481       IF( INFO .NE. 0 ) THEN
00482          CALL XERBLA( 'SORCSD', -INFO )
00483          RETURN
00484       ELSE IF( LQUERY ) THEN
00485          RETURN
00486       END IF
00487 *
00488 *     Transform to bidiagonal block form
00489 *
00490       CALL SORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21,
00491      $             LDX21, X22, LDX22, THETA, WORK(IPHI), WORK(ITAUP1),
00492      $             WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2),
00493      $             WORK(IORBDB), LORBDBWORK, CHILDINFO )
00494 *
00495 *     Accumulate Householder reflectors
00496 *
00497       IF( COLMAJOR ) THEN
00498          IF( WANTU1 .AND. P .GT. 0 ) THEN
00499             CALL SLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
00500             CALL SORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
00501      $                   LORGQRWORK, INFO)
00502          END IF
00503          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
00504             CALL SLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
00505             CALL SORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
00506      $                   WORK(IORGQR), LORGQRWORK, INFO )
00507          END IF
00508          IF( WANTV1T .AND. Q .GT. 0 ) THEN
00509             CALL SLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2),
00510      $                   LDV1T )
00511             V1T(1, 1) = ONE
00512             DO J = 2, Q
00513                V1T(1,J) = ZERO
00514                V1T(J,1) = ZERO
00515             END DO
00516             CALL SORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
00517      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00518          END IF
00519          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
00520             CALL SLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T )
00521             CALL SLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
00522      $                   V2T(P+1,P+1), LDV2T )
00523             CALL SORGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
00524      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00525          END IF
00526       ELSE
00527          IF( WANTU1 .AND. P .GT. 0 ) THEN
00528             CALL SLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 )
00529             CALL SORGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ),
00530      $                   LORGLQWORK, INFO)
00531          END IF
00532          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
00533             CALL SLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 )
00534             CALL SORGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
00535      $                   WORK(IORGLQ), LORGLQWORK, INFO )
00536          END IF
00537          IF( WANTV1T .AND. Q .GT. 0 ) THEN
00538             CALL SLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2),
00539      $                   LDV1T )
00540             V1T(1, 1) = ONE
00541             DO J = 2, Q
00542                V1T(1,J) = ZERO
00543                V1T(J,1) = ZERO
00544             END DO
00545             CALL SORGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
00546      $                   WORK(IORGQR), LORGQRWORK, INFO )
00547          END IF
00548          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
00549             CALL SLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T )
00550             CALL SLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22,
00551      $                   V2T(P+1,P+1), LDV2T )
00552             CALL SORGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
00553      $                   WORK(IORGQR), LORGQRWORK, INFO )
00554          END IF
00555       END IF
00556 *
00557 *     Compute the CSD of the matrix in bidiagonal-block form
00558 *
00559       CALL SBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA,
00560      $             WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
00561      $             LDV2T, WORK(IB11D), WORK(IB11E), WORK(IB12D),
00562      $             WORK(IB12E), WORK(IB21D), WORK(IB21E), WORK(IB22D),
00563      $             WORK(IB22E), WORK(IBBCSD), LBBCSDWORK, INFO )
00564 *
00565 *     Permute rows and columns to place identity submatrices in top-
00566 *     left corner of (1,1)-block and/or bottom-right corner of (1,2)-
00567 *     block and/or bottom-right corner of (2,1)-block and/or top-left
00568 *     corner of (2,2)-block 
00569 *
00570       IF( Q .GT. 0 .AND. WANTU2 ) THEN
00571          DO I = 1, Q
00572             IWORK(I) = M - P - Q + I
00573          END DO
00574          DO I = Q + 1, M - P
00575             IWORK(I) = I - Q
00576          END DO
00577          IF( COLMAJOR ) THEN
00578             CALL SLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
00579          ELSE
00580             CALL SLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK )
00581          END IF
00582       END IF
00583       IF( M .GT. 0 .AND. WANTV2T ) THEN
00584          DO I = 1, P
00585             IWORK(I) = M - P - Q + I
00586          END DO
00587          DO I = P + 1, M - Q
00588             IWORK(I) = I - P
00589          END DO
00590          IF( .NOT. COLMAJOR ) THEN
00591             CALL SLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
00592          ELSE
00593             CALL SLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
00594          END IF
00595       END IF
00596 *
00597       RETURN
00598 *
00599 *     End SORCSD
00600 *
00601       END
00602 
 All Files Functions