LAPACK: SRC/zunbdb.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
00002      $                   X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
00003      $                   TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
00004       IMPLICIT NONE
00005 *
00006 *  -- LAPACK routine ((version 3.3.0)) --
00007 *
00008 *  -- Contributed by Brian Sutton of the Randolph-Macon College --
00009 *  -- November 2010
00010 *
00011 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00012 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--     
00013 *
00014 *     .. Scalar Arguments ..
00015       CHARACTER          SIGNS, TRANS
00016       INTEGER            INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
00017      $                   Q
00018 *     ..
00019 *     .. Array Arguments ..
00020       DOUBLE PRECISION   PHI( * ), THETA( * )
00021       COMPLEX*16         TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
00022      $                   WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
00023      $                   X21( LDX21, * ), X22( LDX22, * )
00024 *     ..
00025 *
00026 *  Purpose
00027 *  =======
00028 *
00029 *  ZUNBDB simultaneously bidiagonalizes the blocks of an M-by-M
00030 *  partitioned unitary matrix X:
00031 *
00032 *                                  [ B11 | B12 0  0 ]
00033 *      [ X11 | X12 ]   [ P1 |    ] [  0  |  0 -I  0 ] [ Q1 |    ]**H
00034 *  X = [-----------] = [---------] [----------------] [---------]   .
00035 *      [ X21 | X22 ]   [    | P2 ] [ B21 | B22 0  0 ] [    | Q2 ]
00036 *                                  [  0  |  0  0  I ]
00037 *
00038 *  X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
00039 *  not the case, then X must be transposed and/or permuted. This can be
00040 *  done in constant time using the TRANS and SIGNS options. See ZUNCSD
00041 *  for details.)
00042 *
00043 *  The unitary matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
00044 *  (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
00045 *  represented implicitly by Householder vectors.
00046 *
00047 *  B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
00048 *  implicitly by angles THETA, PHI.
00049 *
00050 *  Arguments
00051 *  =========
00052 *
00053 *  TRANS   (input) CHARACTER
00054 *          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
00055 *                      order;
00056 *          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
00057 *                      major order.
00058 *
00059 *  SIGNS   (input) CHARACTER
00060 *          = 'O':      The lower-left block is made nonpositive (the
00061 *                      "other" convention);
00062 *          otherwise:  The upper-right block is made nonpositive (the
00063 *                      "default" convention).
00064 *
00065 *  M       (input) INTEGER
00066 *          The number of rows and columns in X.
00067 *
00068 *  P       (input) INTEGER
00069 *          The number of rows in X11 and X12. 0 <= P <= M.
00070 *
00071 *  Q       (input) INTEGER
00072 *          The number of columns in X11 and X21. 0 <= Q <=
00073 *          MIN(P,M-P,M-Q).
00074 *
00075 *  X11     (input/output) COMPLEX*16 array, dimension (LDX11,Q)
00076 *          On entry, the top-left block of the unitary matrix to be
00077 *          reduced. On exit, the form depends on TRANS:
00078 *          If TRANS = 'N', then
00079 *             the columns of tril(X11) specify reflectors for P1,
00080 *             the rows of triu(X11,1) specify reflectors for Q1;
00081 *          else TRANS = 'T', and
00082 *             the rows of triu(X11) specify reflectors for P1,
00083 *             the columns of tril(X11,-1) specify reflectors for Q1.
00084 *
00085 *  LDX11   (input) INTEGER
00086 *          The leading dimension of X11. If TRANS = 'N', then LDX11 >=
00087 *          P; else LDX11 >= Q.
00088 *
00089 *  X12     (input/output) COMPLEX*16 array, dimension (LDX12,M-Q)
00090 *          On entry, the top-right block of the unitary matrix to
00091 *          be reduced. On exit, the form depends on TRANS:
00092 *          If TRANS = 'N', then
00093 *             the rows of triu(X12) specify the first P reflectors for
00094 *             Q2;
00095 *          else TRANS = 'T', and
00096 *             the columns of tril(X12) specify the first P reflectors
00097 *             for Q2.
00098 *
00099 *  LDX12   (input) INTEGER
00100 *          The leading dimension of X12. If TRANS = 'N', then LDX12 >=
00101 *          P; else LDX11 >= M-Q.
00102 *
00103 *  X21     (input/output) COMPLEX*16 array, dimension (LDX21,Q)
00104 *          On entry, the bottom-left block of the unitary matrix to
00105 *          be reduced. On exit, the form depends on TRANS:
00106 *          If TRANS = 'N', then
00107 *             the columns of tril(X21) specify reflectors for P2;
00108 *          else TRANS = 'T', and
00109 *             the rows of triu(X21) specify reflectors for P2.
00110 *
00111 *  LDX21   (input) INTEGER
00112 *          The leading dimension of X21. If TRANS = 'N', then LDX21 >=
00113 *          M-P; else LDX21 >= Q.
00114 *
00115 *  X22     (input/output) COMPLEX*16 array, dimension (LDX22,M-Q)
00116 *          On entry, the bottom-right block of the unitary matrix to
00117 *          be reduced. On exit, the form depends on TRANS:
00118 *          If TRANS = 'N', then
00119 *             the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
00120 *             M-P-Q reflectors for Q2,
00121 *          else TRANS = 'T', and
00122 *             the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
00123 *             M-P-Q reflectors for P2.
00124 *
00125 *  LDX22   (input) INTEGER
00126 *          The leading dimension of X22. If TRANS = 'N', then LDX22 >=
00127 *          M-P; else LDX22 >= M-Q.
00128 *
00129 *  THETA   (output) DOUBLE PRECISION array, dimension (Q)
00130 *          The entries of the bidiagonal blocks B11, B12, B21, B22 can
00131 *          be computed from the angles THETA and PHI. See Further
00132 *          Details.
00133 *
00134 *  PHI     (output) DOUBLE PRECISION array, dimension (Q-1)
00135 *          The entries of the bidiagonal blocks B11, B12, B21, B22 can
00136 *          be computed from the angles THETA and PHI. See Further
00137 *          Details.
00138 *
00139 *  TAUP1   (output) COMPLEX*16 array, dimension (P)
00140 *          The scalar factors of the elementary reflectors that define
00141 *          P1.
00142 *
00143 *  TAUP2   (output) COMPLEX*16 array, dimension (M-P)
00144 *          The scalar factors of the elementary reflectors that define
00145 *          P2.
00146 *
00147 *  TAUQ1   (output) COMPLEX*16 array, dimension (Q)
00148 *          The scalar factors of the elementary reflectors that define
00149 *          Q1.
00150 *
00151 *  TAUQ2   (output) COMPLEX*16 array, dimension (M-Q)
00152 *          The scalar factors of the elementary reflectors that define
00153 *          Q2.
00154 *
00155 *  WORK    (workspace) COMPLEX*16 array, dimension (LWORK)
00156 *
00157 *  LWORK   (input) INTEGER
00158 *          The dimension of the array WORK. LWORK >= M-Q.
00159 *
00160 *          If LWORK = -1, then a workspace query is assumed; the routine
00161 *          only calculates the optimal size of the WORK array, returns
00162 *          this value as the first entry of the WORK array, and no error
00163 *          message related to LWORK is issued by XERBLA.
00164 *
00165 *  INFO    (output) INTEGER
00166 *          = 0:  successful exit.
00167 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00168 *
00169 *  Further Details
00170 *  ===============
00171 *
00172 *  The bidiagonal blocks B11, B12, B21, and B22 are represented
00173 *  implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
00174 *  PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
00175 *  lower bidiagonal. Every entry in each bidiagonal band is a product
00176 *  of a sine or cosine of a THETA with a sine or cosine of a PHI. See
00177 *  [1] or ZUNCSD for details.
00178 *
00179 *  P1, P2, Q1, and Q2 are represented as products of elementary
00180 *  reflectors. See ZUNCSD for details on generating P1, P2, Q1, and Q2
00181 *  using ZUNGQR and ZUNGLQ.
00182 *
00183 *  Reference
00184 *  =========
00185 *
00186 *  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
00187 *      Algorithms, 50(1):33-65, 2009.
00188 *
00189 *  ====================================================================
00190 *
00191 *     .. Parameters ..
00192       DOUBLE PRECISION   REALONE
00193       PARAMETER          ( REALONE = 1.0D0 )
00194       COMPLEX*16         NEGONE, ONE
00195       PARAMETER          ( NEGONE = (-1.0D0,0.0D0),
00196      $                     ONE = (1.0D0,0.0D0) )
00197 *     ..
00198 *     .. Local Scalars ..
00199       LOGICAL            COLMAJOR, LQUERY
00200       INTEGER            I, LWORKMIN, LWORKOPT
00201       DOUBLE PRECISION   Z1, Z2, Z3, Z4
00202 *     ..
00203 *     .. External Subroutines ..
00204       EXTERNAL           ZAXPY, ZLARF, ZLARFGP, ZSCAL, XERBLA
00205       EXTERNAL           ZLACGV
00206 *
00207 *     ..
00208 *     .. External Functions ..
00209       DOUBLE PRECISION   DZNRM2
00210       LOGICAL            LSAME
00211       EXTERNAL           DZNRM2, LSAME
00212 *     ..
00213 *     .. Intrinsic Functions
00214       INTRINSIC          ATAN2, COS, MAX, MIN, SIN
00215       INTRINSIC          DCMPLX, DCONJG
00216 *     ..
00217 *     .. Executable Statements ..
00218 *
00219 *     Test input arguments
00220 *
00221       INFO = 0
00222       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
00223       IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN
00224          Z1 = REALONE
00225          Z2 = REALONE
00226          Z3 = REALONE
00227          Z4 = REALONE
00228       ELSE
00229          Z1 = REALONE
00230          Z2 = -REALONE
00231          Z3 = REALONE
00232          Z4 = -REALONE
00233       END IF
00234       LQUERY = LWORK .EQ. -1
00235 *
00236       IF( M .LT. 0 ) THEN
00237          INFO = -3
00238       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
00239          INFO = -4
00240       ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR.
00241      $         Q .GT. M-Q ) THEN
00242          INFO = -5
00243       ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN
00244          INFO = -7
00245       ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
00246          INFO = -7
00247       ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
00248          INFO = -9
00249       ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
00250          INFO = -9
00251       ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
00252          INFO = -11
00253       ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
00254          INFO = -11
00255       ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
00256          INFO = -13
00257       ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
00258          INFO = -13
00259       END IF
00260 *
00261 *     Compute workspace
00262 *
00263       IF( INFO .EQ. 0 ) THEN
00264          LWORKOPT = M - Q
00265          LWORKMIN = M - Q
00266          WORK(1) = LWORKOPT
00267          IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
00268             INFO = -21
00269          END IF
00270       END IF
00271       IF( INFO .NE. 0 ) THEN
00272          CALL XERBLA( 'xORBDB', -INFO )
00273          RETURN
00274       ELSE IF( LQUERY ) THEN
00275          RETURN
00276       END IF
00277 *
00278 *     Handle column-major and row-major separately
00279 *
00280       IF( COLMAJOR ) THEN
00281 *
00282 *        Reduce columns 1, ..., Q of X11, X12, X21, and X22 
00283 *
00284          DO I = 1, Q
00285 *
00286             IF( I .EQ. 1 ) THEN
00287                CALL ZSCAL( P-I+1, DCMPLX( Z1, 0.0D0 ), X11(I,I), 1 )
00288             ELSE
00289                CALL ZSCAL( P-I+1, DCMPLX( Z1*COS(PHI(I-1)), 0.0D0 ),
00290      $                     X11(I,I), 1 )
00291                CALL ZAXPY( P-I+1, DCMPLX( -Z1*Z3*Z4*SIN(PHI(I-1)),
00292      $                     0.0D0 ), X12(I,I-1), 1, X11(I,I), 1 )
00293             END IF
00294             IF( I .EQ. 1 ) THEN
00295                CALL ZSCAL( M-P-I+1, DCMPLX( Z2, 0.0D0 ), X21(I,I), 1 )
00296             ELSE
00297                CALL ZSCAL( M-P-I+1, DCMPLX( Z2*COS(PHI(I-1)), 0.0D0 ),
00298      $                     X21(I,I), 1 )
00299                CALL ZAXPY( M-P-I+1, DCMPLX( -Z2*Z3*Z4*SIN(PHI(I-1)),
00300      $                     0.0D0 ), X22(I,I-1), 1, X21(I,I), 1 )
00301             END IF
00302 *
00303             THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), 1 ),
00304      $                 DZNRM2( P-I+1, X11(I,I), 1 ) )
00305 *
00306             CALL ZLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
00307             X11(I,I) = ONE
00308             CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
00309             X21(I,I) = ONE
00310 *
00311             CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, DCONJG(TAUP1(I)),
00312      $                  X11(I,I+1), LDX11, WORK )
00313             CALL ZLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1,
00314      $                  DCONJG(TAUP1(I)), X12(I,I), LDX12, WORK )
00315             CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1,
00316      $                  DCONJG(TAUP2(I)), X21(I,I+1), LDX21, WORK )
00317             CALL ZLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1,
00318      $                  DCONJG(TAUP2(I)), X22(I,I), LDX22, WORK )
00319 *
00320             IF( I .LT. Q ) THEN
00321                CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ),
00322      $                     X11(I,I+1), LDX11 )
00323                CALL ZAXPY( Q-I, DCMPLX( Z2*Z3*COS(THETA(I)), 0.0D0 ),
00324      $                     X21(I,I+1), LDX21, X11(I,I+1), LDX11 )
00325             END IF
00326             CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4*SIN(THETA(I)), 0.0D0 ),
00327      $                  X12(I,I), LDX12 )
00328             CALL ZAXPY( M-Q-I+1, DCMPLX( Z2*Z4*COS(THETA(I)), 0.0D0 ),
00329      $                  X22(I,I), LDX22, X12(I,I), LDX12 )
00330 *
00331             IF( I .LT. Q )
00332      $         PHI(I) = ATAN2( DZNRM2( Q-I, X11(I,I+1), LDX11 ),
00333      $                  DZNRM2( M-Q-I+1, X12(I,I), LDX12 ) )
00334 *
00335             IF( I .LT. Q ) THEN
00336                CALL ZLACGV( Q-I, X11(I,I+1), LDX11 )
00337                CALL ZLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
00338      $                       TAUQ1(I) )
00339                X11(I,I+1) = ONE
00340             END IF
00341             CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
00342             CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
00343      $                    TAUQ2(I) )
00344             X12(I,I) = ONE
00345 *
00346             IF( I .LT. Q ) THEN
00347                CALL ZLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
00348      $                     X11(I+1,I+1), LDX11, WORK )
00349                CALL ZLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
00350      $                     X21(I+1,I+1), LDX21, WORK )
00351             END IF
00352             CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
00353      $                  X12(I+1,I), LDX12, WORK )
00354             CALL ZLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
00355      $                  X22(I+1,I), LDX22, WORK )
00356 *
00357             IF( I .LT. Q )
00358      $         CALL ZLACGV( Q-I, X11(I,I+1), LDX11 )
00359             CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
00360 *
00361          END DO
00362 *
00363 *        Reduce columns Q + 1, ..., P of X12, X22
00364 *
00365          DO I = Q + 1, P
00366 *
00367             CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I),
00368      $                  LDX12 )
00369             CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
00370             CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
00371      $                    TAUQ2(I) )
00372             X12(I,I) = ONE
00373 *
00374             CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
00375      $                  X12(I+1,I), LDX12, WORK )
00376             IF( M-P-Q .GE. 1 )
00377      $         CALL ZLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
00378      $                     TAUQ2(I), X22(Q+1,I), LDX22, WORK )
00379 *
00380             CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
00381 *
00382          END DO
00383 *
00384 *        Reduce columns P + 1, ..., M - Q of X12, X22
00385 *
00386          DO I = 1, M - P - Q
00387 *
00388             CALL ZSCAL( M-P-Q-I+1, DCMPLX( Z2*Z4, 0.0D0 ),
00389      $                  X22(Q+I,P+I), LDX22 )
00390             CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 )
00391             CALL ZLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
00392      $                    LDX22, TAUQ2(P+I) )
00393             X22(Q+I,P+I) = ONE
00394             CALL ZLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
00395      $                  TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
00396 *
00397             CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 )
00398 *
00399          END DO
00400 *
00401       ELSE
00402 *
00403 *        Reduce columns 1, ..., Q of X11, X12, X21, X22
00404 *
00405          DO I = 1, Q
00406 *
00407             IF( I .EQ. 1 ) THEN
00408                CALL ZSCAL( P-I+1, DCMPLX( Z1, 0.0D0 ), X11(I,I),
00409      $                     LDX11 )
00410             ELSE
00411                CALL ZSCAL( P-I+1, DCMPLX( Z1*COS(PHI(I-1)), 0.0D0 ),
00412      $                     X11(I,I), LDX11 )
00413                CALL ZAXPY( P-I+1, DCMPLX( -Z1*Z3*Z4*SIN(PHI(I-1)),
00414      $                     0.0D0 ), X12(I-1,I), LDX12, X11(I,I), LDX11 )
00415             END IF
00416             IF( I .EQ. 1 ) THEN
00417                CALL ZSCAL( M-P-I+1, DCMPLX( Z2, 0.0D0 ), X21(I,I),
00418      $                     LDX21 )
00419             ELSE
00420                CALL ZSCAL( M-P-I+1, DCMPLX( Z2*COS(PHI(I-1)), 0.0D0 ),
00421      $                     X21(I,I), LDX21 )
00422                CALL ZAXPY( M-P-I+1, DCMPLX( -Z2*Z3*Z4*SIN(PHI(I-1)),
00423      $                     0.0D0 ), X22(I-1,I), LDX22, X21(I,I), LDX21 )
00424             END IF
00425 *
00426             THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), LDX21 ),
00427      $                 DZNRM2( P-I+1, X11(I,I), LDX11 ) )
00428 *
00429             CALL ZLACGV( P-I+1, X11(I,I), LDX11 )
00430             CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 )
00431 *
00432             CALL ZLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
00433             X11(I,I) = ONE
00434             CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
00435      $                    TAUP2(I) )
00436             X21(I,I) = ONE
00437 *
00438             CALL ZLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
00439      $                  X11(I+1,I), LDX11, WORK )
00440             CALL ZLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I),
00441      $                  X12(I,I), LDX12, WORK )
00442             CALL ZLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
00443      $                  X21(I+1,I), LDX21, WORK )
00444             CALL ZLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
00445      $                  TAUP2(I), X22(I,I), LDX22, WORK )
00446 *
00447             CALL ZLACGV( P-I+1, X11(I,I), LDX11 )
00448             CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 )
00449 *
00450             IF( I .LT. Q ) THEN
00451                CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ),
00452      $                     X11(I+1,I), 1 )
00453                CALL ZAXPY( Q-I, DCMPLX( Z2*Z3*COS(THETA(I)), 0.0D0 ),
00454      $                     X21(I+1,I), 1, X11(I+1,I), 1 )
00455             END IF
00456             CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4*SIN(THETA(I)), 0.0D0 ),
00457      $                  X12(I,I), 1 )
00458             CALL ZAXPY( M-Q-I+1, DCMPLX( Z2*Z4*COS(THETA(I)), 0.0D0 ),
00459      $                  X22(I,I), 1, X12(I,I), 1 )
00460 *
00461             IF( I .LT. Q )
00462      $         PHI(I) = ATAN2( DZNRM2( Q-I, X11(I+1,I), 1 ),
00463      $                  DZNRM2( M-Q-I+1, X12(I,I), 1 ) )
00464 *
00465             IF( I .LT. Q ) THEN
00466                CALL ZLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) )
00467                X11(I+1,I) = ONE
00468             END IF
00469             CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
00470             X12(I,I) = ONE
00471 *
00472             IF( I .LT. Q ) THEN
00473                CALL ZLARF( 'L', Q-I, P-I, X11(I+1,I), 1,
00474      $                     DCONJG(TAUQ1(I)), X11(I+1,I+1), LDX11, WORK )
00475                CALL ZLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1,
00476      $                     DCONJG(TAUQ1(I)), X21(I+1,I+1), LDX21, WORK )
00477             END IF
00478             CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1,
00479      $                  DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK )
00480             CALL ZLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1,
00481      $                  DCONJG(TAUQ2(I)), X22(I,I+1), LDX22, WORK )
00482 *
00483          END DO
00484 *
00485 *        Reduce columns Q + 1, ..., P of X12, X22
00486 *
00487          DO I = Q + 1, P
00488 *
00489             CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I), 1 )
00490             CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
00491             X12(I,I) = ONE
00492 *
00493             CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1,
00494      $                  DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK )
00495             IF( M-P-Q .GE. 1 )
00496      $         CALL ZLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1,
00497      $                     DCONJG(TAUQ2(I)), X22(I,Q+1), LDX22, WORK )
00498 *
00499          END DO
00500 *
00501 *        Reduce columns P + 1, ..., M - Q of X12, X22
00502 *
00503          DO I = 1, M - P - Q
00504 *
00505             CALL ZSCAL( M-P-Q-I+1, DCMPLX( Z2*Z4, 0.0D0 ),
00506      $                  X22(P+I,Q+I), 1 )
00507             CALL ZLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
00508      $                    TAUQ2(P+I) )
00509             X22(P+I,Q+I) = ONE
00510 *
00511             CALL ZLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
00512      $                  DCONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), LDX22,
00513      $                  WORK )
00514 *
00515          END DO
00516 *
00517       END IF
00518 *
00519       RETURN
00520 *
00521 *     End of ZUNBDB
00522 *
00523       END
00524