LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cuncsd.f
Go to the documentation of this file.
1*> \brief \b CUNCSD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CUNCSD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cuncsd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cuncsd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cuncsd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* RECURSIVE SUBROUTINE CUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
22* SIGNS, M, P, Q, X11, LDX11, X12,
23* LDX12, X21, LDX21, X22, LDX22, THETA,
24* U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
25* LDV2T, WORK, LWORK, RWORK, LRWORK,
26* IWORK, INFO )
27*
28* .. Scalar Arguments ..
29* CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
30* INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
31* $ LDX21, LDX22, LRWORK, LWORK, M, P, Q
32* ..
33* .. Array Arguments ..
34* INTEGER IWORK( * )
35* REAL THETA( * )
36* REAL RWORK( * )
37* COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
38* $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
39* $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
40* $ * )
41* ..
42*
43*
44*> \par Purpose:
45* =============
46*>
47*> \verbatim
48*>
49*> CUNCSD computes the CS decomposition of an M-by-M partitioned
50*> unitary matrix X:
51*>
52*> [ I 0 0 | 0 0 0 ]
53*> [ 0 C 0 | 0 -S 0 ]
54*> [ X11 | X12 ] [ U1 | ] [ 0 0 0 | 0 0 -I ] [ V1 | ]**H
55*> X = [-----------] = [---------] [---------------------] [---------] .
56*> [ X21 | X22 ] [ | U2 ] [ 0 0 0 | I 0 0 ] [ | V2 ]
57*> [ 0 S 0 | 0 C 0 ]
58*> [ 0 0 I | 0 0 0 ]
59*>
60*> X11 is P-by-Q. The unitary matrices U1, U2, V1, and V2 are P-by-P,
61*> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
62*> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
63*> which R = MIN(P,M-P,Q,M-Q).
64*> \endverbatim
65*
66* Arguments:
67* ==========
68*
69*> \param[in] JOBU1
70*> \verbatim
71*> JOBU1 is CHARACTER
72*> = 'Y': U1 is computed;
73*> otherwise: U1 is not computed.
74*> \endverbatim
75*>
76*> \param[in] JOBU2
77*> \verbatim
78*> JOBU2 is CHARACTER
79*> = 'Y': U2 is computed;
80*> otherwise: U2 is not computed.
81*> \endverbatim
82*>
83*> \param[in] JOBV1T
84*> \verbatim
85*> JOBV1T is CHARACTER
86*> = 'Y': V1T is computed;
87*> otherwise: V1T is not computed.
88*> \endverbatim
89*>
90*> \param[in] JOBV2T
91*> \verbatim
92*> JOBV2T is CHARACTER
93*> = 'Y': V2T is computed;
94*> otherwise: V2T is not computed.
95*> \endverbatim
96*>
97*> \param[in] TRANS
98*> \verbatim
99*> TRANS is CHARACTER
100*> = 'T': X, U1, U2, V1T, and V2T are stored in row-major
101*> order;
102*> otherwise: X, U1, U2, V1T, and V2T are stored in column-
103*> major order.
104*> \endverbatim
105*>
106*> \param[in] SIGNS
107*> \verbatim
108*> SIGNS is CHARACTER
109*> = 'O': The lower-left block is made nonpositive (the
110*> "other" convention);
111*> otherwise: The upper-right block is made nonpositive (the
112*> "default" convention).
113*> \endverbatim
114*>
115*> \param[in] M
116*> \verbatim
117*> M is INTEGER
118*> The number of rows and columns in X.
119*> \endverbatim
120*>
121*> \param[in] P
122*> \verbatim
123*> P is INTEGER
124*> The number of rows in X11 and X12. 0 <= P <= M.
125*> \endverbatim
126*>
127*> \param[in] Q
128*> \verbatim
129*> Q is INTEGER
130*> The number of columns in X11 and X21. 0 <= Q <= M.
131*> \endverbatim
132*>
133*> \param[in,out] X11
134*> \verbatim
135*> X11 is COMPLEX array, dimension (LDX11,Q)
136*> On entry, part of the unitary matrix whose CSD is desired.
137*> \endverbatim
138*>
139*> \param[in] LDX11
140*> \verbatim
141*> LDX11 is INTEGER
142*> The leading dimension of X11. LDX11 >= MAX(1,P).
143*> \endverbatim
144*>
145*> \param[in,out] X12
146*> \verbatim
147*> X12 is COMPLEX array, dimension (LDX12,M-Q)
148*> On entry, part of the unitary matrix whose CSD is desired.
149*> \endverbatim
150*>
151*> \param[in] LDX12
152*> \verbatim
153*> LDX12 is INTEGER
154*> The leading dimension of X12. LDX12 >= MAX(1,P).
155*> \endverbatim
156*>
157*> \param[in,out] X21
158*> \verbatim
159*> X21 is COMPLEX array, dimension (LDX21,Q)
160*> On entry, part of the unitary matrix whose CSD is desired.
161*> \endverbatim
162*>
163*> \param[in] LDX21
164*> \verbatim
165*> LDX21 is INTEGER
166*> The leading dimension of X11. LDX21 >= MAX(1,M-P).
167*> \endverbatim
168*>
169*> \param[in,out] X22
170*> \verbatim
171*> X22 is COMPLEX array, dimension (LDX22,M-Q)
172*> On entry, part of the unitary matrix whose CSD is desired.
173*> \endverbatim
174*>
175*> \param[in] LDX22
176*> \verbatim
177*> LDX22 is INTEGER
178*> The leading dimension of X11. LDX22 >= MAX(1,M-P).
179*> \endverbatim
180*>
181*> \param[out] THETA
182*> \verbatim
183*> THETA is REAL array, dimension (R), in which R =
184*> MIN(P,M-P,Q,M-Q).
185*> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
186*> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
187*> \endverbatim
188*>
189*> \param[out] U1
190*> \verbatim
191*> U1 is COMPLEX array, dimension (LDU1,P)
192*> If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1.
193*> \endverbatim
194*>
195*> \param[in] LDU1
196*> \verbatim
197*> LDU1 is INTEGER
198*> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
199*> MAX(1,P).
200*> \endverbatim
201*>
202*> \param[out] U2
203*> \verbatim
204*> U2 is COMPLEX array, dimension (LDU2,M-P)
205*> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
206*> matrix U2.
207*> \endverbatim
208*>
209*> \param[in] LDU2
210*> \verbatim
211*> LDU2 is INTEGER
212*> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
213*> MAX(1,M-P).
214*> \endverbatim
215*>
216*> \param[out] V1T
217*> \verbatim
218*> V1T is COMPLEX array, dimension (LDV1T,Q)
219*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
220*> matrix V1**H.
221*> \endverbatim
222*>
223*> \param[in] LDV1T
224*> \verbatim
225*> LDV1T is INTEGER
226*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
227*> MAX(1,Q).
228*> \endverbatim
229*>
230*> \param[out] V2T
231*> \verbatim
232*> V2T is COMPLEX array, dimension (LDV2T,M-Q)
233*> If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) unitary
234*> matrix V2**H.
235*> \endverbatim
236*>
237*> \param[in] LDV2T
238*> \verbatim
239*> LDV2T is INTEGER
240*> The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
241*> MAX(1,M-Q).
242*> \endverbatim
243*>
244*> \param[out] WORK
245*> \verbatim
246*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
247*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
248*> \endverbatim
249*>
250*> \param[in] LWORK
251*> \verbatim
252*> LWORK is INTEGER
253*> The dimension of the array WORK.
254*>
255*> If LWORK = -1, then a workspace query is assumed; the routine
256*> only calculates the optimal size of the WORK array, returns
257*> this value as the first entry of the work array, and no error
258*> message related to LWORK is issued by XERBLA.
259*> \endverbatim
260*>
261*> \param[out] RWORK
262*> \verbatim
263*> RWORK is REAL array, dimension MAX(1,LRWORK)
264*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
265*> If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
266*> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
267*> define the matrix in intermediate bidiagonal-block form
268*> remaining after nonconvergence. INFO specifies the number
269*> of nonzero PHI's.
270*> \endverbatim
271*>
272*> \param[in] LRWORK
273*> \verbatim
274*> LRWORK is INTEGER
275*> The dimension of the array RWORK.
276*>
277*> If LRWORK = -1, then a workspace query is assumed; the routine
278*> only calculates the optimal size of the RWORK array, returns
279*> this value as the first entry of the work array, and no error
280*> message related to LRWORK is issued by XERBLA.
281*> \endverbatim
282*>
283*> \param[out] IWORK
284*> \verbatim
285*> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
286*> \endverbatim
287*>
288*> \param[out] INFO
289*> \verbatim
290*> INFO is INTEGER
291*> = 0: successful exit.
292*> < 0: if INFO = -i, the i-th argument had an illegal value.
293*> > 0: CBBCSD did not converge. See the description of RWORK
294*> above for details.
295*> \endverbatim
296*
297*> \par References:
298* ================
299*>
300*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
301*> Algorithms, 50(1):33-65, 2009.
302*
303* Authors:
304* ========
305*
306*> \author Univ. of Tennessee
307*> \author Univ. of California Berkeley
308*> \author Univ. of Colorado Denver
309*> \author NAG Ltd.
310*
311*> \ingroup uncsd
312*
313* =====================================================================
314 RECURSIVE SUBROUTINE cuncsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
315 $ SIGNS, M, P, Q, X11, LDX11, X12,
316 $ LDX12, X21, LDX21, X22, LDX22, THETA,
317 $ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
318 $ LDV2T, WORK, LWORK, RWORK, LRWORK,
319 $ IWORK, INFO )
320*
321* -- LAPACK computational routine --
322* -- LAPACK is a software package provided by Univ. of Tennessee, --
323* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
324*
325* .. Scalar Arguments ..
326 CHARACTER jobu1, jobu2, jobv1t, jobv2t, signs, trans
327 INTEGER info, ldu1, ldu2, ldv1t, ldv2t, ldx11, ldx12,
328 $ ldx21, ldx22, lrwork, lwork, m, p, q
329* ..
330* .. Array Arguments ..
331 INTEGER iwork( * )
332 REAL theta( * )
333 REAL rwork( * )
334 COMPLEX u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
335 $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
336 $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
337 $ * )
338* ..
339*
340* ===================================================================
341*
342* .. Parameters ..
343 COMPLEX one, zero
344 PARAMETER ( one = (1.0e0,0.0e0),
345 $ zero = (0.0e0,0.0e0) )
346* ..
347* .. Local Scalars ..
348 CHARACTER transt, signst
349 INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
350 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
351 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
352 $ itauq2, j, lbbcsdwork, lbbcsdworkmin,
353 $ lbbcsdworkopt, lorbdbwork, lorbdbworkmin,
354 $ lorbdbworkopt, lorglqwork, lorglqworkmin,
355 $ lorglqworkopt, lorgqrwork, lorgqrworkmin,
356 $ lorgqrworkopt, lworkmin, lworkopt, p1, q1
357 LOGICAL colmajor, defaultsigns, lquery, wantu1, wantu2,
358 $ wantv1t, wantv2t
359 INTEGER lrworkmin, lrworkopt
360 LOGICAL lrquery
361* ..
362* .. External Subroutines ..
363 EXTERNAL xerbla, cbbcsd, clacpy, clapmr, clapmt,
365* ..
366* .. External Functions ..
367 LOGICAL lsame
368 REAL sroundup_lwork
369 EXTERNAL lsame, sroundup_lwork
370* ..
371* .. Intrinsic Functions
372 INTRINSIC int, max, min
373* ..
374* .. Executable Statements ..
375*
376* Test input arguments
377*
378 info = 0
379 wantu1 = lsame( jobu1, 'Y' )
380 wantu2 = lsame( jobu2, 'Y' )
381 wantv1t = lsame( jobv1t, 'Y' )
382 wantv2t = lsame( jobv2t, 'Y' )
383 colmajor = .NOT. lsame( trans, 'T' )
384 defaultsigns = .NOT. lsame( signs, 'O' )
385 lquery = lwork .EQ. -1
386 lrquery = lrwork .EQ. -1
387 IF( m .LT. 0 ) THEN
388 info = -7
389 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
390 info = -8
391 ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
392 info = -9
393 ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
394 info = -11
395 ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
396 info = -11
397 ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
398 info = -13
399 ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
400 info = -13
401 ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
402 info = -15
403 ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
404 info = -15
405 ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
406 info = -17
407 ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
408 info = -17
409 ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
410 info = -20
411 ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
412 info = -22
413 ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
414 info = -24
415 ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
416 info = -26
417 END IF
418*
419* Work with transpose if convenient
420*
421 IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) ) THEN
422 IF( colmajor ) THEN
423 transt = 'T'
424 ELSE
425 transt = 'N'
426 END IF
427 IF( defaultsigns ) THEN
428 signst = 'O'
429 ELSE
430 signst = 'D'
431 END IF
432 CALL cuncsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
433 $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
434 $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
435 $ u2, ldu2, work, lwork, rwork, lrwork, iwork,
436 $ info )
437 RETURN
438 END IF
439*
440* Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
441* convenient
442*
443 IF( info .EQ. 0 .AND. m-q .LT. q ) THEN
444 IF( defaultsigns ) THEN
445 signst = 'O'
446 ELSE
447 signst = 'D'
448 END IF
449 CALL cuncsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
450 $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
451 $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
452 $ ldv1t, work, lwork, rwork, lrwork, iwork, info )
453 RETURN
454 END IF
455*
456* Compute workspace
457*
458 IF( info .EQ. 0 ) THEN
459*
460* Real workspace
461*
462 iphi = 2
463 ib11d = iphi + max( 1, q - 1 )
464 ib11e = ib11d + max( 1, q )
465 ib12d = ib11e + max( 1, q - 1 )
466 ib12e = ib12d + max( 1, q )
467 ib21d = ib12e + max( 1, q - 1 )
468 ib21e = ib21d + max( 1, q )
469 ib22d = ib21e + max( 1, q - 1 )
470 ib22e = ib22d + max( 1, q )
471 ibbcsd = ib22e + max( 1, q - 1 )
472 CALL cbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
473 $ theta, theta, u1, ldu1, u2, ldu2, v1t, ldv1t,
474 $ v2t, ldv2t, theta, theta, theta, theta, theta,
475 $ theta, theta, theta, rwork, -1, childinfo )
476 lbbcsdworkopt = int( rwork(1) )
477 lbbcsdworkmin = lbbcsdworkopt
478 lrworkopt = ibbcsd + lbbcsdworkopt - 1
479 lrworkmin = ibbcsd + lbbcsdworkmin - 1
480 rwork(1) = lrworkopt
481*
482* Complex workspace
483*
484 itaup1 = 2
485 itaup2 = itaup1 + max( 1, p )
486 itauq1 = itaup2 + max( 1, m - p )
487 itauq2 = itauq1 + max( 1, q )
488 iorgqr = itauq2 + max( 1, m - q )
489 CALL cungqr( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
490 $ childinfo )
491 lorgqrworkopt = int( work(1) )
492 lorgqrworkmin = max( 1, m - q )
493 iorglq = itauq2 + max( 1, m - q )
494 CALL cunglq( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
495 $ childinfo )
496 lorglqworkopt = int( work(1) )
497 lorglqworkmin = max( 1, m - q )
498 iorbdb = itauq2 + max( 1, m - q )
499 CALL cunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
500 $ x21, ldx21, x22, ldx22, theta, theta, u1, u2,
501 $ v1t, v2t, work, -1, childinfo )
502 lorbdbworkopt = int( work(1) )
503 lorbdbworkmin = lorbdbworkopt
504 lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
505 $ iorbdb + lorbdbworkopt ) - 1
506 lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
507 $ iorbdb + lorbdbworkmin ) - 1
508 lworkopt = max(lworkopt,lworkmin)
509 work(1) = sroundup_lwork(lworkopt)
510*
511 IF( lwork .LT. lworkmin
512 $ .AND. .NOT. ( lquery .OR. lrquery ) ) THEN
513 info = -22
514 ELSE IF( lrwork .LT. lrworkmin
515 $ .AND. .NOT. ( lquery .OR. lrquery ) ) THEN
516 info = -24
517 ELSE
518 lorgqrwork = lwork - iorgqr + 1
519 lorglqwork = lwork - iorglq + 1
520 lorbdbwork = lwork - iorbdb + 1
521 lbbcsdwork = lrwork - ibbcsd + 1
522 END IF
523 END IF
524*
525* Abort if any illegal arguments
526*
527 IF( info .NE. 0 ) THEN
528 CALL xerbla( 'CUNCSD', -info )
529 RETURN
530 ELSE IF( lquery .OR. lrquery ) THEN
531 RETURN
532 END IF
533*
534* Transform to bidiagonal block form
535*
536 CALL cunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
537 $ ldx21, x22, ldx22, theta, rwork(iphi), work(itaup1),
538 $ work(itaup2), work(itauq1), work(itauq2),
539 $ work(iorbdb), lorbdbwork, childinfo )
540*
541* Accumulate Householder reflectors
542*
543 IF( colmajor ) THEN
544 IF( wantu1 .AND. p .GT. 0 ) THEN
545 CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
546 CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
547 $ lorgqrwork, info)
548 END IF
549 IF( wantu2 .AND. m-p .GT. 0 ) THEN
550 CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
551 CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
552 $ work(iorgqr), lorgqrwork, info )
553 END IF
554 IF( wantv1t .AND. q .GT. 0 ) THEN
555 CALL clacpy( 'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
556 $ ldv1t )
557 v1t(1, 1) = one
558 DO j = 2, q
559 v1t(1,j) = zero
560 v1t(j,1) = zero
561 END DO
562 CALL cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
563 $ work(iorglq), lorglqwork, info )
564 END IF
565 IF( wantv2t .AND. m-q .GT. 0 ) THEN
566 CALL clacpy( 'U', p, m-q, x12, ldx12, v2t, ldv2t )
567 IF( m-p .GT. q ) THEN
568 CALL clacpy( 'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
569 $ v2t(p+1,p+1), ldv2t )
570 END IF
571 IF( m .GT. q ) THEN
572 CALL cunglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
573 $ work(iorglq), lorglqwork, info )
574 END IF
575 END IF
576 ELSE
577 IF( wantu1 .AND. p .GT. 0 ) THEN
578 CALL clacpy( 'U', q, p, x11, ldx11, u1, ldu1 )
579 CALL cunglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
580 $ lorglqwork, info)
581 END IF
582 IF( wantu2 .AND. m-p .GT. 0 ) THEN
583 CALL clacpy( 'U', q, m-p, x21, ldx21, u2, ldu2 )
584 CALL cunglq( m-p, m-p, q, u2, ldu2, work(itaup2),
585 $ work(iorglq), lorglqwork, info )
586 END IF
587 IF( wantv1t .AND. q .GT. 0 ) THEN
588 CALL clacpy( 'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
589 $ ldv1t )
590 v1t(1, 1) = one
591 DO j = 2, q
592 v1t(1,j) = zero
593 v1t(j,1) = zero
594 END DO
595 CALL cungqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
596 $ work(iorgqr), lorgqrwork, info )
597 END IF
598 IF( wantv2t .AND. m-q .GT. 0 ) THEN
599 p1 = min( p+1, m )
600 q1 = min( q+1, m )
601 CALL clacpy( 'L', m-q, p, x12, ldx12, v2t, ldv2t )
602 IF ( m .GT. p+q ) THEN
603 CALL clacpy( 'L', m-p-q, m-p-q, x22(p1,q1), ldx22,
604 $ v2t(p+1,p+1), ldv2t )
605 END IF
606 CALL cungqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
607 $ work(iorgqr), lorgqrwork, info )
608 END IF
609 END IF
610*
611* Compute the CSD of the matrix in bidiagonal-block form
612*
613 CALL cbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
614 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
615 $ ldv2t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
616 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
617 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
618 $ lbbcsdwork, info )
619*
620* Permute rows and columns to place identity submatrices in top-
621* left corner of (1,1)-block and/or bottom-right corner of (1,2)-
622* block and/or bottom-right corner of (2,1)-block and/or top-left
623* corner of (2,2)-block
624*
625 IF( q .GT. 0 .AND. wantu2 ) THEN
626 DO i = 1, q
627 iwork(i) = m - p - q + i
628 END DO
629 DO i = q + 1, m - p
630 iwork(i) = i - q
631 END DO
632 IF( colmajor ) THEN
633 CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
634 ELSE
635 CALL clapmr( .false., m-p, m-p, u2, ldu2, iwork )
636 END IF
637 END IF
638 IF( m .GT. 0 .AND. wantv2t ) THEN
639 DO i = 1, p
640 iwork(i) = m - p - q + i
641 END DO
642 DO i = p + 1, m - q
643 iwork(i) = i - p
644 END DO
645 IF( .NOT. colmajor ) THEN
646 CALL clapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
647 ELSE
648 CALL clapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
649 END IF
650 END IF
651*
652 RETURN
653*
654* End CUNCSD
655*
656 END
657
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cbbcsd(jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta, phi, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t, b11d, b11e, b12d, b12e, b21d, b21e, b22d, b22e, rwork, lrwork, info)
CBBCSD
Definition cbbcsd.f:332
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
subroutine clapmr(forwrd, m, n, x, ldx, k)
CLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition clapmr.f:104
subroutine clapmt(forwrd, m, n, x, ldx, k)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition clapmt.f:104
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
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
recursive subroutine cuncsd(jobu1, jobu2, jobv1t, jobv2t, trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21, ldx21, x22, ldx22, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t, work, lwork, rwork, lrwork, iwork, info)
CUNCSD
Definition cuncsd.f:320
subroutine cunglq(m, n, k, a, lda, tau, work, lwork, info)
CUNGLQ
Definition cunglq.f:127
subroutine cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
Definition cungqr.f:128