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