1 SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
2 $ X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
3 $ TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
4 IMPLICIT NONE
5 *
6 * -- LAPACK routine ((version 3.3.0)) --
7 *
8 * -- Contributed by Brian Sutton of the Randolph-Macon College --
9 * -- November 2010
10 *
11 * -- LAPACK is a software package provided by Univ. of Tennessee, --
12 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
13 *
14 * .. Scalar Arguments ..
15 CHARACTER SIGNS, TRANS
16 INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
17 $ Q
18 * ..
19 * .. Array Arguments ..
20 DOUBLE PRECISION PHI( * ), THETA( * )
21 COMPLEX*16 TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
22 $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
23 $ X21( LDX21, * ), X22( LDX22, * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * ZUNBDB simultaneously bidiagonalizes the blocks of an M-by-M
30 * partitioned unitary matrix X:
31 *
32 * [ B11 | B12 0 0 ]
33 * [ X11 | X12 ] [ P1 | ] [ 0 | 0 -I 0 ] [ Q1 | ]**H
34 * X = [-----------] = [---------] [----------------] [---------] .
35 * [ X21 | X22 ] [ | P2 ] [ B21 | B22 0 0 ] [ | Q2 ]
36 * [ 0 | 0 0 I ]
37 *
38 * X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
39 * not the case, then X must be transposed and/or permuted. This can be
40 * done in constant time using the TRANS and SIGNS options. See ZUNCSD
41 * for details.)
42 *
43 * The unitary matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
44 * (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
45 * represented implicitly by Householder vectors.
46 *
47 * B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
48 * implicitly by angles THETA, PHI.
49 *
50 * Arguments
51 * =========
52 *
53 * TRANS (input) CHARACTER
54 * = 'T': X, U1, U2, V1T, and V2T are stored in row-major
55 * order;
56 * otherwise: X, U1, U2, V1T, and V2T are stored in column-
57 * major order.
58 *
59 * SIGNS (input) CHARACTER
60 * = 'O': The lower-left block is made nonpositive (the
61 * "other" convention);
62 * otherwise: The upper-right block is made nonpositive (the
63 * "default" convention).
64 *
65 * M (input) INTEGER
66 * The number of rows and columns in X.
67 *
68 * P (input) INTEGER
69 * The number of rows in X11 and X12. 0 <= P <= M.
70 *
71 * Q (input) INTEGER
72 * The number of columns in X11 and X21. 0 <= Q <=
73 * MIN(P,M-P,M-Q).
74 *
75 * X11 (input/output) COMPLEX*16 array, dimension (LDX11,Q)
76 * On entry, the top-left block of the unitary matrix to be
77 * reduced. On exit, the form depends on TRANS:
78 * If TRANS = 'N', then
79 * the columns of tril(X11) specify reflectors for P1,
80 * the rows of triu(X11,1) specify reflectors for Q1;
81 * else TRANS = 'T', and
82 * the rows of triu(X11) specify reflectors for P1,
83 * the columns of tril(X11,-1) specify reflectors for Q1.
84 *
85 * LDX11 (input) INTEGER
86 * The leading dimension of X11. If TRANS = 'N', then LDX11 >=
87 * P; else LDX11 >= Q.
88 *
89 * X12 (input/output) COMPLEX*16 array, dimension (LDX12,M-Q)
90 * On entry, the top-right block of the unitary matrix to
91 * be reduced. On exit, the form depends on TRANS:
92 * If TRANS = 'N', then
93 * the rows of triu(X12) specify the first P reflectors for
94 * Q2;
95 * else TRANS = 'T', and
96 * the columns of tril(X12) specify the first P reflectors
97 * for Q2.
98 *
99 * LDX12 (input) INTEGER
100 * The leading dimension of X12. If TRANS = 'N', then LDX12 >=
101 * P; else LDX11 >= M-Q.
102 *
103 * X21 (input/output) COMPLEX*16 array, dimension (LDX21,Q)
104 * On entry, the bottom-left block of the unitary matrix to
105 * be reduced. On exit, the form depends on TRANS:
106 * If TRANS = 'N', then
107 * the columns of tril(X21) specify reflectors for P2;
108 * else TRANS = 'T', and
109 * the rows of triu(X21) specify reflectors for P2.
110 *
111 * LDX21 (input) INTEGER
112 * The leading dimension of X21. If TRANS = 'N', then LDX21 >=
113 * M-P; else LDX21 >= Q.
114 *
115 * X22 (input/output) COMPLEX*16 array, dimension (LDX22,M-Q)
116 * On entry, the bottom-right block of the unitary matrix to
117 * be reduced. On exit, the form depends on TRANS:
118 * If TRANS = 'N', then
119 * the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
120 * M-P-Q reflectors for Q2,
121 * else TRANS = 'T', and
122 * the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
123 * M-P-Q reflectors for P2.
124 *
125 * LDX22 (input) INTEGER
126 * The leading dimension of X22. If TRANS = 'N', then LDX22 >=
127 * M-P; else LDX22 >= M-Q.
128 *
129 * THETA (output) DOUBLE PRECISION array, dimension (Q)
130 * The entries of the bidiagonal blocks B11, B12, B21, B22 can
131 * be computed from the angles THETA and PHI. See Further
132 * Details.
133 *
134 * PHI (output) DOUBLE PRECISION array, dimension (Q-1)
135 * The entries of the bidiagonal blocks B11, B12, B21, B22 can
136 * be computed from the angles THETA and PHI. See Further
137 * Details.
138 *
139 * TAUP1 (output) COMPLEX*16 array, dimension (P)
140 * The scalar factors of the elementary reflectors that define
141 * P1.
142 *
143 * TAUP2 (output) COMPLEX*16 array, dimension (M-P)
144 * The scalar factors of the elementary reflectors that define
145 * P2.
146 *
147 * TAUQ1 (output) COMPLEX*16 array, dimension (Q)
148 * The scalar factors of the elementary reflectors that define
149 * Q1.
150 *
151 * TAUQ2 (output) COMPLEX*16 array, dimension (M-Q)
152 * The scalar factors of the elementary reflectors that define
153 * Q2.
154 *
155 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
156 *
157 * LWORK (input) INTEGER
158 * The dimension of the array WORK. LWORK >= M-Q.
159 *
160 * If LWORK = -1, then a workspace query is assumed; the routine
161 * only calculates the optimal size of the WORK array, returns
162 * this value as the first entry of the WORK array, and no error
163 * message related to LWORK is issued by XERBLA.
164 *
165 * INFO (output) INTEGER
166 * = 0: successful exit.
167 * < 0: if INFO = -i, the i-th argument had an illegal value.
168 *
169 * Further Details
170 * ===============
171 *
172 * The bidiagonal blocks B11, B12, B21, and B22 are represented
173 * implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
174 * PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
175 * lower bidiagonal. Every entry in each bidiagonal band is a product
176 * of a sine or cosine of a THETA with a sine or cosine of a PHI. See
177 * [1] or ZUNCSD for details.
178 *
179 * P1, P2, Q1, and Q2 are represented as products of elementary
180 * reflectors. See ZUNCSD for details on generating P1, P2, Q1, and Q2
181 * using ZUNGQR and ZUNGLQ.
182 *
183 * Reference
184 * =========
185 *
186 * [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
187 * Algorithms, 50(1):33-65, 2009.
188 *
189 * ====================================================================
190 *
191 * .. Parameters ..
192 DOUBLE PRECISION REALONE
193 PARAMETER ( REALONE = 1.0D0 )
194 COMPLEX*16 NEGONE, ONE
195 PARAMETER ( NEGONE = (-1.0D0,0.0D0),
196 $ ONE = (1.0D0,0.0D0) )
197 * ..
198 * .. Local Scalars ..
199 LOGICAL COLMAJOR, LQUERY
200 INTEGER I, LWORKMIN, LWORKOPT
201 DOUBLE PRECISION Z1, Z2, Z3, Z4
202 * ..
203 * .. External Subroutines ..
204 EXTERNAL ZAXPY, ZLARF, ZLARFGP, ZSCAL, XERBLA
205 EXTERNAL ZLACGV
206 *
207 * ..
208 * .. External Functions ..
209 DOUBLE PRECISION DZNRM2
210 LOGICAL LSAME
211 EXTERNAL DZNRM2, LSAME
212 * ..
213 * .. Intrinsic Functions
214 INTRINSIC ATAN2, COS, MAX, MIN, SIN
215 INTRINSIC DCMPLX, DCONJG
216 * ..
217 * .. Executable Statements ..
218 *
219 * Test input arguments
220 *
221 INFO = 0
222 COLMAJOR = .NOT. LSAME( TRANS, 'T' )
223 IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN
224 Z1 = REALONE
225 Z2 = REALONE
226 Z3 = REALONE
227 Z4 = REALONE
228 ELSE
229 Z1 = REALONE
230 Z2 = -REALONE
231 Z3 = REALONE
232 Z4 = -REALONE
233 END IF
234 LQUERY = LWORK .EQ. -1
235 *
236 IF( M .LT. 0 ) THEN
237 INFO = -3
238 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
239 INFO = -4
240 ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR.
241 $ Q .GT. M-Q ) THEN
242 INFO = -5
243 ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN
244 INFO = -7
245 ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
246 INFO = -7
247 ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
248 INFO = -9
249 ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
250 INFO = -9
251 ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
252 INFO = -11
253 ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
254 INFO = -11
255 ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
256 INFO = -13
257 ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
258 INFO = -13
259 END IF
260 *
261 * Compute workspace
262 *
263 IF( INFO .EQ. 0 ) THEN
264 LWORKOPT = M - Q
265 LWORKMIN = M - Q
266 WORK(1) = LWORKOPT
267 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
268 INFO = -21
269 END IF
270 END IF
271 IF( INFO .NE. 0 ) THEN
272 CALL XERBLA( 'xORBDB', -INFO )
273 RETURN
274 ELSE IF( LQUERY ) THEN
275 RETURN
276 END IF
277 *
278 * Handle column-major and row-major separately
279 *
280 IF( COLMAJOR ) THEN
281 *
282 * Reduce columns 1, ..., Q of X11, X12, X21, and X22
283 *
284 DO I = 1, Q
285 *
286 IF( I .EQ. 1 ) THEN
287 CALL ZSCAL( P-I+1, DCMPLX( Z1, 0.0D0 ), X11(I,I), 1 )
288 ELSE
289 CALL ZSCAL( P-I+1, DCMPLX( Z1*COS(PHI(I-1)), 0.0D0 ),
290 $ X11(I,I), 1 )
291 CALL ZAXPY( P-I+1, DCMPLX( -Z1*Z3*Z4*SIN(PHI(I-1)),
292 $ 0.0D0 ), X12(I,I-1), 1, X11(I,I), 1 )
293 END IF
294 IF( I .EQ. 1 ) THEN
295 CALL ZSCAL( M-P-I+1, DCMPLX( Z2, 0.0D0 ), X21(I,I), 1 )
296 ELSE
297 CALL ZSCAL( M-P-I+1, DCMPLX( Z2*COS(PHI(I-1)), 0.0D0 ),
298 $ X21(I,I), 1 )
299 CALL ZAXPY( M-P-I+1, DCMPLX( -Z2*Z3*Z4*SIN(PHI(I-1)),
300 $ 0.0D0 ), X22(I,I-1), 1, X21(I,I), 1 )
301 END IF
302 *
303 THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), 1 ),
304 $ DZNRM2( P-I+1, X11(I,I), 1 ) )
305 *
306 CALL ZLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
307 X11(I,I) = ONE
308 CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
309 X21(I,I) = ONE
310 *
311 CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, DCONJG(TAUP1(I)),
312 $ X11(I,I+1), LDX11, WORK )
313 CALL ZLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1,
314 $ DCONJG(TAUP1(I)), X12(I,I), LDX12, WORK )
315 CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1,
316 $ DCONJG(TAUP2(I)), X21(I,I+1), LDX21, WORK )
317 CALL ZLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1,
318 $ DCONJG(TAUP2(I)), X22(I,I), LDX22, WORK )
319 *
320 IF( I .LT. Q ) THEN
321 CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ),
322 $ X11(I,I+1), LDX11 )
323 CALL ZAXPY( Q-I, DCMPLX( Z2*Z3*COS(THETA(I)), 0.0D0 ),
324 $ X21(I,I+1), LDX21, X11(I,I+1), LDX11 )
325 END IF
326 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4*SIN(THETA(I)), 0.0D0 ),
327 $ X12(I,I), LDX12 )
328 CALL ZAXPY( M-Q-I+1, DCMPLX( Z2*Z4*COS(THETA(I)), 0.0D0 ),
329 $ X22(I,I), LDX22, X12(I,I), LDX12 )
330 *
331 IF( I .LT. Q )
332 $ PHI(I) = ATAN2( DZNRM2( Q-I, X11(I,I+1), LDX11 ),
333 $ DZNRM2( M-Q-I+1, X12(I,I), LDX12 ) )
334 *
335 IF( I .LT. Q ) THEN
336 CALL ZLACGV( Q-I, X11(I,I+1), LDX11 )
337 CALL ZLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
338 $ TAUQ1(I) )
339 X11(I,I+1) = ONE
340 END IF
341 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
342 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
343 $ TAUQ2(I) )
344 X12(I,I) = ONE
345 *
346 IF( I .LT. Q ) THEN
347 CALL ZLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
348 $ X11(I+1,I+1), LDX11, WORK )
349 CALL ZLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
350 $ X21(I+1,I+1), LDX21, WORK )
351 END IF
352 CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
353 $ X12(I+1,I), LDX12, WORK )
354 CALL ZLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
355 $ X22(I+1,I), LDX22, WORK )
356 *
357 IF( I .LT. Q )
358 $ CALL ZLACGV( Q-I, X11(I,I+1), LDX11 )
359 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
360 *
361 END DO
362 *
363 * Reduce columns Q + 1, ..., P of X12, X22
364 *
365 DO I = Q + 1, P
366 *
367 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I),
368 $ LDX12 )
369 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
370 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
371 $ TAUQ2(I) )
372 X12(I,I) = ONE
373 *
374 CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
375 $ X12(I+1,I), LDX12, WORK )
376 IF( M-P-Q .GE. 1 )
377 $ CALL ZLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
378 $ TAUQ2(I), X22(Q+1,I), LDX22, WORK )
379 *
380 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
381 *
382 END DO
383 *
384 * Reduce columns P + 1, ..., M - Q of X12, X22
385 *
386 DO I = 1, M - P - Q
387 *
388 CALL ZSCAL( M-P-Q-I+1, DCMPLX( Z2*Z4, 0.0D0 ),
389 $ X22(Q+I,P+I), LDX22 )
390 CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 )
391 CALL ZLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
392 $ LDX22, TAUQ2(P+I) )
393 X22(Q+I,P+I) = ONE
394 CALL ZLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
395 $ TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
396 *
397 CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 )
398 *
399 END DO
400 *
401 ELSE
402 *
403 * Reduce columns 1, ..., Q of X11, X12, X21, X22
404 *
405 DO I = 1, Q
406 *
407 IF( I .EQ. 1 ) THEN
408 CALL ZSCAL( P-I+1, DCMPLX( Z1, 0.0D0 ), X11(I,I),
409 $ LDX11 )
410 ELSE
411 CALL ZSCAL( P-I+1, DCMPLX( Z1*COS(PHI(I-1)), 0.0D0 ),
412 $ X11(I,I), LDX11 )
413 CALL ZAXPY( P-I+1, DCMPLX( -Z1*Z3*Z4*SIN(PHI(I-1)),
414 $ 0.0D0 ), X12(I-1,I), LDX12, X11(I,I), LDX11 )
415 END IF
416 IF( I .EQ. 1 ) THEN
417 CALL ZSCAL( M-P-I+1, DCMPLX( Z2, 0.0D0 ), X21(I,I),
418 $ LDX21 )
419 ELSE
420 CALL ZSCAL( M-P-I+1, DCMPLX( Z2*COS(PHI(I-1)), 0.0D0 ),
421 $ X21(I,I), LDX21 )
422 CALL ZAXPY( M-P-I+1, DCMPLX( -Z2*Z3*Z4*SIN(PHI(I-1)),
423 $ 0.0D0 ), X22(I-1,I), LDX22, X21(I,I), LDX21 )
424 END IF
425 *
426 THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), LDX21 ),
427 $ DZNRM2( P-I+1, X11(I,I), LDX11 ) )
428 *
429 CALL ZLACGV( P-I+1, X11(I,I), LDX11 )
430 CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 )
431 *
432 CALL ZLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
433 X11(I,I) = ONE
434 CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
435 $ TAUP2(I) )
436 X21(I,I) = ONE
437 *
438 CALL ZLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
439 $ X11(I+1,I), LDX11, WORK )
440 CALL ZLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I),
441 $ X12(I,I), LDX12, WORK )
442 CALL ZLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
443 $ X21(I+1,I), LDX21, WORK )
444 CALL ZLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
445 $ TAUP2(I), X22(I,I), LDX22, WORK )
446 *
447 CALL ZLACGV( P-I+1, X11(I,I), LDX11 )
448 CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 )
449 *
450 IF( I .LT. Q ) THEN
451 CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ),
452 $ X11(I+1,I), 1 )
453 CALL ZAXPY( Q-I, DCMPLX( Z2*Z3*COS(THETA(I)), 0.0D0 ),
454 $ X21(I+1,I), 1, X11(I+1,I), 1 )
455 END IF
456 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4*SIN(THETA(I)), 0.0D0 ),
457 $ X12(I,I), 1 )
458 CALL ZAXPY( M-Q-I+1, DCMPLX( Z2*Z4*COS(THETA(I)), 0.0D0 ),
459 $ X22(I,I), 1, X12(I,I), 1 )
460 *
461 IF( I .LT. Q )
462 $ PHI(I) = ATAN2( DZNRM2( Q-I, X11(I+1,I), 1 ),
463 $ DZNRM2( M-Q-I+1, X12(I,I), 1 ) )
464 *
465 IF( I .LT. Q ) THEN
466 CALL ZLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) )
467 X11(I+1,I) = ONE
468 END IF
469 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
470 X12(I,I) = ONE
471 *
472 IF( I .LT. Q ) THEN
473 CALL ZLARF( 'L', Q-I, P-I, X11(I+1,I), 1,
474 $ DCONJG(TAUQ1(I)), X11(I+1,I+1), LDX11, WORK )
475 CALL ZLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1,
476 $ DCONJG(TAUQ1(I)), X21(I+1,I+1), LDX21, WORK )
477 END IF
478 CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1,
479 $ DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK )
480 CALL ZLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1,
481 $ DCONJG(TAUQ2(I)), X22(I,I+1), LDX22, WORK )
482 *
483 END DO
484 *
485 * Reduce columns Q + 1, ..., P of X12, X22
486 *
487 DO I = Q + 1, P
488 *
489 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I), 1 )
490 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
491 X12(I,I) = ONE
492 *
493 CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1,
494 $ DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK )
495 IF( M-P-Q .GE. 1 )
496 $ CALL ZLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1,
497 $ DCONJG(TAUQ2(I)), X22(I,Q+1), LDX22, WORK )
498 *
499 END DO
500 *
501 * Reduce columns P + 1, ..., M - Q of X12, X22
502 *
503 DO I = 1, M - P - Q
504 *
505 CALL ZSCAL( M-P-Q-I+1, DCMPLX( Z2*Z4, 0.0D0 ),
506 $ X22(P+I,Q+I), 1 )
507 CALL ZLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
508 $ TAUQ2(P+I) )
509 X22(P+I,Q+I) = ONE
510 *
511 CALL ZLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
512 $ DCONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), LDX22,
513 $ WORK )
514 *
515 END DO
516 *
517 END IF
518 *
519 RETURN
520 *
521 * End of ZUNBDB
522 *
523 END
524
2 $ X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
3 $ TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
4 IMPLICIT NONE
5 *
6 * -- LAPACK routine ((version 3.3.0)) --
7 *
8 * -- Contributed by Brian Sutton of the Randolph-Macon College --
9 * -- November 2010
10 *
11 * -- LAPACK is a software package provided by Univ. of Tennessee, --
12 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
13 *
14 * .. Scalar Arguments ..
15 CHARACTER SIGNS, TRANS
16 INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
17 $ Q
18 * ..
19 * .. Array Arguments ..
20 DOUBLE PRECISION PHI( * ), THETA( * )
21 COMPLEX*16 TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
22 $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
23 $ X21( LDX21, * ), X22( LDX22, * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * ZUNBDB simultaneously bidiagonalizes the blocks of an M-by-M
30 * partitioned unitary matrix X:
31 *
32 * [ B11 | B12 0 0 ]
33 * [ X11 | X12 ] [ P1 | ] [ 0 | 0 -I 0 ] [ Q1 | ]**H
34 * X = [-----------] = [---------] [----------------] [---------] .
35 * [ X21 | X22 ] [ | P2 ] [ B21 | B22 0 0 ] [ | Q2 ]
36 * [ 0 | 0 0 I ]
37 *
38 * X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
39 * not the case, then X must be transposed and/or permuted. This can be
40 * done in constant time using the TRANS and SIGNS options. See ZUNCSD
41 * for details.)
42 *
43 * The unitary matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
44 * (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
45 * represented implicitly by Householder vectors.
46 *
47 * B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
48 * implicitly by angles THETA, PHI.
49 *
50 * Arguments
51 * =========
52 *
53 * TRANS (input) CHARACTER
54 * = 'T': X, U1, U2, V1T, and V2T are stored in row-major
55 * order;
56 * otherwise: X, U1, U2, V1T, and V2T are stored in column-
57 * major order.
58 *
59 * SIGNS (input) CHARACTER
60 * = 'O': The lower-left block is made nonpositive (the
61 * "other" convention);
62 * otherwise: The upper-right block is made nonpositive (the
63 * "default" convention).
64 *
65 * M (input) INTEGER
66 * The number of rows and columns in X.
67 *
68 * P (input) INTEGER
69 * The number of rows in X11 and X12. 0 <= P <= M.
70 *
71 * Q (input) INTEGER
72 * The number of columns in X11 and X21. 0 <= Q <=
73 * MIN(P,M-P,M-Q).
74 *
75 * X11 (input/output) COMPLEX*16 array, dimension (LDX11,Q)
76 * On entry, the top-left block of the unitary matrix to be
77 * reduced. On exit, the form depends on TRANS:
78 * If TRANS = 'N', then
79 * the columns of tril(X11) specify reflectors for P1,
80 * the rows of triu(X11,1) specify reflectors for Q1;
81 * else TRANS = 'T', and
82 * the rows of triu(X11) specify reflectors for P1,
83 * the columns of tril(X11,-1) specify reflectors for Q1.
84 *
85 * LDX11 (input) INTEGER
86 * The leading dimension of X11. If TRANS = 'N', then LDX11 >=
87 * P; else LDX11 >= Q.
88 *
89 * X12 (input/output) COMPLEX*16 array, dimension (LDX12,M-Q)
90 * On entry, the top-right block of the unitary matrix to
91 * be reduced. On exit, the form depends on TRANS:
92 * If TRANS = 'N', then
93 * the rows of triu(X12) specify the first P reflectors for
94 * Q2;
95 * else TRANS = 'T', and
96 * the columns of tril(X12) specify the first P reflectors
97 * for Q2.
98 *
99 * LDX12 (input) INTEGER
100 * The leading dimension of X12. If TRANS = 'N', then LDX12 >=
101 * P; else LDX11 >= M-Q.
102 *
103 * X21 (input/output) COMPLEX*16 array, dimension (LDX21,Q)
104 * On entry, the bottom-left block of the unitary matrix to
105 * be reduced. On exit, the form depends on TRANS:
106 * If TRANS = 'N', then
107 * the columns of tril(X21) specify reflectors for P2;
108 * else TRANS = 'T', and
109 * the rows of triu(X21) specify reflectors for P2.
110 *
111 * LDX21 (input) INTEGER
112 * The leading dimension of X21. If TRANS = 'N', then LDX21 >=
113 * M-P; else LDX21 >= Q.
114 *
115 * X22 (input/output) COMPLEX*16 array, dimension (LDX22,M-Q)
116 * On entry, the bottom-right block of the unitary matrix to
117 * be reduced. On exit, the form depends on TRANS:
118 * If TRANS = 'N', then
119 * the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
120 * M-P-Q reflectors for Q2,
121 * else TRANS = 'T', and
122 * the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
123 * M-P-Q reflectors for P2.
124 *
125 * LDX22 (input) INTEGER
126 * The leading dimension of X22. If TRANS = 'N', then LDX22 >=
127 * M-P; else LDX22 >= M-Q.
128 *
129 * THETA (output) DOUBLE PRECISION array, dimension (Q)
130 * The entries of the bidiagonal blocks B11, B12, B21, B22 can
131 * be computed from the angles THETA and PHI. See Further
132 * Details.
133 *
134 * PHI (output) DOUBLE PRECISION array, dimension (Q-1)
135 * The entries of the bidiagonal blocks B11, B12, B21, B22 can
136 * be computed from the angles THETA and PHI. See Further
137 * Details.
138 *
139 * TAUP1 (output) COMPLEX*16 array, dimension (P)
140 * The scalar factors of the elementary reflectors that define
141 * P1.
142 *
143 * TAUP2 (output) COMPLEX*16 array, dimension (M-P)
144 * The scalar factors of the elementary reflectors that define
145 * P2.
146 *
147 * TAUQ1 (output) COMPLEX*16 array, dimension (Q)
148 * The scalar factors of the elementary reflectors that define
149 * Q1.
150 *
151 * TAUQ2 (output) COMPLEX*16 array, dimension (M-Q)
152 * The scalar factors of the elementary reflectors that define
153 * Q2.
154 *
155 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
156 *
157 * LWORK (input) INTEGER
158 * The dimension of the array WORK. LWORK >= M-Q.
159 *
160 * If LWORK = -1, then a workspace query is assumed; the routine
161 * only calculates the optimal size of the WORK array, returns
162 * this value as the first entry of the WORK array, and no error
163 * message related to LWORK is issued by XERBLA.
164 *
165 * INFO (output) INTEGER
166 * = 0: successful exit.
167 * < 0: if INFO = -i, the i-th argument had an illegal value.
168 *
169 * Further Details
170 * ===============
171 *
172 * The bidiagonal blocks B11, B12, B21, and B22 are represented
173 * implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
174 * PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
175 * lower bidiagonal. Every entry in each bidiagonal band is a product
176 * of a sine or cosine of a THETA with a sine or cosine of a PHI. See
177 * [1] or ZUNCSD for details.
178 *
179 * P1, P2, Q1, and Q2 are represented as products of elementary
180 * reflectors. See ZUNCSD for details on generating P1, P2, Q1, and Q2
181 * using ZUNGQR and ZUNGLQ.
182 *
183 * Reference
184 * =========
185 *
186 * [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
187 * Algorithms, 50(1):33-65, 2009.
188 *
189 * ====================================================================
190 *
191 * .. Parameters ..
192 DOUBLE PRECISION REALONE
193 PARAMETER ( REALONE = 1.0D0 )
194 COMPLEX*16 NEGONE, ONE
195 PARAMETER ( NEGONE = (-1.0D0,0.0D0),
196 $ ONE = (1.0D0,0.0D0) )
197 * ..
198 * .. Local Scalars ..
199 LOGICAL COLMAJOR, LQUERY
200 INTEGER I, LWORKMIN, LWORKOPT
201 DOUBLE PRECISION Z1, Z2, Z3, Z4
202 * ..
203 * .. External Subroutines ..
204 EXTERNAL ZAXPY, ZLARF, ZLARFGP, ZSCAL, XERBLA
205 EXTERNAL ZLACGV
206 *
207 * ..
208 * .. External Functions ..
209 DOUBLE PRECISION DZNRM2
210 LOGICAL LSAME
211 EXTERNAL DZNRM2, LSAME
212 * ..
213 * .. Intrinsic Functions
214 INTRINSIC ATAN2, COS, MAX, MIN, SIN
215 INTRINSIC DCMPLX, DCONJG
216 * ..
217 * .. Executable Statements ..
218 *
219 * Test input arguments
220 *
221 INFO = 0
222 COLMAJOR = .NOT. LSAME( TRANS, 'T' )
223 IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN
224 Z1 = REALONE
225 Z2 = REALONE
226 Z3 = REALONE
227 Z4 = REALONE
228 ELSE
229 Z1 = REALONE
230 Z2 = -REALONE
231 Z3 = REALONE
232 Z4 = -REALONE
233 END IF
234 LQUERY = LWORK .EQ. -1
235 *
236 IF( M .LT. 0 ) THEN
237 INFO = -3
238 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
239 INFO = -4
240 ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR.
241 $ Q .GT. M-Q ) THEN
242 INFO = -5
243 ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN
244 INFO = -7
245 ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
246 INFO = -7
247 ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
248 INFO = -9
249 ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
250 INFO = -9
251 ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
252 INFO = -11
253 ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
254 INFO = -11
255 ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
256 INFO = -13
257 ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
258 INFO = -13
259 END IF
260 *
261 * Compute workspace
262 *
263 IF( INFO .EQ. 0 ) THEN
264 LWORKOPT = M - Q
265 LWORKMIN = M - Q
266 WORK(1) = LWORKOPT
267 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
268 INFO = -21
269 END IF
270 END IF
271 IF( INFO .NE. 0 ) THEN
272 CALL XERBLA( 'xORBDB', -INFO )
273 RETURN
274 ELSE IF( LQUERY ) THEN
275 RETURN
276 END IF
277 *
278 * Handle column-major and row-major separately
279 *
280 IF( COLMAJOR ) THEN
281 *
282 * Reduce columns 1, ..., Q of X11, X12, X21, and X22
283 *
284 DO I = 1, Q
285 *
286 IF( I .EQ. 1 ) THEN
287 CALL ZSCAL( P-I+1, DCMPLX( Z1, 0.0D0 ), X11(I,I), 1 )
288 ELSE
289 CALL ZSCAL( P-I+1, DCMPLX( Z1*COS(PHI(I-1)), 0.0D0 ),
290 $ X11(I,I), 1 )
291 CALL ZAXPY( P-I+1, DCMPLX( -Z1*Z3*Z4*SIN(PHI(I-1)),
292 $ 0.0D0 ), X12(I,I-1), 1, X11(I,I), 1 )
293 END IF
294 IF( I .EQ. 1 ) THEN
295 CALL ZSCAL( M-P-I+1, DCMPLX( Z2, 0.0D0 ), X21(I,I), 1 )
296 ELSE
297 CALL ZSCAL( M-P-I+1, DCMPLX( Z2*COS(PHI(I-1)), 0.0D0 ),
298 $ X21(I,I), 1 )
299 CALL ZAXPY( M-P-I+1, DCMPLX( -Z2*Z3*Z4*SIN(PHI(I-1)),
300 $ 0.0D0 ), X22(I,I-1), 1, X21(I,I), 1 )
301 END IF
302 *
303 THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), 1 ),
304 $ DZNRM2( P-I+1, X11(I,I), 1 ) )
305 *
306 CALL ZLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
307 X11(I,I) = ONE
308 CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
309 X21(I,I) = ONE
310 *
311 CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, DCONJG(TAUP1(I)),
312 $ X11(I,I+1), LDX11, WORK )
313 CALL ZLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1,
314 $ DCONJG(TAUP1(I)), X12(I,I), LDX12, WORK )
315 CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1,
316 $ DCONJG(TAUP2(I)), X21(I,I+1), LDX21, WORK )
317 CALL ZLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1,
318 $ DCONJG(TAUP2(I)), X22(I,I), LDX22, WORK )
319 *
320 IF( I .LT. Q ) THEN
321 CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ),
322 $ X11(I,I+1), LDX11 )
323 CALL ZAXPY( Q-I, DCMPLX( Z2*Z3*COS(THETA(I)), 0.0D0 ),
324 $ X21(I,I+1), LDX21, X11(I,I+1), LDX11 )
325 END IF
326 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4*SIN(THETA(I)), 0.0D0 ),
327 $ X12(I,I), LDX12 )
328 CALL ZAXPY( M-Q-I+1, DCMPLX( Z2*Z4*COS(THETA(I)), 0.0D0 ),
329 $ X22(I,I), LDX22, X12(I,I), LDX12 )
330 *
331 IF( I .LT. Q )
332 $ PHI(I) = ATAN2( DZNRM2( Q-I, X11(I,I+1), LDX11 ),
333 $ DZNRM2( M-Q-I+1, X12(I,I), LDX12 ) )
334 *
335 IF( I .LT. Q ) THEN
336 CALL ZLACGV( Q-I, X11(I,I+1), LDX11 )
337 CALL ZLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
338 $ TAUQ1(I) )
339 X11(I,I+1) = ONE
340 END IF
341 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
342 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
343 $ TAUQ2(I) )
344 X12(I,I) = ONE
345 *
346 IF( I .LT. Q ) THEN
347 CALL ZLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
348 $ X11(I+1,I+1), LDX11, WORK )
349 CALL ZLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
350 $ X21(I+1,I+1), LDX21, WORK )
351 END IF
352 CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
353 $ X12(I+1,I), LDX12, WORK )
354 CALL ZLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
355 $ X22(I+1,I), LDX22, WORK )
356 *
357 IF( I .LT. Q )
358 $ CALL ZLACGV( Q-I, X11(I,I+1), LDX11 )
359 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
360 *
361 END DO
362 *
363 * Reduce columns Q + 1, ..., P of X12, X22
364 *
365 DO I = Q + 1, P
366 *
367 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I),
368 $ LDX12 )
369 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
370 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
371 $ TAUQ2(I) )
372 X12(I,I) = ONE
373 *
374 CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
375 $ X12(I+1,I), LDX12, WORK )
376 IF( M-P-Q .GE. 1 )
377 $ CALL ZLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
378 $ TAUQ2(I), X22(Q+1,I), LDX22, WORK )
379 *
380 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
381 *
382 END DO
383 *
384 * Reduce columns P + 1, ..., M - Q of X12, X22
385 *
386 DO I = 1, M - P - Q
387 *
388 CALL ZSCAL( M-P-Q-I+1, DCMPLX( Z2*Z4, 0.0D0 ),
389 $ X22(Q+I,P+I), LDX22 )
390 CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 )
391 CALL ZLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
392 $ LDX22, TAUQ2(P+I) )
393 X22(Q+I,P+I) = ONE
394 CALL ZLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
395 $ TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
396 *
397 CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 )
398 *
399 END DO
400 *
401 ELSE
402 *
403 * Reduce columns 1, ..., Q of X11, X12, X21, X22
404 *
405 DO I = 1, Q
406 *
407 IF( I .EQ. 1 ) THEN
408 CALL ZSCAL( P-I+1, DCMPLX( Z1, 0.0D0 ), X11(I,I),
409 $ LDX11 )
410 ELSE
411 CALL ZSCAL( P-I+1, DCMPLX( Z1*COS(PHI(I-1)), 0.0D0 ),
412 $ X11(I,I), LDX11 )
413 CALL ZAXPY( P-I+1, DCMPLX( -Z1*Z3*Z4*SIN(PHI(I-1)),
414 $ 0.0D0 ), X12(I-1,I), LDX12, X11(I,I), LDX11 )
415 END IF
416 IF( I .EQ. 1 ) THEN
417 CALL ZSCAL( M-P-I+1, DCMPLX( Z2, 0.0D0 ), X21(I,I),
418 $ LDX21 )
419 ELSE
420 CALL ZSCAL( M-P-I+1, DCMPLX( Z2*COS(PHI(I-1)), 0.0D0 ),
421 $ X21(I,I), LDX21 )
422 CALL ZAXPY( M-P-I+1, DCMPLX( -Z2*Z3*Z4*SIN(PHI(I-1)),
423 $ 0.0D0 ), X22(I-1,I), LDX22, X21(I,I), LDX21 )
424 END IF
425 *
426 THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), LDX21 ),
427 $ DZNRM2( P-I+1, X11(I,I), LDX11 ) )
428 *
429 CALL ZLACGV( P-I+1, X11(I,I), LDX11 )
430 CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 )
431 *
432 CALL ZLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
433 X11(I,I) = ONE
434 CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
435 $ TAUP2(I) )
436 X21(I,I) = ONE
437 *
438 CALL ZLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
439 $ X11(I+1,I), LDX11, WORK )
440 CALL ZLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I),
441 $ X12(I,I), LDX12, WORK )
442 CALL ZLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
443 $ X21(I+1,I), LDX21, WORK )
444 CALL ZLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
445 $ TAUP2(I), X22(I,I), LDX22, WORK )
446 *
447 CALL ZLACGV( P-I+1, X11(I,I), LDX11 )
448 CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 )
449 *
450 IF( I .LT. Q ) THEN
451 CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ),
452 $ X11(I+1,I), 1 )
453 CALL ZAXPY( Q-I, DCMPLX( Z2*Z3*COS(THETA(I)), 0.0D0 ),
454 $ X21(I+1,I), 1, X11(I+1,I), 1 )
455 END IF
456 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4*SIN(THETA(I)), 0.0D0 ),
457 $ X12(I,I), 1 )
458 CALL ZAXPY( M-Q-I+1, DCMPLX( Z2*Z4*COS(THETA(I)), 0.0D0 ),
459 $ X22(I,I), 1, X12(I,I), 1 )
460 *
461 IF( I .LT. Q )
462 $ PHI(I) = ATAN2( DZNRM2( Q-I, X11(I+1,I), 1 ),
463 $ DZNRM2( M-Q-I+1, X12(I,I), 1 ) )
464 *
465 IF( I .LT. Q ) THEN
466 CALL ZLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) )
467 X11(I+1,I) = ONE
468 END IF
469 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
470 X12(I,I) = ONE
471 *
472 IF( I .LT. Q ) THEN
473 CALL ZLARF( 'L', Q-I, P-I, X11(I+1,I), 1,
474 $ DCONJG(TAUQ1(I)), X11(I+1,I+1), LDX11, WORK )
475 CALL ZLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1,
476 $ DCONJG(TAUQ1(I)), X21(I+1,I+1), LDX21, WORK )
477 END IF
478 CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1,
479 $ DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK )
480 CALL ZLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1,
481 $ DCONJG(TAUQ2(I)), X22(I,I+1), LDX22, WORK )
482 *
483 END DO
484 *
485 * Reduce columns Q + 1, ..., P of X12, X22
486 *
487 DO I = Q + 1, P
488 *
489 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I), 1 )
490 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
491 X12(I,I) = ONE
492 *
493 CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1,
494 $ DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK )
495 IF( M-P-Q .GE. 1 )
496 $ CALL ZLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1,
497 $ DCONJG(TAUQ2(I)), X22(I,Q+1), LDX22, WORK )
498 *
499 END DO
500 *
501 * Reduce columns P + 1, ..., M - Q of X12, X22
502 *
503 DO I = 1, M - P - Q
504 *
505 CALL ZSCAL( M-P-Q-I+1, DCMPLX( Z2*Z4, 0.0D0 ),
506 $ X22(P+I,Q+I), 1 )
507 CALL ZLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
508 $ TAUQ2(P+I) )
509 X22(P+I,Q+I) = ONE
510 *
511 CALL ZLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
512 $ DCONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), LDX22,
513 $ WORK )
514 *
515 END DO
516 *
517 END IF
518 *
519 RETURN
520 *
521 * End of ZUNBDB
522 *
523 END
524