Annotation of rpl/lapack/lapack/zuncsd.f, revision 1.7

1.4       bertrand    1: *> \brief \b ZUNCSD
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
                      5: * Online html documentation available at 
                      6: *            http://www.netlib.org/lapack/explore-html/ 
                      7: *
                      8: *> \htmlonly
                      9: *> Download ZUNCSD + dependencies 
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zuncsd.f"> 
                     11: *> [TGZ]</a> 
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zuncsd.f"> 
                     13: *> [ZIP]</a> 
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zuncsd.f"> 
                     15: *> [TXT]</a>
                     16: *> \endhtmlonly 
                     17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       RECURSIVE SUBROUTINE ZUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
                     22: *                                    SIGNS, M, P, Q, X11, LDX11, X12,
                     23: *                                    LDX12, X21, LDX21, X22, LDX22, THETA,
                     24: *                                    U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
                     25: *                                    LDV2T, WORK, LWORK, RWORK, LRWORK,
                     26: *                                    IWORK, INFO )
                     27: * 
                     28: *       .. Scalar Arguments ..
                     29: *       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
                     30: *       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
                     31: *      $                   LDX21, LDX22, LRWORK, LWORK, M, P, Q
                     32: *       ..
                     33: *       .. Array Arguments ..
                     34: *       INTEGER            IWORK( * )
                     35: *       DOUBLE PRECISION   THETA( * )
                     36: *       DOUBLE PRECISION   RWORK( * )
                     37: *       COMPLEX*16         U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
                     38: *      $                   V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
                     39: *      $                   X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
                     40: *      $                   * )
                     41: *       ..
                     42: *  
                     43: *
                     44: *> \par Purpose:
                     45: *  =============
                     46: *>
                     47: *> \verbatim
                     48: *>
                     49: *> ZUNCSD computes the CS decomposition of an M-by-M partitioned
                     50: *> unitary matrix X:
                     51: *>
                     52: *>                                 [  I  0  0 |  0  0  0 ]
                     53: *>                                 [  0  C  0 |  0 -S  0 ]
                     54: *>     [ X11 | X12 ]   [ U1 |    ] [  0  0  0 |  0  0 -I ] [ V1 |    ]**H
                     55: *> X = [-----------] = [---------] [---------------------] [---------]   .
                     56: *>     [ X21 | X22 ]   [    | U2 ] [  0  0  0 |  I  0  0 ] [    | V2 ]
                     57: *>                                 [  0  S  0 |  0  C  0 ]
                     58: *>                                 [  0  0  I |  0  0  0 ]
                     59: *>
                     60: *> X11 is P-by-Q. The unitary matrices U1, U2, V1, and V2 are P-by-P,
                     61: *> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
                     62: *> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
                     63: *> which R = MIN(P,M-P,Q,M-Q).
                     64: *> \endverbatim
                     65: *
                     66: *  Arguments:
                     67: *  ==========
                     68: *
                     69: *> \param[in] JOBU1
                     70: *> \verbatim
                     71: *>          JOBU1 is CHARACTER
                     72: *>          = 'Y':      U1 is computed;
                     73: *>          otherwise:  U1 is not computed.
                     74: *> \endverbatim
                     75: *>
                     76: *> \param[in] JOBU2
                     77: *> \verbatim
                     78: *>          JOBU2 is CHARACTER
                     79: *>          = 'Y':      U2 is computed;
                     80: *>          otherwise:  U2 is not computed.
                     81: *> \endverbatim
                     82: *>
                     83: *> \param[in] JOBV1T
                     84: *> \verbatim
                     85: *>          JOBV1T is CHARACTER
                     86: *>          = 'Y':      V1T is computed;
                     87: *>          otherwise:  V1T is not computed.
                     88: *> \endverbatim
                     89: *>
                     90: *> \param[in] JOBV2T
                     91: *> \verbatim
                     92: *>          JOBV2T is CHARACTER
                     93: *>          = 'Y':      V2T is computed;
                     94: *>          otherwise:  V2T is not computed.
                     95: *> \endverbatim
                     96: *>
                     97: *> \param[in] TRANS
                     98: *> \verbatim
                     99: *>          TRANS is CHARACTER
                    100: *>          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
                    101: *>                      order;
                    102: *>          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
                    103: *>                      major order.
                    104: *> \endverbatim
                    105: *>
                    106: *> \param[in] SIGNS
                    107: *> \verbatim
                    108: *>          SIGNS is CHARACTER
                    109: *>          = 'O':      The lower-left block is made nonpositive (the
                    110: *>                      "other" convention);
                    111: *>          otherwise:  The upper-right block is made nonpositive (the
                    112: *>                      "default" convention).
                    113: *> \endverbatim
                    114: *>
                    115: *> \param[in] M
                    116: *> \verbatim
                    117: *>          M is INTEGER
                    118: *>          The number of rows and columns in X.
                    119: *> \endverbatim
                    120: *>
                    121: *> \param[in] P
                    122: *> \verbatim
                    123: *>          P is INTEGER
                    124: *>          The number of rows in X11 and X12. 0 <= P <= M.
                    125: *> \endverbatim
                    126: *>
                    127: *> \param[in] Q
                    128: *> \verbatim
                    129: *>          Q is INTEGER
                    130: *>          The number of columns in X11 and X21. 0 <= Q <= M.
                    131: *> \endverbatim
                    132: *>
                    133: *> \param[in,out] X11
                    134: *> \verbatim
                    135: *>          X11 is COMPLEX*16 array, dimension (LDX11,Q)
                    136: *>          On entry, part of the unitary matrix whose CSD is desired.
                    137: *> \endverbatim
                    138: *>
                    139: *> \param[in] LDX11
                    140: *> \verbatim
                    141: *>          LDX11 is INTEGER
                    142: *>          The leading dimension of X11. LDX11 >= MAX(1,P).
                    143: *> \endverbatim
                    144: *>
                    145: *> \param[in,out] X12
                    146: *> \verbatim
                    147: *>          X12 is COMPLEX*16 array, dimension (LDX12,M-Q)
                    148: *>          On entry, part of the unitary matrix whose CSD is desired.
                    149: *> \endverbatim
                    150: *>
                    151: *> \param[in] LDX12
                    152: *> \verbatim
                    153: *>          LDX12 is INTEGER
                    154: *>          The leading dimension of X12. LDX12 >= MAX(1,P).
                    155: *> \endverbatim
                    156: *>
                    157: *> \param[in,out] X21
                    158: *> \verbatim
                    159: *>          X21 is COMPLEX*16 array, dimension (LDX21,Q)
                    160: *>          On entry, part of the unitary matrix whose CSD is desired.
                    161: *> \endverbatim
                    162: *>
                    163: *> \param[in] LDX21
                    164: *> \verbatim
                    165: *>          LDX21 is INTEGER
                    166: *>          The leading dimension of X11. LDX21 >= MAX(1,M-P).
                    167: *> \endverbatim
                    168: *>
                    169: *> \param[in,out] X22
                    170: *> \verbatim
                    171: *>          X22 is COMPLEX*16 array, dimension (LDX22,M-Q)
                    172: *>          On entry, part of the unitary matrix whose CSD is desired.
                    173: *> \endverbatim
                    174: *>
                    175: *> \param[in] LDX22
                    176: *> \verbatim
                    177: *>          LDX22 is INTEGER
                    178: *>          The leading dimension of X11. LDX22 >= MAX(1,M-P).
                    179: *> \endverbatim
                    180: *>
                    181: *> \param[out] THETA
                    182: *> \verbatim
                    183: *>          THETA is DOUBLE PRECISION array, dimension (R), in which R =
                    184: *>          MIN(P,M-P,Q,M-Q).
                    185: *>          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
                    186: *>          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
                    187: *> \endverbatim
                    188: *>
                    189: *> \param[out] U1
                    190: *> \verbatim
                    191: *>          U1 is COMPLEX*16 array, dimension (P)
                    192: *>          If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1.
                    193: *> \endverbatim
                    194: *>
                    195: *> \param[in] LDU1
                    196: *> \verbatim
                    197: *>          LDU1 is INTEGER
                    198: *>          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
                    199: *>          MAX(1,P).
                    200: *> \endverbatim
                    201: *>
                    202: *> \param[out] U2
                    203: *> \verbatim
                    204: *>          U2 is COMPLEX*16 array, dimension (M-P)
                    205: *>          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
                    206: *>          matrix U2.
                    207: *> \endverbatim
                    208: *>
                    209: *> \param[in] LDU2
                    210: *> \verbatim
                    211: *>          LDU2 is INTEGER
                    212: *>          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
                    213: *>          MAX(1,M-P).
                    214: *> \endverbatim
                    215: *>
                    216: *> \param[out] V1T
                    217: *> \verbatim
                    218: *>          V1T is COMPLEX*16 array, dimension (Q)
                    219: *>          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
                    220: *>          matrix V1**H.
                    221: *> \endverbatim
                    222: *>
                    223: *> \param[in] LDV1T
                    224: *> \verbatim
                    225: *>          LDV1T is INTEGER
                    226: *>          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
                    227: *>          MAX(1,Q).
                    228: *> \endverbatim
                    229: *>
                    230: *> \param[out] V2T
                    231: *> \verbatim
                    232: *>          V2T is COMPLEX*16 array, dimension (M-Q)
                    233: *>          If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) unitary
                    234: *>          matrix V2**H.
                    235: *> \endverbatim
                    236: *>
                    237: *> \param[in] LDV2T
                    238: *> \verbatim
                    239: *>          LDV2T is INTEGER
                    240: *>          The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
                    241: *>          MAX(1,M-Q).
                    242: *> \endverbatim
                    243: *>
                    244: *> \param[out] WORK
                    245: *> \verbatim
                    246: *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                    247: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
                    248: *> \endverbatim
                    249: *>
                    250: *> \param[in] LWORK
                    251: *> \verbatim
                    252: *>          LWORK is INTEGER
                    253: *>          The dimension of the array WORK.
                    254: *>
                    255: *>          If LWORK = -1, then a workspace query is assumed; the routine
                    256: *>          only calculates the optimal size of the WORK array, returns
                    257: *>          this value as the first entry of the work array, and no error
                    258: *>          message related to LWORK is issued by XERBLA.
                    259: *> \endverbatim
                    260: *>
                    261: *> \param[out] RWORK
                    262: *> \verbatim
                    263: *>          RWORK is DOUBLE PRECISION array, dimension MAX(1,LRWORK)
                    264: *>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
                    265: *>          If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
                    266: *>          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
                    267: *>          define the matrix in intermediate bidiagonal-block form
                    268: *>          remaining after nonconvergence. INFO specifies the number
                    269: *>          of nonzero PHI's.
                    270: *> \endverbatim
                    271: *>
                    272: *> \param[in] LRWORK
                    273: *> \verbatim
                    274: *>          LRWORK is INTEGER
                    275: *>          The dimension of the array RWORK.
                    276: *>
                    277: *>          If LRWORK = -1, then a workspace query is assumed; the routine
                    278: *>          only calculates the optimal size of the RWORK array, returns
                    279: *>          this value as the first entry of the work array, and no error
                    280: *>          message related to LRWORK is issued by XERBLA.
                    281: *> \endverbatim
                    282: *>
                    283: *> \param[out] IWORK
                    284: *> \verbatim
                    285: *>          IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
                    286: *> \endverbatim
                    287: *>
                    288: *> \param[out] INFO
                    289: *> \verbatim
                    290: *>          INFO is INTEGER
                    291: *>          = 0:  successful exit.
                    292: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
                    293: *>          > 0:  ZBBCSD did not converge. See the description of RWORK
                    294: *>                above for details.
                    295: *> \endverbatim
                    296: *
                    297: *> \par References:
                    298: *  ================
                    299: *>
                    300: *>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
                    301: *>      Algorithms, 50(1):33-65, 2009.
                    302: *
                    303: *  Authors:
                    304: *  ========
                    305: *
                    306: *> \author Univ. of Tennessee 
                    307: *> \author Univ. of California Berkeley 
                    308: *> \author Univ. of Colorado Denver 
                    309: *> \author NAG Ltd. 
                    310: *
                    311: *> \date November 2011
                    312: *
                    313: *> \ingroup complex16OTHERcomputational
                    314: *
                    315: *  =====================================================================
1.1       bertrand  316:       RECURSIVE SUBROUTINE ZUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
                    317:      $                             SIGNS, M, P, Q, X11, LDX11, X12,
                    318:      $                             LDX12, X21, LDX21, X22, LDX22, THETA,
                    319:      $                             U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
                    320:      $                             LDV2T, WORK, LWORK, RWORK, LRWORK,
                    321:      $                             IWORK, INFO )
                    322: *
1.4       bertrand  323: *  -- LAPACK computational routine (version 3.4.0) --
1.1       bertrand  324: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
1.4       bertrand  325: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    326: *     November 2011
1.3       bertrand  327: *
1.1       bertrand  328: *     .. Scalar Arguments ..
                    329:       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
                    330:       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
                    331:      $                   LDX21, LDX22, LRWORK, LWORK, M, P, Q
                    332: *     ..
                    333: *     .. Array Arguments ..
                    334:       INTEGER            IWORK( * )
                    335:       DOUBLE PRECISION   THETA( * )
                    336:       DOUBLE PRECISION   RWORK( * )
                    337:       COMPLEX*16         U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
                    338:      $                   V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
                    339:      $                   X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
                    340:      $                   * )
                    341: *     ..
                    342: *
                    343: *  ===================================================================
                    344: *
                    345: *     .. Parameters ..
1.6       bertrand  346:       COMPLEX*16         ONE, ZERO
                    347:       PARAMETER          ( ONE = (1.0D0,0.0D0),
1.1       bertrand  348:      $                     ZERO = (0.0D0,0.0D0) )
                    349: *     ..
                    350: *     .. Local Scalars ..
                    351:       CHARACTER          TRANST, SIGNST
                    352:       INTEGER            CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
                    353:      $                   IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
                    354:      $                   IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
                    355:      $                   ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN,
                    356:      $                   LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN,
                    357:      $                   LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN,
                    358:      $                   LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN,
                    359:      $                   LORGQRWORKOPT, LWORKMIN, LWORKOPT
                    360:       LOGICAL            COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
                    361:      $                   WANTV1T, WANTV2T
                    362:       INTEGER            LRWORKMIN, LRWORKOPT
                    363:       LOGICAL            LRQUERY
                    364: *     ..
                    365: *     .. External Subroutines ..
                    366:       EXTERNAL           XERBLA, ZBBCSD, ZLACPY, ZLAPMR, ZLAPMT, ZLASCL,
                    367:      $                   ZLASET, ZUNBDB, ZUNGLQ, ZUNGQR
                    368: *     ..
                    369: *     .. External Functions ..
                    370:       LOGICAL            LSAME
                    371:       EXTERNAL           LSAME
                    372: *     ..
                    373: *     .. Intrinsic Functions
                    374:       INTRINSIC          COS, INT, MAX, MIN, SIN
                    375: *     ..
                    376: *     .. Executable Statements ..
                    377: *
                    378: *     Test input arguments
                    379: *
                    380:       INFO = 0
                    381:       WANTU1 = LSAME( JOBU1, 'Y' )
                    382:       WANTU2 = LSAME( JOBU2, 'Y' )
                    383:       WANTV1T = LSAME( JOBV1T, 'Y' )
                    384:       WANTV2T = LSAME( JOBV2T, 'Y' )
                    385:       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
                    386:       DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' )
                    387:       LQUERY = LWORK .EQ. -1
                    388:       LRQUERY = LRWORK .EQ. -1
                    389:       IF( M .LT. 0 ) THEN
                    390:          INFO = -7
                    391:       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
                    392:          INFO = -8
                    393:       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
                    394:          INFO = -9
                    395:       ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR.
                    396:      $         ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN
                    397:          INFO = -11
                    398:       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
1.4       bertrand  399:          INFO = -20
1.1       bertrand  400:       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
1.4       bertrand  401:          INFO = -22
1.1       bertrand  402:       ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
1.4       bertrand  403:          INFO = -24
1.1       bertrand  404:       ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
1.4       bertrand  405:          INFO = -26
1.1       bertrand  406:       END IF
                    407: *
                    408: *     Work with transpose if convenient
                    409: *
                    410:       IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN
                    411:          IF( COLMAJOR ) THEN
                    412:             TRANST = 'T'
                    413:          ELSE
                    414:             TRANST = 'N'
                    415:          END IF
                    416:          IF( DEFAULTSIGNS ) THEN
                    417:             SIGNST = 'O'
                    418:          ELSE
                    419:             SIGNST = 'D'
                    420:          END IF
                    421:          CALL ZUNCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M,
                    422:      $                Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22,
                    423:      $                LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1,
                    424:      $                U2, LDU2, WORK, LWORK, RWORK, LRWORK, IWORK,
                    425:      $                INFO )
                    426:          RETURN
                    427:       END IF
                    428: *
                    429: *     Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
                    430: *     convenient
                    431: *
                    432:       IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN
                    433:          IF( DEFAULTSIGNS ) THEN
                    434:             SIGNST = 'O'
                    435:          ELSE
                    436:             SIGNST = 'D'
                    437:          END IF
                    438:          CALL ZUNCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M,
                    439:      $                M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11,
                    440:      $                LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T,
                    441:      $                LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO )
                    442:          RETURN
                    443:       END IF
                    444: *
                    445: *     Compute workspace
                    446: *
                    447:       IF( INFO .EQ. 0 ) THEN
                    448: *
                    449: *        Real workspace
                    450: *
                    451:          IPHI = 2
                    452:          IB11D = IPHI + MAX( 1, Q - 1 )
                    453:          IB11E = IB11D + MAX( 1, Q )
                    454:          IB12D = IB11E + MAX( 1, Q - 1 )
                    455:          IB12E = IB12D + MAX( 1, Q )
                    456:          IB21D = IB12E + MAX( 1, Q - 1 )
                    457:          IB21E = IB21D + MAX( 1, Q )
                    458:          IB22D = IB21E + MAX( 1, Q - 1 )
                    459:          IB22E = IB22D + MAX( 1, Q )
                    460:          IBBCSD = IB22E + MAX( 1, Q - 1 )
                    461:          CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0,
                    462:      $                0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0,
                    463:      $                0, 0, 0, 0, 0, 0, 0, RWORK, -1, CHILDINFO )
                    464:          LBBCSDWORKOPT = INT( RWORK(1) )
                    465:          LBBCSDWORKMIN = LBBCSDWORKOPT
                    466:          LRWORKOPT = IBBCSD + LBBCSDWORKOPT - 1
                    467:          LRWORKMIN = IBBCSD + LBBCSDWORKMIN - 1
                    468:          RWORK(1) = LRWORKOPT
                    469: *
                    470: *        Complex workspace
                    471: *
                    472:          ITAUP1 = 2
                    473:          ITAUP2 = ITAUP1 + MAX( 1, P )
                    474:          ITAUQ1 = ITAUP2 + MAX( 1, M - P )
                    475:          ITAUQ2 = ITAUQ1 + MAX( 1, Q )
                    476:          IORGQR = ITAUQ2 + MAX( 1, M - Q )
                    477:          CALL ZUNGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
                    478:      $                CHILDINFO )
                    479:          LORGQRWORKOPT = INT( WORK(1) )
                    480:          LORGQRWORKMIN = MAX( 1, M - Q )
                    481:          IORGLQ = ITAUQ2 + MAX( 1, M - Q )
                    482:          CALL ZUNGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
                    483:      $                CHILDINFO )
                    484:          LORGLQWORKOPT = INT( WORK(1) )
                    485:          LORGLQWORKMIN = MAX( 1, M - Q )
                    486:          IORBDB = ITAUQ2 + MAX( 1, M - Q )
                    487:          CALL ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
                    488:      $                X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK,
                    489:      $                -1, CHILDINFO )
                    490:          LORBDBWORKOPT = INT( WORK(1) )
                    491:          LORBDBWORKMIN = LORBDBWORKOPT
                    492:          LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT,
                    493:      $              IORBDB + LORBDBWORKOPT ) - 1
                    494:          LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN,
                    495:      $              IORBDB + LORBDBWORKMIN ) - 1
1.3       bertrand  496:          WORK(1) = MAX(LWORKOPT,LWORKMIN)
1.1       bertrand  497: *
                    498:          IF( LWORK .LT. LWORKMIN
                    499:      $       .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN
                    500:             INFO = -22
                    501:          ELSE IF( LRWORK .LT. LRWORKMIN
                    502:      $            .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN
                    503:             INFO = -24
                    504:          ELSE
                    505:             LORGQRWORK = LWORK - IORGQR + 1
                    506:             LORGLQWORK = LWORK - IORGLQ + 1
                    507:             LORBDBWORK = LWORK - IORBDB + 1
                    508:             LBBCSDWORK = LRWORK - IBBCSD + 1
                    509:          END IF
                    510:       END IF
                    511: *
                    512: *     Abort if any illegal arguments
                    513: *
                    514:       IF( INFO .NE. 0 ) THEN
                    515:          CALL XERBLA( 'ZUNCSD', -INFO )
                    516:          RETURN
                    517:       ELSE IF( LQUERY .OR. LRQUERY ) THEN
                    518:          RETURN
                    519:       END IF
                    520: *
                    521: *     Transform to bidiagonal block form
                    522: *
                    523:       CALL ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21,
                    524:      $             LDX21, X22, LDX22, THETA, RWORK(IPHI), WORK(ITAUP1),
                    525:      $             WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2),
                    526:      $             WORK(IORBDB), LORBDBWORK, CHILDINFO )
                    527: *
                    528: *     Accumulate Householder reflectors
                    529: *
                    530:       IF( COLMAJOR ) THEN
                    531:          IF( WANTU1 .AND. P .GT. 0 ) THEN
                    532:             CALL ZLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
                    533:             CALL ZUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
                    534:      $                   LORGQRWORK, INFO)
                    535:          END IF
                    536:          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
                    537:             CALL ZLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
                    538:             CALL ZUNGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
                    539:      $                   WORK(IORGQR), LORGQRWORK, INFO )
                    540:          END IF
                    541:          IF( WANTV1T .AND. Q .GT. 0 ) THEN
                    542:             CALL ZLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2),
                    543:      $                   LDV1T )
                    544:             V1T(1, 1) = ONE
                    545:             DO J = 2, Q
                    546:                V1T(1,J) = ZERO
                    547:                V1T(J,1) = ZERO
                    548:             END DO
                    549:             CALL ZUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
                    550:      $                   WORK(IORGLQ), LORGLQWORK, INFO )
                    551:          END IF
                    552:          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
                    553:             CALL ZLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T )
                    554:             CALL ZLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
                    555:      $                   V2T(P+1,P+1), LDV2T )
                    556:             CALL ZUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
                    557:      $                   WORK(IORGLQ), LORGLQWORK, INFO )
                    558:          END IF
                    559:       ELSE
                    560:          IF( WANTU1 .AND. P .GT. 0 ) THEN
                    561:             CALL ZLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 )
                    562:             CALL ZUNGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ),
                    563:      $                   LORGLQWORK, INFO)
                    564:          END IF
                    565:          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
                    566:             CALL ZLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 )
                    567:             CALL ZUNGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
                    568:      $                   WORK(IORGLQ), LORGLQWORK, INFO )
                    569:          END IF
                    570:          IF( WANTV1T .AND. Q .GT. 0 ) THEN
                    571:             CALL ZLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2),
                    572:      $                   LDV1T )
                    573:             V1T(1, 1) = ONE
                    574:             DO J = 2, Q
                    575:                V1T(1,J) = ZERO
                    576:                V1T(J,1) = ZERO
                    577:             END DO
                    578:             CALL ZUNGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
                    579:      $                   WORK(IORGQR), LORGQRWORK, INFO )
                    580:          END IF
                    581:          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
                    582:             CALL ZLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T )
                    583:             CALL ZLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22,
                    584:      $                   V2T(P+1,P+1), LDV2T )
                    585:             CALL ZUNGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
                    586:      $                   WORK(IORGQR), LORGQRWORK, INFO )
                    587:          END IF
                    588:       END IF
                    589: *
                    590: *     Compute the CSD of the matrix in bidiagonal-block form
                    591: *
                    592:       CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA,
                    593:      $             RWORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
                    594:      $             LDV2T, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D),
                    595:      $             RWORK(IB12E), RWORK(IB21D), RWORK(IB21E),
                    596:      $             RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD),
                    597:      $             LBBCSDWORK, INFO )
                    598: *
                    599: *     Permute rows and columns to place identity submatrices in top-
                    600: *     left corner of (1,1)-block and/or bottom-right corner of (1,2)-
                    601: *     block and/or bottom-right corner of (2,1)-block and/or top-left
                    602: *     corner of (2,2)-block 
                    603: *
                    604:       IF( Q .GT. 0 .AND. WANTU2 ) THEN
                    605:          DO I = 1, Q
                    606:             IWORK(I) = M - P - Q + I
                    607:          END DO
                    608:          DO I = Q + 1, M - P
                    609:             IWORK(I) = I - Q
                    610:          END DO
                    611:          IF( COLMAJOR ) THEN
                    612:             CALL ZLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
                    613:          ELSE
                    614:             CALL ZLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK )
                    615:          END IF
                    616:       END IF
                    617:       IF( M .GT. 0 .AND. WANTV2T ) THEN
                    618:          DO I = 1, P
                    619:             IWORK(I) = M - P - Q + I
                    620:          END DO
                    621:          DO I = P + 1, M - Q
                    622:             IWORK(I) = I - P
                    623:          END DO
                    624:          IF( .NOT. COLMAJOR ) THEN
                    625:             CALL ZLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
                    626:          ELSE
                    627:             CALL ZLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
                    628:          END IF
                    629:       END IF
                    630: *
                    631:       RETURN
                    632: *
                    633: *     End ZUNCSD
                    634: *
                    635:       END
                    636: 

CVSweb interface <joel.bertrand@systella.fr>