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