LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cunbdb.f
Go to the documentation of this file.
1*> \brief \b CUNBDB
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CUNBDB + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cunbdb.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cunbdb.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cunbdb.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CUNBDB( 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* REAL PHI( * ), THETA( * )
32* COMPLEX 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*> CUNBDB 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 CUNCSD
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 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 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 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 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 REAL 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 REAL 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 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 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 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 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 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 unbdb
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 CUNCSD for details.
271*>
272*> P1, P2, Q1, and Q2 are represented as products of elementary
273*> reflectors. See CUNCSD for details on generating P1, P2, Q1, and Q2
274*> using CUNGQR and CUNGLQ.
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 cunbdb( 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 REAL PHI( * ), THETA( * )
299 COMPLEX TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
300 $ work( * ), x11( ldx11, * ), x12( ldx12, * ),
301 $ x21( ldx21, * ), x22( ldx22, * )
302* ..
303*
304* ====================================================================
305*
306* .. Parameters ..
307 REAL REALONE
308 PARAMETER ( REALONE = 1.0e0 )
309 COMPLEX ONE
310 parameter( one = (1.0e0,0.0e0) )
311* ..
312* .. Local Scalars ..
313 LOGICAL COLMAJOR, LQUERY
314 INTEGER I, LWORKMIN, LWORKOPT
315 REAL Z1, Z2, Z3, Z4
316* ..
317* .. External Subroutines ..
318 EXTERNAL caxpy, clarf, clarfgp, cscal, xerbla
319 EXTERNAL clacgv
320*
321* ..
322* .. External Functions ..
323 REAL SCNRM2, SROUNDUP_LWORK
324 LOGICAL LSAME
325 EXTERNAL scnrm2, sroundup_lwork, lsame
326* ..
327* .. Intrinsic Functions
328 INTRINSIC atan2, cos, max, min, sin
329 INTRINSIC cmplx, conjg
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) = sroundup_lwork(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 cscal( p-i+1, cmplx( z1, 0.0e0 ), x11(i,i), 1 )
402 ELSE
403 CALL cscal( p-i+1, cmplx( z1*cos(phi(i-1)), 0.0e0 ),
404 $ x11(i,i), 1 )
405 CALL caxpy( p-i+1, cmplx( -z1*z3*z4*sin(phi(i-1)),
406 $ 0.0e0 ), x12(i,i-1), 1, x11(i,i), 1 )
407 END IF
408 IF( i .EQ. 1 ) THEN
409 CALL cscal( m-p-i+1, cmplx( z2, 0.0e0 ), x21(i,i), 1 )
410 ELSE
411 CALL cscal( m-p-i+1, cmplx( z2*cos(phi(i-1)), 0.0e0 ),
412 $ x21(i,i), 1 )
413 CALL caxpy( m-p-i+1, cmplx( -z2*z3*z4*sin(phi(i-1)),
414 $ 0.0e0 ), x22(i,i-1), 1, x21(i,i), 1 )
415 END IF
416*
417 theta(i) = atan2( scnrm2( m-p-i+1, x21(i,i), 1 ),
418 $ scnrm2( p-i+1, x11(i,i), 1 ) )
419*
420 IF( p .GT. i ) THEN
421 CALL clarfgp( p-i+1, x11(i,i), x11(i+1,i), 1, taup1(i) )
422 ELSE IF ( p .EQ. i ) THEN
423 CALL clarfgp( 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 clarfgp( 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 clarfgp( 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 clarf( 'L', p-i+1, q-i, x11(i,i), 1,
437 $ conjg(taup1(i)), x11(i,i+1), ldx11, work )
438 CALL clarf( 'L', m-p-i+1, q-i, x21(i,i), 1,
439 $ conjg(taup2(i)), x21(i,i+1), ldx21, work )
440 END IF
441 IF ( m-q+1 .GT. i ) THEN
442 CALL clarf( 'L', p-i+1, m-q-i+1, x11(i,i), 1,
443 $ conjg(taup1(i)), x12(i,i), ldx12, work )
444 CALL clarf( 'L', m-p-i+1, m-q-i+1, x21(i,i), 1,
445 $ conjg(taup2(i)), x22(i,i), ldx22, work )
446 END IF
447*
448 IF( i .LT. q ) THEN
449 CALL cscal( q-i, cmplx( -z1*z3*sin(theta(i)), 0.0e0 ),
450 $ x11(i,i+1), ldx11 )
451 CALL caxpy( q-i, cmplx( z2*z3*cos(theta(i)), 0.0e0 ),
452 $ x21(i,i+1), ldx21, x11(i,i+1), ldx11 )
453 END IF
454 CALL cscal( m-q-i+1, cmplx( -z1*z4*sin(theta(i)), 0.0e0 ),
455 $ x12(i,i), ldx12 )
456 CALL caxpy( m-q-i+1, cmplx( z2*z4*cos(theta(i)), 0.0e0 ),
457 $ x22(i,i), ldx22, x12(i,i), ldx12 )
458*
459 IF( i .LT. q )
460 $ phi(i) = atan2( scnrm2( q-i, x11(i,i+1), ldx11 ),
461 $ scnrm2( m-q-i+1, x12(i,i), ldx12 ) )
462*
463 IF( i .LT. q ) THEN
464 CALL clacgv( q-i, x11(i,i+1), ldx11 )
465 IF ( i .EQ. q-1 ) THEN
466 CALL clarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
467 $ tauq1(i) )
468 ELSE
469 CALL clarfgp( 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 clacgv( m-q-i+1, x12(i,i), ldx12 )
476 IF ( m-q .EQ. i ) THEN
477 CALL clarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
478 $ tauq2(i) )
479 ELSE
480 CALL clarfgp( 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 clarf( 'R', p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
488 $ x11(i+1,i+1), ldx11, work )
489 CALL clarf( '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 clarf( '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 clarf( '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 clacgv( q-i, x11(i,i+1), ldx11 )
503 CALL clacgv( 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 cscal( m-q-i+1, cmplx( -z1*z4, 0.0e0 ), x12(i,i),
512 $ ldx12 )
513 CALL clacgv( m-q-i+1, x12(i,i), ldx12 )
514 IF ( i .GE. m-q ) THEN
515 CALL clarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
516 $ tauq2(i) )
517 ELSE
518 CALL clarfgp( 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 clarf( '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 clarf( 'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
529 $ tauq2(i), x22(q+1,i), ldx22, work )
530*
531 CALL clacgv( 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 cscal( m-p-q-i+1, cmplx( z2*z4, 0.0e0 ),
540 $ x22(q+i,p+i), ldx22 )
541 CALL clacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
542 CALL clarfgp( 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 clarf( '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 clacgv( 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 cscal( p-i+1, cmplx( z1, 0.0e0 ), x11(i,i),
560 $ ldx11 )
561 ELSE
562 CALL cscal( p-i+1, cmplx( z1*cos(phi(i-1)), 0.0e0 ),
563 $ x11(i,i), ldx11 )
564 CALL caxpy( p-i+1, cmplx( -z1*z3*z4*sin(phi(i-1)),
565 $ 0.0e0 ), x12(i-1,i), ldx12, x11(i,i), ldx11 )
566 END IF
567 IF( i .EQ. 1 ) THEN
568 CALL cscal( m-p-i+1, cmplx( z2, 0.0e0 ), x21(i,i),
569 $ ldx21 )
570 ELSE
571 CALL cscal( m-p-i+1, cmplx( z2*cos(phi(i-1)), 0.0e0 ),
572 $ x21(i,i), ldx21 )
573 CALL caxpy( m-p-i+1, cmplx( -z2*z3*z4*sin(phi(i-1)),
574 $ 0.0e0 ), x22(i-1,i), ldx22, x21(i,i), ldx21 )
575 END IF
576*
577 theta(i) = atan2( scnrm2( m-p-i+1, x21(i,i), ldx21 ),
578 $ scnrm2( p-i+1, x11(i,i), ldx11 ) )
579*
580 CALL clacgv( p-i+1, x11(i,i), ldx11 )
581 CALL clacgv( m-p-i+1, x21(i,i), ldx21 )
582*
583 CALL clarfgp( 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 clarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
587 $ taup2(i) )
588 ELSE
589 CALL clarfgp( 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 clarf( 'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
595 $ x11(i+1,i), ldx11, work )
596 CALL clarf( 'R', m-q-i+1, p-i+1, x11(i,i), ldx11, taup1(i),
597 $ x12(i,i), ldx12, work )
598 CALL clarf( 'R', q-i, m-p-i+1, x21(i,i), ldx21, taup2(i),
599 $ x21(i+1,i), ldx21, work )
600 CALL clarf( 'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
601 $ taup2(i), x22(i,i), ldx22, work )
602*
603 CALL clacgv( p-i+1, x11(i,i), ldx11 )
604 CALL clacgv( m-p-i+1, x21(i,i), ldx21 )
605*
606 IF( i .LT. q ) THEN
607 CALL cscal( q-i, cmplx( -z1*z3*sin(theta(i)), 0.0e0 ),
608 $ x11(i+1,i), 1 )
609 CALL caxpy( q-i, cmplx( z2*z3*cos(theta(i)), 0.0e0 ),
610 $ x21(i+1,i), 1, x11(i+1,i), 1 )
611 END IF
612 CALL cscal( m-q-i+1, cmplx( -z1*z4*sin(theta(i)), 0.0e0 ),
613 $ x12(i,i), 1 )
614 CALL caxpy( m-q-i+1, cmplx( z2*z4*cos(theta(i)), 0.0e0 ),
615 $ x22(i,i), 1, x12(i,i), 1 )
616*
617 IF( i .LT. q )
618 $ phi(i) = atan2( scnrm2( q-i, x11(i+1,i), 1 ),
619 $ scnrm2( m-q-i+1, x12(i,i), 1 ) )
620*
621 IF( i .LT. q ) THEN
622 CALL clarfgp( q-i, x11(i+1,i), x11(i+2,i), 1, tauq1(i) )
623 x11(i+1,i) = one
624 END IF
625 CALL clarfgp( 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 clarf( 'L', q-i, p-i, x11(i+1,i), 1,
630 $ conjg(tauq1(i)), x11(i+1,i+1), ldx11, work )
631 CALL clarf( 'L', q-i, m-p-i, x11(i+1,i), 1,
632 $ conjg(tauq1(i)), x21(i+1,i+1), ldx21, work )
633 END IF
634 CALL clarf( 'L', m-q-i+1, p-i, x12(i,i), 1, conjg(tauq2(i)),
635 $ x12(i,i+1), ldx12, work )
636
637 IF ( m-p .GT. i ) THEN
638 CALL clarf( 'L', m-q-i+1, m-p-i, x12(i,i), 1,
639 $ conjg(tauq2(i)), x22(i,i+1), ldx22, work )
640 END IF
641 END DO
642*
643* Reduce columns Q + 1, ..., P of X12, X22
644*
645 DO i = q + 1, p
646*
647 CALL cscal( m-q-i+1, cmplx( -z1*z4, 0.0e0 ), x12(i,i), 1 )
648 CALL clarfgp( 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 clarf( 'L', m-q-i+1, p-i, x12(i,i), 1,
653 $ conjg(tauq2(i)), x12(i,i+1), ldx12, work )
654 END IF
655 IF( m-p-q .GE. 1 )
656 $ CALL clarf( 'L', m-q-i+1, m-p-q, x12(i,i), 1,
657 $ conjg(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 cscal( m-p-q-i+1, cmplx( z2*z4, 0.0e0 ),
666 $ x22(p+i,q+i), 1 )
667 CALL clarfgp( 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 IF ( m-p-q .NE. i ) THEN
671 CALL clarf( 'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i), 1,
672 $ conjg(tauq2(p+i)), x22(p+i,q+i+1), ldx22,
673 $ work )
674 END IF
675 END DO
676*
677 END IF
678*
679 RETURN
680*
681* End of CUNBDB
682*
683 END
684
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
subroutine clarf(side, m, n, v, incv, tau, c, ldc, work)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition clarf.f:128
subroutine clarfgp(n, alpha, x, incx, tau)
CLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition clarfgp.f:104
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
subroutine cunbdb(trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21, ldx21, x22, ldx22, theta, phi, taup1, taup2, tauq1, tauq2, work, lwork, info)
CUNBDB
Definition cunbdb.f:287