LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zunbdb.f
Go to the documentation of this file.
1 *> \brief \b ZUNBDB
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZUNBDB + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zunbdb.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zunbdb.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zunbdb.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
22 * X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
23 * TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER SIGNS, TRANS
27 * INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
28 * $ Q
29 * ..
30 * .. Array Arguments ..
31 * DOUBLE PRECISION PHI( * ), THETA( * )
32 * COMPLEX*16 TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
33 * $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
34 * $ X21( LDX21, * ), X22( LDX22, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> ZUNBDB simultaneously bidiagonalizes the blocks of an M-by-M
44 *> partitioned unitary matrix X:
45 *>
46 *> [ B11 | B12 0 0 ]
47 *> [ X11 | X12 ] [ P1 | ] [ 0 | 0 -I 0 ] [ Q1 | ]**H
48 *> X = [-----------] = [---------] [----------------] [---------] .
49 *> [ X21 | X22 ] [ | P2 ] [ B21 | B22 0 0 ] [ | Q2 ]
50 *> [ 0 | 0 0 I ]
51 *>
52 *> X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
53 *> not the case, then X must be transposed and/or permuted. This can be
54 *> done in constant time using the TRANS and SIGNS options. See ZUNCSD
55 *> for details.)
56 *>
57 *> The unitary matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
58 *> (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
59 *> represented implicitly by Householder vectors.
60 *>
61 *> B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
62 *> implicitly by angles THETA, PHI.
63 *> \endverbatim
64 *
65 * Arguments:
66 * ==========
67 *
68 *> \param[in] TRANS
69 *> \verbatim
70 *> TRANS is CHARACTER
71 *> = 'T': X, U1, U2, V1T, and V2T are stored in row-major
72 *> order;
73 *> otherwise: X, U1, U2, V1T, and V2T are stored in column-
74 *> major order.
75 *> \endverbatim
76 *>
77 *> \param[in] SIGNS
78 *> \verbatim
79 *> SIGNS is CHARACTER
80 *> = 'O': The lower-left block is made nonpositive (the
81 *> "other" convention);
82 *> otherwise: The upper-right block is made nonpositive (the
83 *> "default" convention).
84 *> \endverbatim
85 *>
86 *> \param[in] M
87 *> \verbatim
88 *> M is INTEGER
89 *> The number of rows and columns in X.
90 *> \endverbatim
91 *>
92 *> \param[in] P
93 *> \verbatim
94 *> P is INTEGER
95 *> The number of rows in X11 and X12. 0 <= P <= M.
96 *> \endverbatim
97 *>
98 *> \param[in] Q
99 *> \verbatim
100 *> Q is INTEGER
101 *> The number of columns in X11 and X21. 0 <= Q <=
102 *> MIN(P,M-P,M-Q).
103 *> \endverbatim
104 *>
105 *> \param[in,out] X11
106 *> \verbatim
107 *> X11 is COMPLEX*16 array, dimension (LDX11,Q)
108 *> On entry, the top-left block of the unitary matrix to be
109 *> reduced. On exit, the form depends on TRANS:
110 *> If TRANS = 'N', then
111 *> the columns of tril(X11) specify reflectors for P1,
112 *> the rows of triu(X11,1) specify reflectors for Q1;
113 *> else TRANS = 'T', and
114 *> the rows of triu(X11) specify reflectors for P1,
115 *> the columns of tril(X11,-1) specify reflectors for Q1.
116 *> \endverbatim
117 *>
118 *> \param[in] LDX11
119 *> \verbatim
120 *> LDX11 is INTEGER
121 *> The leading dimension of X11. If TRANS = 'N', then LDX11 >=
122 *> P; else LDX11 >= Q.
123 *> \endverbatim
124 *>
125 *> \param[in,out] X12
126 *> \verbatim
127 *> X12 is COMPLEX*16 array, dimension (LDX12,M-Q)
128 *> On entry, the top-right block of the unitary matrix to
129 *> be reduced. On exit, the form depends on TRANS:
130 *> If TRANS = 'N', then
131 *> the rows of triu(X12) specify the first P reflectors for
132 *> Q2;
133 *> else TRANS = 'T', and
134 *> the columns of tril(X12) specify the first P reflectors
135 *> for Q2.
136 *> \endverbatim
137 *>
138 *> \param[in] LDX12
139 *> \verbatim
140 *> LDX12 is INTEGER
141 *> The leading dimension of X12. If TRANS = 'N', then LDX12 >=
142 *> P; else LDX11 >= M-Q.
143 *> \endverbatim
144 *>
145 *> \param[in,out] X21
146 *> \verbatim
147 *> X21 is COMPLEX*16 array, dimension (LDX21,Q)
148 *> On entry, the bottom-left block of the unitary matrix to
149 *> be reduced. On exit, the form depends on TRANS:
150 *> If TRANS = 'N', then
151 *> the columns of tril(X21) specify reflectors for P2;
152 *> else TRANS = 'T', and
153 *> the rows of triu(X21) specify reflectors for P2.
154 *> \endverbatim
155 *>
156 *> \param[in] LDX21
157 *> \verbatim
158 *> LDX21 is INTEGER
159 *> The leading dimension of X21. If TRANS = 'N', then LDX21 >=
160 *> M-P; else LDX21 >= Q.
161 *> \endverbatim
162 *>
163 *> \param[in,out] X22
164 *> \verbatim
165 *> X22 is COMPLEX*16 array, dimension (LDX22,M-Q)
166 *> On entry, the bottom-right block of the unitary matrix to
167 *> be reduced. On exit, the form depends on TRANS:
168 *> If TRANS = 'N', then
169 *> the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
170 *> M-P-Q reflectors for Q2,
171 *> else TRANS = 'T', and
172 *> the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
173 *> M-P-Q reflectors for P2.
174 *> \endverbatim
175 *>
176 *> \param[in] LDX22
177 *> \verbatim
178 *> LDX22 is INTEGER
179 *> The leading dimension of X22. If TRANS = 'N', then LDX22 >=
180 *> M-P; else LDX22 >= M-Q.
181 *> \endverbatim
182 *>
183 *> \param[out] THETA
184 *> \verbatim
185 *> THETA is DOUBLE PRECISION array, dimension (Q)
186 *> The entries of the bidiagonal blocks B11, B12, B21, B22 can
187 *> be computed from the angles THETA and PHI. See Further
188 *> Details.
189 *> \endverbatim
190 *>
191 *> \param[out] PHI
192 *> \verbatim
193 *> PHI is DOUBLE PRECISION array, dimension (Q-1)
194 *> The entries of the bidiagonal blocks B11, B12, B21, B22 can
195 *> be computed from the angles THETA and PHI. See Further
196 *> Details.
197 *> \endverbatim
198 *>
199 *> \param[out] TAUP1
200 *> \verbatim
201 *> TAUP1 is COMPLEX*16 array, dimension (P)
202 *> The scalar factors of the elementary reflectors that define
203 *> P1.
204 *> \endverbatim
205 *>
206 *> \param[out] TAUP2
207 *> \verbatim
208 *> TAUP2 is COMPLEX*16 array, dimension (M-P)
209 *> The scalar factors of the elementary reflectors that define
210 *> P2.
211 *> \endverbatim
212 *>
213 *> \param[out] TAUQ1
214 *> \verbatim
215 *> TAUQ1 is COMPLEX*16 array, dimension (Q)
216 *> The scalar factors of the elementary reflectors that define
217 *> Q1.
218 *> \endverbatim
219 *>
220 *> \param[out] TAUQ2
221 *> \verbatim
222 *> TAUQ2 is COMPLEX*16 array, dimension (M-Q)
223 *> The scalar factors of the elementary reflectors that define
224 *> Q2.
225 *> \endverbatim
226 *>
227 *> \param[out] WORK
228 *> \verbatim
229 *> WORK is COMPLEX*16 array, dimension (LWORK)
230 *> \endverbatim
231 *>
232 *> \param[in] LWORK
233 *> \verbatim
234 *> LWORK is INTEGER
235 *> The dimension of the array WORK. LWORK >= M-Q.
236 *>
237 *> If LWORK = -1, then a workspace query is assumed; the routine
238 *> only calculates the optimal size of the WORK array, returns
239 *> this value as the first entry of the WORK array, and no error
240 *> message related to LWORK is issued by XERBLA.
241 *> \endverbatim
242 *>
243 *> \param[out] INFO
244 *> \verbatim
245 *> INFO is INTEGER
246 *> = 0: successful exit.
247 *> < 0: if INFO = -i, the i-th argument had an illegal value.
248 *> \endverbatim
249 *
250 * Authors:
251 * ========
252 *
253 *> \author Univ. of Tennessee
254 *> \author Univ. of California Berkeley
255 *> \author Univ. of Colorado Denver
256 *> \author NAG Ltd.
257 *
258 *> \ingroup complex16OTHERcomputational
259 *
260 *> \par Further Details:
261 * =====================
262 *>
263 *> \verbatim
264 *>
265 *> The bidiagonal blocks B11, B12, B21, and B22 are represented
266 *> implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
267 *> PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
268 *> lower bidiagonal. Every entry in each bidiagonal band is a product
269 *> of a sine or cosine of a THETA with a sine or cosine of a PHI. See
270 *> [1] or ZUNCSD for details.
271 *>
272 *> P1, P2, Q1, and Q2 are represented as products of elementary
273 *> reflectors. See ZUNCSD for details on generating P1, P2, Q1, and Q2
274 *> using ZUNGQR and ZUNGLQ.
275 *> \endverbatim
276 *
277 *> \par References:
278 * ================
279 *>
280 *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
281 *> Algorithms, 50(1):33-65, 2009.
282 *>
283 * =====================================================================
284  SUBROUTINE zunbdb( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
285  $ X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
286  $ TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
287 *
288 * -- LAPACK computational routine --
289 * -- LAPACK is a software package provided by Univ. of Tennessee, --
290 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
291 *
292 * .. Scalar Arguments ..
293  CHARACTER SIGNS, TRANS
294  INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
295  $ q
296 * ..
297 * .. Array Arguments ..
298  DOUBLE PRECISION PHI( * ), THETA( * )
299  COMPLEX*16 TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
300  $ work( * ), x11( ldx11, * ), x12( ldx12, * ),
301  $ x21( ldx21, * ), x22( ldx22, * )
302 * ..
303 *
304 * ====================================================================
305 *
306 * .. Parameters ..
307  DOUBLE PRECISION REALONE
308  PARAMETER ( REALONE = 1.0d0 )
309  COMPLEX*16 ONE
310  parameter( one = (1.0d0,0.0d0) )
311 * ..
312 * .. Local Scalars ..
313  LOGICAL COLMAJOR, LQUERY
314  INTEGER I, LWORKMIN, LWORKOPT
315  DOUBLE PRECISION Z1, Z2, Z3, Z4
316 * ..
317 * .. External Subroutines ..
318  EXTERNAL zaxpy, zlarf, zlarfgp, zscal, xerbla
319  EXTERNAL zlacgv
320 *
321 * ..
322 * .. External Functions ..
323  DOUBLE PRECISION DZNRM2
324  LOGICAL LSAME
325  EXTERNAL dznrm2, lsame
326 * ..
327 * .. Intrinsic Functions
328  INTRINSIC atan2, cos, max, min, sin
329  INTRINSIC dcmplx, dconjg
330 * ..
331 * .. Executable Statements ..
332 *
333 * Test input arguments
334 *
335  info = 0
336  colmajor = .NOT. lsame( trans, 'T' )
337  IF( .NOT. lsame( signs, 'O' ) ) THEN
338  z1 = realone
339  z2 = realone
340  z3 = realone
341  z4 = realone
342  ELSE
343  z1 = realone
344  z2 = -realone
345  z3 = realone
346  z4 = -realone
347  END IF
348  lquery = lwork .EQ. -1
349 *
350  IF( m .LT. 0 ) THEN
351  info = -3
352  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
353  info = -4
354  ELSE IF( q .LT. 0 .OR. q .GT. p .OR. q .GT. m-p .OR.
355  $ q .GT. m-q ) THEN
356  info = -5
357  ELSE IF( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
358  info = -7
359  ELSE IF( .NOT.colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
360  info = -7
361  ELSE IF( colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
362  info = -9
363  ELSE IF( .NOT.colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
364  info = -9
365  ELSE IF( colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
366  info = -11
367  ELSE IF( .NOT.colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
368  info = -11
369  ELSE IF( colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
370  info = -13
371  ELSE IF( .NOT.colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
372  info = -13
373  END IF
374 *
375 * Compute workspace
376 *
377  IF( info .EQ. 0 ) THEN
378  lworkopt = m - q
379  lworkmin = m - q
380  work(1) = lworkopt
381  IF( lwork .LT. lworkmin .AND. .NOT. lquery ) THEN
382  info = -21
383  END IF
384  END IF
385  IF( info .NE. 0 ) THEN
386  CALL xerbla( 'xORBDB', -info )
387  RETURN
388  ELSE IF( lquery ) THEN
389  RETURN
390  END IF
391 *
392 * Handle column-major and row-major separately
393 *
394  IF( colmajor ) THEN
395 *
396 * Reduce columns 1, ..., Q of X11, X12, X21, and X22
397 *
398  DO i = 1, q
399 *
400  IF( i .EQ. 1 ) THEN
401  CALL zscal( p-i+1, dcmplx( z1, 0.0d0 ), x11(i,i), 1 )
402  ELSE
403  CALL zscal( p-i+1, dcmplx( z1*cos(phi(i-1)), 0.0d0 ),
404  $ x11(i,i), 1 )
405  CALL zaxpy( p-i+1, dcmplx( -z1*z3*z4*sin(phi(i-1)),
406  $ 0.0d0 ), x12(i,i-1), 1, x11(i,i), 1 )
407  END IF
408  IF( i .EQ. 1 ) THEN
409  CALL zscal( m-p-i+1, dcmplx( z2, 0.0d0 ), x21(i,i), 1 )
410  ELSE
411  CALL zscal( m-p-i+1, dcmplx( z2*cos(phi(i-1)), 0.0d0 ),
412  $ x21(i,i), 1 )
413  CALL zaxpy( m-p-i+1, dcmplx( -z2*z3*z4*sin(phi(i-1)),
414  $ 0.0d0 ), x22(i,i-1), 1, x21(i,i), 1 )
415  END IF
416 *
417  theta(i) = atan2( dznrm2( m-p-i+1, x21(i,i), 1 ),
418  $ dznrm2( p-i+1, x11(i,i), 1 ) )
419 *
420  IF( p .GT. i ) THEN
421  CALL zlarfgp( p-i+1, x11(i,i), x11(i+1,i), 1, taup1(i) )
422  ELSE IF ( p .EQ. i ) THEN
423  CALL zlarfgp( p-i+1, x11(i,i), x11(i,i), 1, taup1(i) )
424  END IF
425  x11(i,i) = one
426  IF ( m-p .GT. i ) THEN
427  CALL zlarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1,
428  $ taup2(i) )
429  ELSE IF ( m-p .EQ. i ) THEN
430  CALL zlarfgp( m-p-i+1, x21(i,i), x21(i,i), 1,
431  $ taup2(i) )
432  END IF
433  x21(i,i) = one
434 *
435  IF ( q .GT. i ) THEN
436  CALL zlarf( 'L', p-i+1, q-i, x11(i,i), 1,
437  $ dconjg(taup1(i)), x11(i,i+1), ldx11, work )
438  CALL zlarf( 'L', m-p-i+1, q-i, x21(i,i), 1,
439  $ dconjg(taup2(i)), x21(i,i+1), ldx21, work )
440  END IF
441  IF ( m-q+1 .GT. i ) THEN
442  CALL zlarf( 'L', p-i+1, m-q-i+1, x11(i,i), 1,
443  $ dconjg(taup1(i)), x12(i,i), ldx12, work )
444  CALL zlarf( 'L', m-p-i+1, m-q-i+1, x21(i,i), 1,
445  $ dconjg(taup2(i)), x22(i,i), ldx22, work )
446  END IF
447 *
448  IF( i .LT. q ) THEN
449  CALL zscal( q-i, dcmplx( -z1*z3*sin(theta(i)), 0.0d0 ),
450  $ x11(i,i+1), ldx11 )
451  CALL zaxpy( q-i, dcmplx( z2*z3*cos(theta(i)), 0.0d0 ),
452  $ x21(i,i+1), ldx21, x11(i,i+1), ldx11 )
453  END IF
454  CALL zscal( m-q-i+1, dcmplx( -z1*z4*sin(theta(i)), 0.0d0 ),
455  $ x12(i,i), ldx12 )
456  CALL zaxpy( m-q-i+1, dcmplx( z2*z4*cos(theta(i)), 0.0d0 ),
457  $ x22(i,i), ldx22, x12(i,i), ldx12 )
458 *
459  IF( i .LT. q )
460  $ phi(i) = atan2( dznrm2( q-i, x11(i,i+1), ldx11 ),
461  $ dznrm2( m-q-i+1, x12(i,i), ldx12 ) )
462 *
463  IF( i .LT. q ) THEN
464  CALL zlacgv( q-i, x11(i,i+1), ldx11 )
465  IF ( i .EQ. q-1 ) THEN
466  CALL zlarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
467  $ tauq1(i) )
468  ELSE
469  CALL zlarfgp( q-i, x11(i,i+1), x11(i,i+2), ldx11,
470  $ tauq1(i) )
471  END IF
472  x11(i,i+1) = one
473  END IF
474  IF ( m-q+1 .GT. i ) THEN
475  CALL zlacgv( m-q-i+1, x12(i,i), ldx12 )
476  IF ( m-q .EQ. i ) THEN
477  CALL zlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
478  $ tauq2(i) )
479  ELSE
480  CALL zlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
481  $ tauq2(i) )
482  END IF
483  END IF
484  x12(i,i) = one
485 *
486  IF( i .LT. q ) THEN
487  CALL zlarf( 'R', p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
488  $ x11(i+1,i+1), ldx11, work )
489  CALL zlarf( 'R', m-p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
490  $ x21(i+1,i+1), ldx21, work )
491  END IF
492  IF ( p .GT. i ) THEN
493  CALL zlarf( 'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
494  $ x12(i+1,i), ldx12, work )
495  END IF
496  IF ( m-p .GT. i ) THEN
497  CALL zlarf( 'R', m-p-i, m-q-i+1, x12(i,i), ldx12,
498  $ tauq2(i), x22(i+1,i), ldx22, work )
499  END IF
500 *
501  IF( i .LT. q )
502  $ CALL zlacgv( q-i, x11(i,i+1), ldx11 )
503  CALL zlacgv( m-q-i+1, x12(i,i), ldx12 )
504 *
505  END DO
506 *
507 * Reduce columns Q + 1, ..., P of X12, X22
508 *
509  DO i = q + 1, p
510 *
511  CALL zscal( m-q-i+1, dcmplx( -z1*z4, 0.0d0 ), x12(i,i),
512  $ ldx12 )
513  CALL zlacgv( m-q-i+1, x12(i,i), ldx12 )
514  IF ( i .GE. m-q ) THEN
515  CALL zlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
516  $ tauq2(i) )
517  ELSE
518  CALL zlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
519  $ tauq2(i) )
520  END IF
521  x12(i,i) = one
522 *
523  IF ( p .GT. i ) THEN
524  CALL zlarf( 'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
525  $ x12(i+1,i), ldx12, work )
526  END IF
527  IF( m-p-q .GE. 1 )
528  $ CALL zlarf( 'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
529  $ tauq2(i), x22(q+1,i), ldx22, work )
530 *
531  CALL zlacgv( m-q-i+1, x12(i,i), ldx12 )
532 *
533  END DO
534 *
535 * Reduce columns P + 1, ..., M - Q of X12, X22
536 *
537  DO i = 1, m - p - q
538 *
539  CALL zscal( m-p-q-i+1, dcmplx( z2*z4, 0.0d0 ),
540  $ x22(q+i,p+i), ldx22 )
541  CALL zlacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
542  CALL zlarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i+1),
543  $ ldx22, tauq2(p+i) )
544  x22(q+i,p+i) = one
545  CALL zlarf( 'R', m-p-q-i, m-p-q-i+1, x22(q+i,p+i), ldx22,
546  $ tauq2(p+i), x22(q+i+1,p+i), ldx22, work )
547 *
548  CALL zlacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
549 *
550  END DO
551 *
552  ELSE
553 *
554 * Reduce columns 1, ..., Q of X11, X12, X21, X22
555 *
556  DO i = 1, q
557 *
558  IF( i .EQ. 1 ) THEN
559  CALL zscal( p-i+1, dcmplx( z1, 0.0d0 ), x11(i,i),
560  $ ldx11 )
561  ELSE
562  CALL zscal( p-i+1, dcmplx( z1*cos(phi(i-1)), 0.0d0 ),
563  $ x11(i,i), ldx11 )
564  CALL zaxpy( p-i+1, dcmplx( -z1*z3*z4*sin(phi(i-1)),
565  $ 0.0d0 ), x12(i-1,i), ldx12, x11(i,i), ldx11 )
566  END IF
567  IF( i .EQ. 1 ) THEN
568  CALL zscal( m-p-i+1, dcmplx( z2, 0.0d0 ), x21(i,i),
569  $ ldx21 )
570  ELSE
571  CALL zscal( m-p-i+1, dcmplx( z2*cos(phi(i-1)), 0.0d0 ),
572  $ x21(i,i), ldx21 )
573  CALL zaxpy( m-p-i+1, dcmplx( -z2*z3*z4*sin(phi(i-1)),
574  $ 0.0d0 ), x22(i-1,i), ldx22, x21(i,i), ldx21 )
575  END IF
576 *
577  theta(i) = atan2( dznrm2( m-p-i+1, x21(i,i), ldx21 ),
578  $ dznrm2( p-i+1, x11(i,i), ldx11 ) )
579 *
580  CALL zlacgv( p-i+1, x11(i,i), ldx11 )
581  CALL zlacgv( m-p-i+1, x21(i,i), ldx21 )
582 *
583  CALL zlarfgp( p-i+1, x11(i,i), x11(i,i+1), ldx11, taup1(i) )
584  x11(i,i) = one
585  IF ( i .EQ. m-p ) THEN
586  CALL zlarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
587  $ taup2(i) )
588  ELSE
589  CALL zlarfgp( m-p-i+1, x21(i,i), x21(i,i+1), ldx21,
590  $ taup2(i) )
591  END IF
592  x21(i,i) = one
593 *
594  CALL zlarf( 'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
595  $ x11(i+1,i), ldx11, work )
596  CALL zlarf( 'R', m-q-i+1, p-i+1, x11(i,i), ldx11, taup1(i),
597  $ x12(i,i), ldx12, work )
598  CALL zlarf( 'R', q-i, m-p-i+1, x21(i,i), ldx21, taup2(i),
599  $ x21(i+1,i), ldx21, work )
600  CALL zlarf( 'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
601  $ taup2(i), x22(i,i), ldx22, work )
602 *
603  CALL zlacgv( p-i+1, x11(i,i), ldx11 )
604  CALL zlacgv( m-p-i+1, x21(i,i), ldx21 )
605 *
606  IF( i .LT. q ) THEN
607  CALL zscal( q-i, dcmplx( -z1*z3*sin(theta(i)), 0.0d0 ),
608  $ x11(i+1,i), 1 )
609  CALL zaxpy( q-i, dcmplx( z2*z3*cos(theta(i)), 0.0d0 ),
610  $ x21(i+1,i), 1, x11(i+1,i), 1 )
611  END IF
612  CALL zscal( m-q-i+1, dcmplx( -z1*z4*sin(theta(i)), 0.0d0 ),
613  $ x12(i,i), 1 )
614  CALL zaxpy( m-q-i+1, dcmplx( z2*z4*cos(theta(i)), 0.0d0 ),
615  $ x22(i,i), 1, x12(i,i), 1 )
616 *
617  IF( i .LT. q )
618  $ phi(i) = atan2( dznrm2( q-i, x11(i+1,i), 1 ),
619  $ dznrm2( m-q-i+1, x12(i,i), 1 ) )
620 *
621  IF( i .LT. q ) THEN
622  CALL zlarfgp( q-i, x11(i+1,i), x11(i+2,i), 1, tauq1(i) )
623  x11(i+1,i) = one
624  END IF
625  CALL zlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1, tauq2(i) )
626  x12(i,i) = one
627 *
628  IF( i .LT. q ) THEN
629  CALL zlarf( 'L', q-i, p-i, x11(i+1,i), 1,
630  $ dconjg(tauq1(i)), x11(i+1,i+1), ldx11, work )
631  CALL zlarf( 'L', q-i, m-p-i, x11(i+1,i), 1,
632  $ dconjg(tauq1(i)), x21(i+1,i+1), ldx21, work )
633  END IF
634  CALL zlarf( 'L', m-q-i+1, p-i, x12(i,i), 1,
635  $ dconjg(tauq2(i)), x12(i,i+1), ldx12, work )
636  IF ( m-p .GT. i ) THEN
637  CALL zlarf( 'L', m-q-i+1, m-p-i, x12(i,i), 1,
638  $ dconjg(tauq2(i)), x22(i,i+1), ldx22, work )
639  END IF
640 *
641  END DO
642 *
643 * Reduce columns Q + 1, ..., P of X12, X22
644 *
645  DO i = q + 1, p
646 *
647  CALL zscal( m-q-i+1, dcmplx( -z1*z4, 0.0d0 ), x12(i,i), 1 )
648  CALL zlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1, tauq2(i) )
649  x12(i,i) = one
650 *
651  IF ( p .GT. i ) THEN
652  CALL zlarf( 'L', m-q-i+1, p-i, x12(i,i), 1,
653  $ dconjg(tauq2(i)), x12(i,i+1), ldx12, work )
654  END IF
655  IF( m-p-q .GE. 1 )
656  $ CALL zlarf( 'L', m-q-i+1, m-p-q, x12(i,i), 1,
657  $ dconjg(tauq2(i)), x22(i,q+1), ldx22, work )
658 *
659  END DO
660 *
661 * Reduce columns P + 1, ..., M - Q of X12, X22
662 *
663  DO i = 1, m - p - q
664 *
665  CALL zscal( m-p-q-i+1, dcmplx( z2*z4, 0.0d0 ),
666  $ x22(p+i,q+i), 1 )
667  CALL zlarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i+1,q+i), 1,
668  $ tauq2(p+i) )
669  x22(p+i,q+i) = one
670 *
671  IF ( m-p-q .NE. i ) THEN
672  CALL zlarf( 'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i), 1,
673  $ dconjg(tauq2(p+i)), x22(p+i,q+i+1), ldx22,
674  $ work )
675  END IF
676 *
677  END DO
678 *
679  END IF
680 *
681  RETURN
682 *
683 * End of ZUNBDB
684 *
685  END
686 
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zlarfgp(N, ALPHA, X, INCX, TAU)
ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: zlarfgp.f:104
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:74
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:128
subroutine zunbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
ZUNBDB
Definition: zunbdb.f:287