LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zuncsd.f
Go to the documentation of this file.
1*> \brief \b ZUNCSD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZUNCSD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zuncsd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zuncsd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zuncsd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* RECURSIVE SUBROUTINE ZUNCSD( 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* DOUBLE PRECISION THETA( * )
36* DOUBLE PRECISION RWORK( * )
37* COMPLEX*16 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*> ZUNCSD 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*16 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*16 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*16 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*16 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 DOUBLE PRECISION 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*16 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*16 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*16 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*16 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*16 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 DOUBLE PRECISION 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: ZBBCSD 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 complex16OTHERcomputational
312*
313* =====================================================================
314 RECURSIVE SUBROUTINE zuncsd( 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 DOUBLE PRECISION theta( * )
333 DOUBLE PRECISION rwork( * )
334 COMPLEX*16 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*16 one, zero
344 PARAMETER ( one = (1.0d0,0.0d0),
345 $ zero = (0.0d0,0.0d0) )
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, zbbcsd, zlacpy, zlapmr, zlapmt,
365* ..
366* .. External Functions ..
367 LOGICAL lsame
368 EXTERNAL lsame
369* ..
370* .. Intrinsic Functions
371 INTRINSIC int, max, min
372* ..
373* .. Executable Statements ..
374*
375* Test input arguments
376*
377 info = 0
378 wantu1 = lsame( jobu1, 'Y' )
379 wantu2 = lsame( jobu2, 'Y' )
380 wantv1t = lsame( jobv1t, 'Y' )
381 wantv2t = lsame( jobv2t, 'Y' )
382 colmajor = .NOT. lsame( trans, 'T' )
383 defaultsigns = .NOT. lsame( signs, 'O' )
384 lquery = lwork .EQ. -1
385 lrquery = lrwork .EQ. -1
386 IF( m .LT. 0 ) THEN
387 info = -7
388 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
389 info = -8
390 ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
391 info = -9
392 ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
393 info = -11
394 ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
395 info = -11
396 ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
397 info = -13
398 ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
399 info = -13
400 ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
401 info = -15
402 ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
403 info = -15
404 ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
405 info = -17
406 ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
407 info = -17
408 ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
409 info = -20
410 ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
411 info = -22
412 ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
413 info = -24
414 ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
415 info = -26
416 END IF
417*
418* Work with transpose if convenient
419*
420 IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) ) THEN
421 IF( colmajor ) THEN
422 transt = 'T'
423 ELSE
424 transt = 'N'
425 END IF
426 IF( defaultsigns ) THEN
427 signst = 'O'
428 ELSE
429 signst = 'D'
430 END IF
431 CALL zuncsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
432 $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
433 $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
434 $ u2, ldu2, work, lwork, rwork, lrwork, iwork,
435 $ info )
436 RETURN
437 END IF
438*
439* Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
440* convenient
441*
442 IF( info .EQ. 0 .AND. m-q .LT. q ) THEN
443 IF( defaultsigns ) THEN
444 signst = 'O'
445 ELSE
446 signst = 'D'
447 END IF
448 CALL zuncsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
449 $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
450 $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
451 $ ldv1t, work, lwork, rwork, lrwork, iwork, info )
452 RETURN
453 END IF
454*
455* Compute workspace
456*
457 IF( info .EQ. 0 ) THEN
458*
459* Real workspace
460*
461 iphi = 2
462 ib11d = iphi + max( 1, q - 1 )
463 ib11e = ib11d + max( 1, q )
464 ib12d = ib11e + max( 1, q - 1 )
465 ib12e = ib12d + max( 1, q )
466 ib21d = ib12e + max( 1, q - 1 )
467 ib21e = ib21d + max( 1, q )
468 ib22d = ib21e + max( 1, q - 1 )
469 ib22e = ib22d + max( 1, q )
470 ibbcsd = ib22e + max( 1, q - 1 )
471 CALL zbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
472 $ theta, theta, u1, ldu1, u2, ldu2, v1t, ldv1t,
473 $ v2t, ldv2t, theta, theta, theta, theta, theta,
474 $ theta, theta, theta, rwork, -1, childinfo )
475 lbbcsdworkopt = int( rwork(1) )
476 lbbcsdworkmin = lbbcsdworkopt
477 lrworkopt = ibbcsd + lbbcsdworkopt - 1
478 lrworkmin = ibbcsd + lbbcsdworkmin - 1
479 rwork(1) = lrworkopt
480*
481* Complex workspace
482*
483 itaup1 = 2
484 itaup2 = itaup1 + max( 1, p )
485 itauq1 = itaup2 + max( 1, m - p )
486 itauq2 = itauq1 + max( 1, q )
487 iorgqr = itauq2 + max( 1, m - q )
488 CALL zungqr( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
489 $ childinfo )
490 lorgqrworkopt = int( work(1) )
491 lorgqrworkmin = max( 1, m - q )
492 iorglq = itauq2 + max( 1, m - q )
493 CALL zunglq( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
494 $ childinfo )
495 lorglqworkopt = int( work(1) )
496 lorglqworkmin = max( 1, m - q )
497 iorbdb = itauq2 + max( 1, m - q )
498 CALL zunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
499 $ x21, ldx21, x22, ldx22, theta, theta, u1, u2,
500 $ v1t, v2t, work, -1, childinfo )
501 lorbdbworkopt = int( work(1) )
502 lorbdbworkmin = lorbdbworkopt
503 lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
504 $ iorbdb + lorbdbworkopt ) - 1
505 lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
506 $ iorbdb + lorbdbworkmin ) - 1
507 work(1) = max(lworkopt,lworkmin)
508*
509 IF( lwork .LT. lworkmin
510 $ .AND. .NOT. ( lquery .OR. lrquery ) ) THEN
511 info = -22
512 ELSE IF( lrwork .LT. lrworkmin
513 $ .AND. .NOT. ( lquery .OR. lrquery ) ) THEN
514 info = -24
515 ELSE
516 lorgqrwork = lwork - iorgqr + 1
517 lorglqwork = lwork - iorglq + 1
518 lorbdbwork = lwork - iorbdb + 1
519 lbbcsdwork = lrwork - ibbcsd + 1
520 END IF
521 END IF
522*
523* Abort if any illegal arguments
524*
525 IF( info .NE. 0 ) THEN
526 CALL xerbla( 'ZUNCSD', -info )
527 RETURN
528 ELSE IF( lquery .OR. lrquery ) THEN
529 RETURN
530 END IF
531*
532* Transform to bidiagonal block form
533*
534 CALL zunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
535 $ ldx21, x22, ldx22, theta, rwork(iphi), work(itaup1),
536 $ work(itaup2), work(itauq1), work(itauq2),
537 $ work(iorbdb), lorbdbwork, childinfo )
538*
539* Accumulate Householder reflectors
540*
541 IF( colmajor ) THEN
542 IF( wantu1 .AND. p .GT. 0 ) THEN
543 CALL zlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
544 CALL zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
545 $ lorgqrwork, info)
546 END IF
547 IF( wantu2 .AND. m-p .GT. 0 ) THEN
548 CALL zlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
549 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
550 $ work(iorgqr), lorgqrwork, info )
551 END IF
552 IF( wantv1t .AND. q .GT. 0 ) THEN
553 CALL zlacpy( 'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
554 $ ldv1t )
555 v1t(1, 1) = one
556 DO j = 2, q
557 v1t(1,j) = zero
558 v1t(j,1) = zero
559 END DO
560 CALL zunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
561 $ work(iorglq), lorglqwork, info )
562 END IF
563 IF( wantv2t .AND. m-q .GT. 0 ) THEN
564 CALL zlacpy( 'U', p, m-q, x12, ldx12, v2t, ldv2t )
565 IF( m-p .GT. q) THEN
566 CALL zlacpy( 'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
567 $ v2t(p+1,p+1), ldv2t )
568 END IF
569 IF( m .GT. q ) THEN
570 CALL zunglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
571 $ work(iorglq), lorglqwork, info )
572 END IF
573 END IF
574 ELSE
575 IF( wantu1 .AND. p .GT. 0 ) THEN
576 CALL zlacpy( 'U', q, p, x11, ldx11, u1, ldu1 )
577 CALL zunglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
578 $ lorglqwork, info)
579 END IF
580 IF( wantu2 .AND. m-p .GT. 0 ) THEN
581 CALL zlacpy( 'U', q, m-p, x21, ldx21, u2, ldu2 )
582 CALL zunglq( m-p, m-p, q, u2, ldu2, work(itaup2),
583 $ work(iorglq), lorglqwork, info )
584 END IF
585 IF( wantv1t .AND. q .GT. 0 ) THEN
586 CALL zlacpy( 'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
587 $ ldv1t )
588 v1t(1, 1) = one
589 DO j = 2, q
590 v1t(1,j) = zero
591 v1t(j,1) = zero
592 END DO
593 CALL zungqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
594 $ work(iorgqr), lorgqrwork, info )
595 END IF
596 IF( wantv2t .AND. m-q .GT. 0 ) THEN
597 p1 = min( p+1, m )
598 q1 = min( q+1, m )
599 CALL zlacpy( 'L', m-q, p, x12, ldx12, v2t, ldv2t )
600 IF( m .GT. p+q ) THEN
601 CALL zlacpy( 'L', m-p-q, m-p-q, x22(p1,q1), ldx22,
602 $ v2t(p+1,p+1), ldv2t )
603 END IF
604 CALL zungqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
605 $ work(iorgqr), lorgqrwork, info )
606 END IF
607 END IF
608*
609* Compute the CSD of the matrix in bidiagonal-block form
610*
611 CALL zbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
612 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
613 $ ldv2t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
614 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
615 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
616 $ lbbcsdwork, info )
617*
618* Permute rows and columns to place identity submatrices in top-
619* left corner of (1,1)-block and/or bottom-right corner of (1,2)-
620* block and/or bottom-right corner of (2,1)-block and/or top-left
621* corner of (2,2)-block
622*
623 IF( q .GT. 0 .AND. wantu2 ) THEN
624 DO i = 1, q
625 iwork(i) = m - p - q + i
626 END DO
627 DO i = q + 1, m - p
628 iwork(i) = i - q
629 END DO
630 IF( colmajor ) THEN
631 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
632 ELSE
633 CALL zlapmr( .false., m-p, m-p, u2, ldu2, iwork )
634 END IF
635 END IF
636 IF( m .GT. 0 .AND. wantv2t ) THEN
637 DO i = 1, p
638 iwork(i) = m - p - q + i
639 END DO
640 DO i = p + 1, m - q
641 iwork(i) = i - p
642 END DO
643 IF( .NOT. colmajor ) THEN
644 CALL zlapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
645 ELSE
646 CALL zlapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
647 END IF
648 END IF
649*
650 RETURN
651*
652* End ZUNCSD
653*
654 END
655
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zlapmr(FORWRD, M, N, X, LDX, K)
ZLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: zlapmr.f:104
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlapmt(FORWRD, M, N, X, LDX, K)
ZLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: zlapmt.f:104
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
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ
Definition: zunglq.f:127
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:128
recursive subroutine zuncsd(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)
ZUNCSD
Definition: zuncsd.f:320
subroutine zbbcsd(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)
ZBBCSD
Definition: zbbcsd.f:332