LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zuncsd2by1.f
Go to the documentation of this file.
1*> \brief \b ZUNCSD2BY1
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download ZUNCSD2BY1 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zuncsd2by1.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zuncsd2by1.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zuncsd2by1.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
20* X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
21* LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK,
22* INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER JOBU1, JOBU2, JOBV1T
26* INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
27* $ M, P, Q
28* INTEGER LRWORK, LRWORKMIN, LRWORKOPT
29* ..
30* .. Array Arguments ..
31* DOUBLE PRECISION RWORK(*)
32* DOUBLE PRECISION THETA(*)
33* COMPLEX*16 U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
34* $ X11(LDX11,*), X21(LDX21,*)
35* INTEGER IWORK(*)
36* ..
37*
38*
39*> \par Purpose:
40* =============
41*>
42*>\verbatim
43*>
44*> ZUNCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with
45*> orthonormal columns that has been partitioned into a 2-by-1 block
46*> structure:
47*>
48*> [ I1 0 0 ]
49*> [ 0 C 0 ]
50*> [ X11 ] [ U1 | ] [ 0 0 0 ]
51*> X = [-----] = [---------] [----------] V1**T .
52*> [ X21 ] [ | U2 ] [ 0 0 0 ]
53*> [ 0 S 0 ]
54*> [ 0 0 I2]
55*>
56*> X11 is P-by-Q. The unitary matrices U1, U2, and V1 are P-by-P,
57*> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R
58*> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which
59*> R = MIN(P,M-P,Q,M-Q). I1 is a K1-by-K1 identity matrix and I2 is a
60*> K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0).
61*> \endverbatim
62*
63* Arguments:
64* ==========
65*
66*> \param[in] JOBU1
67*> \verbatim
68*> JOBU1 is CHARACTER
69*> = 'Y': U1 is computed;
70*> otherwise: U1 is not computed.
71*> \endverbatim
72*>
73*> \param[in] JOBU2
74*> \verbatim
75*> JOBU2 is CHARACTER
76*> = 'Y': U2 is computed;
77*> otherwise: U2 is not computed.
78*> \endverbatim
79*>
80*> \param[in] JOBV1T
81*> \verbatim
82*> JOBV1T is CHARACTER
83*> = 'Y': V1T is computed;
84*> otherwise: V1T is not computed.
85*> \endverbatim
86*>
87*> \param[in] M
88*> \verbatim
89*> M is INTEGER
90*> The number of rows in X.
91*> \endverbatim
92*>
93*> \param[in] P
94*> \verbatim
95*> P is INTEGER
96*> The number of rows in X11. 0 <= P <= M.
97*> \endverbatim
98*>
99*> \param[in] Q
100*> \verbatim
101*> Q is INTEGER
102*> The number of columns in X11 and X21. 0 <= Q <= M.
103*> \endverbatim
104*>
105*> \param[in,out] X11
106*> \verbatim
107*> X11 is COMPLEX*16 array, dimension (LDX11,Q)
108*> On entry, part of the unitary matrix whose CSD is desired.
109*> \endverbatim
110*>
111*> \param[in] LDX11
112*> \verbatim
113*> LDX11 is INTEGER
114*> The leading dimension of X11. LDX11 >= MAX(1,P).
115*> \endverbatim
116*>
117*> \param[in,out] X21
118*> \verbatim
119*> X21 is COMPLEX*16 array, dimension (LDX21,Q)
120*> On entry, part of the unitary matrix whose CSD is desired.
121*> \endverbatim
122*>
123*> \param[in] LDX21
124*> \verbatim
125*> LDX21 is INTEGER
126*> The leading dimension of X21. LDX21 >= MAX(1,M-P).
127*> \endverbatim
128*>
129*> \param[out] THETA
130*> \verbatim
131*> THETA is DOUBLE PRECISION array, dimension (R), in which R =
132*> MIN(P,M-P,Q,M-Q).
133*> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
134*> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
135*> \endverbatim
136*>
137*> \param[out] U1
138*> \verbatim
139*> U1 is COMPLEX*16 array, dimension (P)
140*> If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1.
141*> \endverbatim
142*>
143*> \param[in] LDU1
144*> \verbatim
145*> LDU1 is INTEGER
146*> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
147*> MAX(1,P).
148*> \endverbatim
149*>
150*> \param[out] U2
151*> \verbatim
152*> U2 is COMPLEX*16 array, dimension (M-P)
153*> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
154*> matrix U2.
155*> \endverbatim
156*>
157*> \param[in] LDU2
158*> \verbatim
159*> LDU2 is INTEGER
160*> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
161*> MAX(1,M-P).
162*> \endverbatim
163*>
164*> \param[out] V1T
165*> \verbatim
166*> V1T is COMPLEX*16 array, dimension (Q)
167*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
168*> matrix V1**T.
169*> \endverbatim
170*>
171*> \param[in] LDV1T
172*> \verbatim
173*> LDV1T is INTEGER
174*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
175*> MAX(1,Q).
176*> \endverbatim
177*>
178*> \param[out] WORK
179*> \verbatim
180*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
181*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182*> \endverbatim
183*>
184*> \param[in] LWORK
185*> \verbatim
186*> LWORK is INTEGER
187*> The dimension of the array WORK.
188*>
189*> If LWORK = -1, then a workspace query is assumed; the routine
190*> only calculates the optimal size of the WORK and RWORK
191*> arrays, returns this value as the first entry of the WORK
192*> and RWORK array, respectively, and no error message related
193*> to LWORK or LRWORK is issued by XERBLA.
194*> \endverbatim
195*>
196*> \param[out] RWORK
197*> \verbatim
198*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
199*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
200*> If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
201*> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
202*> define the matrix in intermediate bidiagonal-block form
203*> remaining after nonconvergence. INFO specifies the number
204*> of nonzero PHI's.
205*> \endverbatim
206*>
207*> \param[in] LRWORK
208*> \verbatim
209*> LRWORK is INTEGER
210*> The dimension of the array RWORK.
211*>
212*> If LRWORK=-1, then a workspace query is assumed; the routine
213*> only calculates the optimal size of the WORK and RWORK
214*> arrays, returns this value as the first entry of the WORK
215*> and RWORK array, respectively, and no error message related
216*> to LWORK or LRWORK is issued by XERBLA.
217*> \endverbatim
218*>
219*> \param[out] IWORK
220*> \verbatim
221*> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
222*> \endverbatim
223*>
224*> \param[out] INFO
225*> \verbatim
226*> INFO is INTEGER
227*> = 0: successful exit.
228*> < 0: if INFO = -i, the i-th argument had an illegal value.
229*> > 0: ZBBCSD did not converge. See the description of WORK
230*> above for details.
231*> \endverbatim
232*
233*> \par References:
234* ================
235*>
236*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
237*> Algorithms, 50(1):33-65, 2009.
238*
239* Authors:
240* ========
241*
242*> \author Univ. of Tennessee
243*> \author Univ. of California Berkeley
244*> \author Univ. of Colorado Denver
245*> \author NAG Ltd.
246*
247*> \ingroup uncsd2by1
248*
249* =====================================================================
250 SUBROUTINE zuncsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11,
251 $ LDX11,
252 $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
253 $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK,
254 $ INFO )
255*
256* -- LAPACK computational routine --
257* -- LAPACK is a software package provided by Univ. of Tennessee, --
258* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
259*
260* .. Scalar Arguments ..
261 CHARACTER JOBU1, JOBU2, JOBV1T
262 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
263 $ M, P, Q
264 INTEGER LRWORK, LRWORKMIN, LRWORKOPT
265* ..
266* .. Array Arguments ..
267 DOUBLE PRECISION RWORK(*)
268 DOUBLE PRECISION THETA(*)
269 COMPLEX*16 U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
270 $ X11(LDX11,*), X21(LDX21,*)
271 INTEGER IWORK(*)
272* ..
273*
274* =====================================================================
275*
276* .. Parameters ..
277 COMPLEX*16 ONE, ZERO
278 PARAMETER ( ONE = (1.0d0,0.0d0), zero = (0.0d0,0.0d0) )
279* ..
280* .. Local Scalars ..
281 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
282 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
283 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
284 $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
285 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
286 $ lworkmin, lworkopt, r
287 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
288* ..
289* .. Local Arrays ..
290 DOUBLE PRECISION DUM( 1 )
291 COMPLEX*16 CDUM( 1, 1 )
292* ..
293* .. External Subroutines ..
294 EXTERNAL zbbcsd, zcopy, zlacpy, zlapmr, zlapmt,
295 $ zunbdb1,
296 $ zunbdb2, zunbdb3, zunbdb4, zunglq, zungqr,
297 $ xerbla
298* ..
299* .. External Functions ..
300 LOGICAL LSAME
301 EXTERNAL LSAME
302* ..
303* .. Intrinsic Function ..
304 INTRINSIC int, max, min
305* ..
306* .. Executable Statements ..
307*
308* Test input arguments
309*
310 info = 0
311 wantu1 = lsame( jobu1, 'Y' )
312 wantu2 = lsame( jobu2, 'Y' )
313 wantv1t = lsame( jobv1t, 'Y' )
314 lquery = ( lwork.EQ.-1 ) .OR. ( lrwork.EQ.-1 )
315*
316 IF( m .LT. 0 ) THEN
317 info = -4
318 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
319 info = -5
320 ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
321 info = -6
322 ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
323 info = -8
324 ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
325 info = -10
326 ELSE IF( wantu1 .AND. ldu1 .LT. max( 1, p ) ) THEN
327 info = -13
328 ELSE IF( wantu2 .AND. ldu2 .LT. max( 1, m - p ) ) THEN
329 info = -15
330 ELSE IF( wantv1t .AND. ldv1t .LT. max( 1, q ) ) THEN
331 info = -17
332 END IF
333*
334 r = min( p, m-p, q, m-q )
335*
336* Compute workspace
337*
338* WORK layout:
339* |-----------------------------------------|
340* | LWORKOPT (1) |
341* |-----------------------------------------|
342* | TAUP1 (MAX(1,P)) |
343* | TAUP2 (MAX(1,M-P)) |
344* | TAUQ1 (MAX(1,Q)) |
345* |-----------------------------------------|
346* | ZUNBDB WORK | ZUNGQR WORK | ZUNGLQ WORK |
347* | | | |
348* | | | |
349* | | | |
350* | | | |
351* |-----------------------------------------|
352* RWORK layout:
353* |------------------|
354* | LRWORKOPT (1) |
355* |------------------|
356* | PHI (MAX(1,R-1)) |
357* |------------------|
358* | B11D (R) |
359* | B11E (R-1) |
360* | B12D (R) |
361* | B12E (R-1) |
362* | B21D (R) |
363* | B21E (R-1) |
364* | B22D (R) |
365* | B22E (R-1) |
366* | ZBBCSD RWORK |
367* |------------------|
368*
369 IF( info .EQ. 0 ) THEN
370 iphi = 2
371 ib11d = iphi + max( 1, r-1 )
372 ib11e = ib11d + max( 1, r )
373 ib12d = ib11e + max( 1, r - 1 )
374 ib12e = ib12d + max( 1, r )
375 ib21d = ib12e + max( 1, r - 1 )
376 ib21e = ib21d + max( 1, r )
377 ib22d = ib21e + max( 1, r - 1 )
378 ib22e = ib22d + max( 1, r )
379 ibbcsd = ib22e + max( 1, r - 1 )
380 itaup1 = 2
381 itaup2 = itaup1 + max( 1, p )
382 itauq1 = itaup2 + max( 1, m-p )
383 iorbdb = itauq1 + max( 1, q )
384 iorgqr = itauq1 + max( 1, q )
385 iorglq = itauq1 + max( 1, q )
386 lorgqrmin = 1
387 lorgqropt = 1
388 lorglqmin = 1
389 lorglqopt = 1
390 IF( r .EQ. q ) THEN
391 CALL zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
392 $ dum,
393 $ cdum, cdum, cdum, work, -1, childinfo )
394 lorbdb = int( work(1) )
395 IF( wantu1 .AND. p .GT. 0 ) THEN
396 CALL zungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
397 $ childinfo )
398 lorgqrmin = max( lorgqrmin, p )
399 lorgqropt = max( lorgqropt, int( work(1) ) )
400 ENDIF
401 IF( wantu2 .AND. m-p .GT. 0 ) THEN
402 CALL zungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
403 $ childinfo )
404 lorgqrmin = max( lorgqrmin, m-p )
405 lorgqropt = max( lorgqropt, int( work(1) ) )
406 END IF
407 IF( wantv1t .AND. q .GT. 0 ) THEN
408 CALL zunglq( q-1, q-1, q-1, v1t, ldv1t,
409 $ cdum, work(1), -1, childinfo )
410 lorglqmin = max( lorglqmin, q-1 )
411 lorglqopt = max( lorglqopt, int( work(1) ) )
412 END IF
413 CALL zbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q,
414 $ theta,
415 $ dum, u1, ldu1, u2, ldu2, v1t, ldv1t, cdum, 1,
416 $ dum, dum, dum, dum, dum, dum, dum, dum,
417 $ rwork(1), -1, childinfo )
418 lbbcsd = int( rwork(1) )
419 ELSE IF( r .EQ. p ) THEN
420 CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
421 $ dum,
422 $ cdum, cdum, cdum, work(1), -1, childinfo )
423 lorbdb = int( work(1) )
424 IF( wantu1 .AND. p .GT. 0 ) THEN
425 CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, cdum,
426 $ work(1),
427 $ -1, childinfo )
428 lorgqrmin = max( lorgqrmin, p-1 )
429 lorgqropt = max( lorgqropt, int( work(1) ) )
430 END IF
431 IF( wantu2 .AND. m-p .GT. 0 ) THEN
432 CALL zungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
433 $ childinfo )
434 lorgqrmin = max( lorgqrmin, m-p )
435 lorgqropt = max( lorgqropt, int( work(1) ) )
436 END IF
437 IF( wantv1t .AND. q .GT. 0 ) THEN
438 CALL zunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
439 $ childinfo )
440 lorglqmin = max( lorglqmin, q )
441 lorglqopt = max( lorglqopt, int( work(1) ) )
442 END IF
443 CALL zbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p,
444 $ theta,
445 $ dum, v1t, ldv1t, cdum, 1, u1, ldu1, u2, ldu2,
446 $ dum, dum, dum, dum, dum, dum, dum, dum,
447 $ rwork(1), -1, childinfo )
448 lbbcsd = int( rwork(1) )
449 ELSE IF( r .EQ. m-p ) THEN
450 CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
451 $ dum,
452 $ cdum, cdum, cdum, work(1), -1, childinfo )
453 lorbdb = int( work(1) )
454 IF( wantu1 .AND. p .GT. 0 ) THEN
455 CALL zungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
456 $ childinfo )
457 lorgqrmin = max( lorgqrmin, p )
458 lorgqropt = max( lorgqropt, int( work(1) ) )
459 END IF
460 IF( wantu2 .AND. m-p .GT. 0 ) THEN
461 CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, cdum,
462 $ work(1), -1, childinfo )
463 lorgqrmin = max( lorgqrmin, m-p-1 )
464 lorgqropt = max( lorgqropt, int( work(1) ) )
465 END IF
466 IF( wantv1t .AND. q .GT. 0 ) THEN
467 CALL zunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
468 $ childinfo )
469 lorglqmin = max( lorglqmin, q )
470 lorglqopt = max( lorglqopt, int( work(1) ) )
471 END IF
472 CALL zbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
473 $ theta, dum, cdum, 1, v1t, ldv1t, u2, ldu2, u1,
474 $ ldu1, dum, dum, dum, dum, dum, dum, dum, dum,
475 $ rwork(1), -1, childinfo )
476 lbbcsd = int( rwork(1) )
477 ELSE
478 CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
479 $ dum,
480 $ cdum, cdum, cdum, cdum, work(1), -1, childinfo
481 $ )
482 lorbdb = m + int( work(1) )
483 IF( wantu1 .AND. p .GT. 0 ) THEN
484 CALL zungqr( p, p, m-q, u1, ldu1, cdum, work(1), -1,
485 $ childinfo )
486 lorgqrmin = max( lorgqrmin, p )
487 lorgqropt = max( lorgqropt, int( work(1) ) )
488 END IF
489 IF( wantu2 .AND. m-p .GT. 0 ) THEN
490 CALL zungqr( m-p, m-p, m-q, u2, ldu2, cdum, work(1),
491 $ -1,
492 $ childinfo )
493 lorgqrmin = max( lorgqrmin, m-p )
494 lorgqropt = max( lorgqropt, int( work(1) ) )
495 END IF
496 IF( wantv1t .AND. q .GT. 0 ) THEN
497 CALL zunglq( q, q, q, v1t, ldv1t, cdum, work(1), -1,
498 $ childinfo )
499 lorglqmin = max( lorglqmin, q )
500 lorglqopt = max( lorglqopt, int( work(1) ) )
501 END IF
502 CALL zbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
503 $ theta, dum, u2, ldu2, u1, ldu1, cdum, 1, v1t,
504 $ ldv1t, dum, dum, dum, dum, dum, dum, dum, dum,
505 $ rwork(1), -1, childinfo )
506 lbbcsd = int( rwork(1) )
507 END IF
508 lrworkmin = ibbcsd+lbbcsd-1
509 lrworkopt = lrworkmin
510 rwork(1) = lrworkopt
511 lworkmin = max( iorbdb+lorbdb-1,
512 $ iorgqr+lorgqrmin-1,
513 $ iorglq+lorglqmin-1 )
514 lworkopt = max( iorbdb+lorbdb-1,
515 $ iorgqr+lorgqropt-1,
516 $ iorglq+lorglqopt-1 )
517 work(1) = lworkopt
518 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
519 info = -19
520 END IF
521 IF( lrwork .LT. lrworkmin .AND. .NOT.lquery ) THEN
522 info = -21
523 END IF
524 END IF
525 IF( info .NE. 0 ) THEN
526 CALL xerbla( 'ZUNCSD2BY1', -info )
527 RETURN
528 ELSE IF( lquery ) THEN
529 RETURN
530 END IF
531 lorgqr = lwork-iorgqr+1
532 lorglq = lwork-iorglq+1
533*
534* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
535* in which R = MIN(P,M-P,Q,M-Q)
536*
537 IF( r .EQ. q ) THEN
538*
539* Case 1: R = Q
540*
541* Simultaneously bidiagonalize X11 and X21
542*
543 CALL zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
544 $ rwork(iphi), work(itaup1), work(itaup2),
545 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
546*
547* Accumulate Householder reflectors
548*
549 IF( wantu1 .AND. p .GT. 0 ) THEN
550 CALL zlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
551 CALL zungqr( p, p, q, u1, ldu1, work(itaup1),
552 $ work(iorgqr),
553 $ lorgqr, childinfo )
554 END IF
555 IF( wantu2 .AND. m-p .GT. 0 ) THEN
556 CALL zlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
557 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
558 $ work(iorgqr), lorgqr, childinfo )
559 END IF
560 IF( wantv1t .AND. q .GT. 0 ) THEN
561 v1t(1,1) = one
562 DO j = 2, q
563 v1t(1,j) = zero
564 v1t(j,1) = zero
565 END DO
566 CALL zlacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
567 $ ldv1t )
568 CALL zunglq( q-1, q-1, q-1, v1t(2,2), ldv1t,
569 $ work(itauq1),
570 $ work(iorglq), lorglq, childinfo )
571 END IF
572*
573* Simultaneously diagonalize X11 and X21.
574*
575 CALL zbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
576 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
577 $ 1, rwork(ib11d), rwork(ib11e), rwork(ib12d),
578 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
579 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
580 $ lrwork-ibbcsd+1, childinfo )
581*
582* Permute rows and columns to place zero submatrices in
583* preferred positions
584*
585 IF( q .GT. 0 .AND. wantu2 ) THEN
586 DO i = 1, q
587 iwork(i) = m - p - q + i
588 END DO
589 DO i = q + 1, m - p
590 iwork(i) = i - q
591 END DO
592 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
593 END IF
594 ELSE IF( r .EQ. p ) THEN
595*
596* Case 2: R = P
597*
598* Simultaneously bidiagonalize X11 and X21
599*
600 CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
601 $ rwork(iphi), work(itaup1), work(itaup2),
602 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
603*
604* Accumulate Householder reflectors
605*
606 IF( wantu1 .AND. p .GT. 0 ) THEN
607 u1(1,1) = one
608 DO j = 2, p
609 u1(1,j) = zero
610 u1(j,1) = zero
611 END DO
612 CALL zlacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2),
613 $ ldu1 )
614 CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
615 $ work(iorgqr), lorgqr, childinfo )
616 END IF
617 IF( wantu2 .AND. m-p .GT. 0 ) THEN
618 CALL zlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
619 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
620 $ work(iorgqr), lorgqr, childinfo )
621 END IF
622 IF( wantv1t .AND. q .GT. 0 ) THEN
623 CALL zlacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
624 CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
625 $ work(iorglq), lorglq, childinfo )
626 END IF
627*
628* Simultaneously diagonalize X11 and X21.
629*
630 CALL zbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
631 $ rwork(iphi), v1t, ldv1t, cdum, 1, u1, ldu1, u2,
632 $ ldu2, rwork(ib11d), rwork(ib11e), rwork(ib12d),
633 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
634 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
635 $ childinfo )
636*
637* Permute rows and columns to place identity submatrices in
638* preferred positions
639*
640 IF( q .GT. 0 .AND. wantu2 ) THEN
641 DO i = 1, q
642 iwork(i) = m - p - q + i
643 END DO
644 DO i = q + 1, m - p
645 iwork(i) = i - q
646 END DO
647 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
648 END IF
649 ELSE IF( r .EQ. m-p ) THEN
650*
651* Case 3: R = M-P
652*
653* Simultaneously bidiagonalize X11 and X21
654*
655 CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
656 $ rwork(iphi), work(itaup1), work(itaup2),
657 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
658*
659* Accumulate Householder reflectors
660*
661 IF( wantu1 .AND. p .GT. 0 ) THEN
662 CALL zlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
663 CALL zungqr( p, p, q, u1, ldu1, work(itaup1),
664 $ work(iorgqr),
665 $ lorgqr, childinfo )
666 END IF
667 IF( wantu2 .AND. m-p .GT. 0 ) THEN
668 u2(1,1) = one
669 DO j = 2, m-p
670 u2(1,j) = zero
671 u2(j,1) = zero
672 END DO
673 CALL zlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
674 $ ldu2 )
675 CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
676 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
677 END IF
678 IF( wantv1t .AND. q .GT. 0 ) THEN
679 CALL zlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
680 CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
681 $ work(iorglq), lorglq, childinfo )
682 END IF
683*
684* Simultaneously diagonalize X11 and X21.
685*
686 CALL zbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
687 $ theta, rwork(iphi), cdum, 1, v1t, ldv1t, u2, ldu2,
688 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
689 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
690 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
691 $ rwork(ibbcsd), lbbcsd, childinfo )
692*
693* Permute rows and columns to place identity submatrices in
694* preferred positions
695*
696 IF( q .GT. r ) THEN
697 DO i = 1, r
698 iwork(i) = q - r + i
699 END DO
700 DO i = r + 1, q
701 iwork(i) = i - r
702 END DO
703 IF( wantu1 ) THEN
704 CALL zlapmt( .false., p, q, u1, ldu1, iwork )
705 END IF
706 IF( wantv1t ) THEN
707 CALL zlapmr( .false., q, q, v1t, ldv1t, iwork )
708 END IF
709 END IF
710 ELSE
711*
712* Case 4: R = M-Q
713*
714* Simultaneously bidiagonalize X11 and X21
715*
716 CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
717 $ rwork(iphi), work(itaup1), work(itaup2),
718 $ work(itauq1), work(iorbdb), work(iorbdb+m),
719 $ lorbdb-m, childinfo )
720*
721* Accumulate Householder reflectors
722*
723 IF( wantu2 .AND. m-p .GT. 0 ) THEN
724 CALL zcopy( m-p, work(iorbdb+p), 1, u2, 1 )
725 END IF
726 IF( wantu1 .AND. p .GT. 0 ) THEN
727 CALL zcopy( p, work(iorbdb), 1, u1, 1 )
728 DO j = 2, p
729 u1(1,j) = zero
730 END DO
731 CALL zlacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
732 $ ldu1 )
733 CALL zungqr( p, p, m-q, u1, ldu1, work(itaup1),
734 $ work(iorgqr), lorgqr, childinfo )
735 END IF
736 IF( wantu2 .AND. m-p .GT. 0 ) THEN
737 DO j = 2, m-p
738 u2(1,j) = zero
739 END DO
740 CALL zlacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
741 $ ldu2 )
742 CALL zungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
743 $ work(iorgqr), lorgqr, childinfo )
744 END IF
745 IF( wantv1t .AND. q .GT. 0 ) THEN
746 CALL zlacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
747 CALL zlacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1),
748 $ ldx11,
749 $ v1t(m-q+1,m-q+1), ldv1t )
750 CALL zlacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
751 $ v1t(p+1,p+1), ldv1t )
752 CALL zunglq( q, q, q, v1t, ldv1t, work(itauq1),
753 $ work(iorglq), lorglq, childinfo )
754 END IF
755*
756* Simultaneously diagonalize X11 and X21.
757*
758 CALL zbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
759 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, cdum, 1,
760 $ v1t, ldv1t, rwork(ib11d), rwork(ib11e),
761 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
762 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
763 $ rwork(ibbcsd), lbbcsd, childinfo )
764*
765* Permute rows and columns to place identity submatrices in
766* preferred positions
767*
768 IF( p .GT. r ) THEN
769 DO i = 1, r
770 iwork(i) = p - r + i
771 END DO
772 DO i = r + 1, p
773 iwork(i) = i - r
774 END DO
775 IF( wantu1 ) THEN
776 CALL zlapmt( .false., p, p, u1, ldu1, iwork )
777 END IF
778 IF( wantv1t ) THEN
779 CALL zlapmr( .false., p, q, v1t, ldv1t, iwork )
780 END IF
781 END IF
782 END IF
783*
784 RETURN
785*
786* End of ZUNCSD2BY1
787*
788 END
789
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zbbcsd
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:331
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
zlapmr
subroutine zlapmr(forwrd, m, n, x, ldx, k)
ZLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition zlapmr.f:102
zlapmt
subroutine zlapmt(forwrd, m, n, x, ldx, k)
ZLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition zlapmt.f:102
zunbdb1
subroutine zunbdb1(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
ZUNBDB1
Definition zunbdb1.f:202
zunbdb2
subroutine zunbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
ZUNBDB2
Definition zunbdb2.f:200
zunbdb3
subroutine zunbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
ZUNBDB3
Definition zunbdb3.f:200
zunbdb4
subroutine zunbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
ZUNBDB4
Definition zunbdb4.f:212
zuncsd2by1
subroutine zuncsd2by1(jobu1, jobu2, jobv1t, m, p, q, x11, ldx11, x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work, lwork, rwork, lrwork, iwork, info)
ZUNCSD2BY1
Definition zuncsd2by1.f:255
zunglq
subroutine zunglq(m, n, k, a, lda, tau, work, lwork, info)
ZUNGLQ
Definition zunglq.f:125
zungqr
subroutine zungqr(m, n, k, a, lda, tau, work, lwork, info)
ZUNGQR
Definition zungqr.f:126