LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
cuncsd2by1.f
Go to the documentation of this file.
1 *> \brief \b CUNCSD2BY1
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CUNCSD2BY1 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cuncsd2by1.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cuncsd2by1.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cuncsd2by1.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CUNCSD2BY1( 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 * REAL RWORK(*)
34 * REAL THETA(*)
35 * COMPLEX 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 *> CUNCSD2BY1 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 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 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 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 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 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 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 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 REAL 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: CBBCSD 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 *
247 * Authors:
248 * ========
249 *
250 *> \author Univ. of Tennessee
251 *> \author Univ. of California Berkeley
252 *> \author Univ. of Colorado Denver
253 *> \author NAG Ltd.
254 *
255 *> \date July 2012
256 *
257 *> \ingroup complexOTHERcomputational
258 *
259 * =====================================================================
260  SUBROUTINE cuncsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
261  $ x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t,
262  $ ldv1t, work, lwork, rwork, lrwork, iwork,
263  $ info )
264 *
265 * -- LAPACK computational routine (version 3.5.0) --
266 * -- LAPACK is a software package provided by Univ. of Tennessee, --
267 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268 * July 2012
269 *
270 * .. Scalar Arguments ..
271  CHARACTER jobu1, jobu2, jobv1t
272  INTEGER info, ldu1, ldu2, ldv1t, lwork, ldx11, ldx21,
273  $ m, p, q
274  INTEGER lrwork, lrworkmin, lrworkopt
275 * ..
276 * .. Array Arguments ..
277  REAL rwork(*)
278  REAL theta(*)
279  COMPLEX u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), work(*),
280  $ x11(ldx11,*), x21(ldx21,*)
281  INTEGER iwork(*)
282 * ..
283 *
284 * =====================================================================
285 *
286 * .. Parameters ..
287  COMPLEX one, zero
288  parameter( one = (1.0e0,0.0e0), zero = (0.0e0,0.0e0) )
289 * ..
290 * .. Local Scalars ..
291  INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
292  $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
293  $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
294  $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
295  $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
296  $ lworkmin, lworkopt, r
297  LOGICAL lquery, wantu1, wantu2, wantv1t
298 * ..
299 * .. External Subroutines ..
300  EXTERNAL cbbcsd, ccopy, clacpy, clapmr, clapmt, cunbdb1,
302  $ xerbla
303 * ..
304 * .. External Functions ..
305  LOGICAL lsame
306  EXTERNAL lsame
307 * ..
308 * .. Intrinsic Function ..
309  INTRINSIC int, max, min
310 * ..
311 * .. Executable Statements ..
312 *
313 * Test input arguments
314 *
315  info = 0
316  wantu1 = lsame( jobu1, 'Y' )
317  wantu2 = lsame( jobu2, 'Y' )
318  wantv1t = lsame( jobv1t, 'Y' )
319  lquery = lwork .EQ. -1
320 *
321  IF( m .LT. 0 ) THEN
322  info = -4
323  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
324  info = -5
325  ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
326  info = -6
327  ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
328  info = -8
329  ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
330  info = -10
331  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
332  info = -13
333  ELSE IF( wantu2 .AND. ldu2 .LT. m - p ) THEN
334  info = -15
335  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
336  info = -17
337  END IF
338 *
339  r = min( p, m-p, q, m-q )
340 *
341 * Compute workspace
342 *
343 * WORK layout:
344 * |-----------------------------------------|
345 * | LWORKOPT (1) |
346 * |-----------------------------------------|
347 * | TAUP1 (MAX(1,P)) |
348 * | TAUP2 (MAX(1,M-P)) |
349 * | TAUQ1 (MAX(1,Q)) |
350 * |-----------------------------------------|
351 * | CUNBDB WORK | CUNGQR WORK | CUNGLQ WORK |
352 * | | | |
353 * | | | |
354 * | | | |
355 * | | | |
356 * |-----------------------------------------|
357 * RWORK layout:
358 * |------------------|
359 * | LRWORKOPT (1) |
360 * |------------------|
361 * | PHI (MAX(1,R-1)) |
362 * |------------------|
363 * | B11D (R) |
364 * | B11E (R-1) |
365 * | B12D (R) |
366 * | B12E (R-1) |
367 * | B21D (R) |
368 * | B21E (R-1) |
369 * | B22D (R) |
370 * | B22E (R-1) |
371 * | CBBCSD RWORK |
372 * |------------------|
373 *
374  IF( info .EQ. 0 ) THEN
375  iphi = 2
376  ib11d = iphi + max( 1, r-1 )
377  ib11e = ib11d + max( 1, r )
378  ib12d = ib11e + max( 1, r - 1 )
379  ib12e = ib12d + max( 1, r )
380  ib21d = ib12e + max( 1, r - 1 )
381  ib21e = ib21d + max( 1, r )
382  ib22d = ib21e + max( 1, r - 1 )
383  ib22e = ib22d + max( 1, r )
384  ibbcsd = ib22e + max( 1, r - 1 )
385  itaup1 = 2
386  itaup2 = itaup1 + max( 1, p )
387  itauq1 = itaup2 + max( 1, m-p )
388  iorbdb = itauq1 + max( 1, q )
389  iorgqr = itauq1 + max( 1, q )
390  iorglq = itauq1 + max( 1, q )
391  IF( r .EQ. q ) THEN
392  CALL cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
393  $ 0, 0, work, -1, childinfo )
394  lorbdb = int( work(1) )
395  IF( p .GE. m-p ) THEN
396  CALL cungqr( p, p, q, u1, ldu1, 0, work(1), -1,
397  $ childinfo )
398  lorgqrmin = max( 1, p )
399  lorgqropt = int( work(1) )
400  ELSE
401  CALL cungqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
402  $ childinfo )
403  lorgqrmin = max( 1, m-p )
404  lorgqropt = int( work(1) )
405  END IF
406  CALL cunglq( max(0,q-1), max(0,q-1), max(0,q-1), v1t, ldv1t,
407  $ 0, work(1), -1, childinfo )
408  lorglqmin = max( 1, q-1 )
409  lorglqopt = int( work(1) )
410  CALL cbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
411  $ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1, 0, 0,
412  $ 0, 0, 0, 0, 0, 0, rwork(1), -1, childinfo )
413  lbbcsd = int( rwork(1) )
414  ELSE IF( r .EQ. p ) THEN
415  CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
416  $ 0, 0, work(1), -1, childinfo )
417  lorbdb = int( work(1) )
418  IF( p-1 .GE. m-p ) THEN
419  CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, 0, work(1),
420  $ -1, childinfo )
421  lorgqrmin = max( 1, p-1 )
422  lorgqropt = int( work(1) )
423  ELSE
424  CALL cungqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
425  $ childinfo )
426  lorgqrmin = max( 1, m-p )
427  lorgqropt = int( work(1) )
428  END IF
429  CALL cunglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
430  $ childinfo )
431  lorglqmin = max( 1, q )
432  lorglqopt = int( work(1) )
433  CALL cbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
434  $ 0, v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2, 0, 0,
435  $ 0, 0, 0, 0, 0, 0, rwork(1), -1, childinfo )
436  lbbcsd = int( rwork(1) )
437  ELSE IF( r .EQ. m-p ) THEN
438  CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
439  $ 0, 0, work(1), -1, childinfo )
440  lorbdb = int( work(1) )
441  IF( p .GE. m-p-1 ) THEN
442  CALL cungqr( p, p, q, u1, ldu1, 0, work(1), -1,
443  $ childinfo )
444  lorgqrmin = max( 1, p )
445  lorgqropt = int( work(1) )
446  ELSE
447  CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, 0,
448  $ work(1), -1, childinfo )
449  lorgqrmin = max( 1, m-p-1 )
450  lorgqropt = int( work(1) )
451  END IF
452  CALL cunglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
453  $ childinfo )
454  lorglqmin = max( 1, q )
455  lorglqopt = int( work(1) )
456  CALL cbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
457  $ theta, 0, 0, 1, v1t, ldv1t, u2, ldu2, u1, ldu1,
458  $ 0, 0, 0, 0, 0, 0, 0, 0, rwork(1), -1,
459  $ childinfo )
460  lbbcsd = int( rwork(1) )
461  ELSE
462  CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
463  $ 0, 0, 0, work(1), -1, childinfo )
464  lorbdb = m + int( work(1) )
465  IF( p .GE. m-p ) THEN
466  CALL cungqr( p, p, m-q, u1, ldu1, 0, work(1), -1,
467  $ childinfo )
468  lorgqrmin = max( 1, p )
469  lorgqropt = int( work(1) )
470  ELSE
471  CALL cungqr( m-p, m-p, m-q, u2, ldu2, 0, work(1), -1,
472  $ childinfo )
473  lorgqrmin = max( 1, m-p )
474  lorgqropt = int( work(1) )
475  END IF
476  CALL cunglq( q, q, q, v1t, ldv1t, 0, work(1), -1,
477  $ childinfo )
478  lorglqmin = max( 1, q )
479  lorglqopt = int( work(1) )
480  CALL cbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
481  $ theta, 0, u2, ldu2, u1, ldu1, 0, 1, v1t, ldv1t,
482  $ 0, 0, 0, 0, 0, 0, 0, 0, rwork(1), -1,
483  $ childinfo )
484  lbbcsd = int( rwork(1) )
485  END IF
486  lrworkmin = ibbcsd+lbbcsd-1
487  lrworkopt = lrworkmin
488  rwork(1) = lrworkopt
489  lworkmin = max( iorbdb+lorbdb-1,
490  $ iorgqr+lorgqrmin-1,
491  $ iorglq+lorglqmin-1 )
492  lworkopt = max( iorbdb+lorbdb-1,
493  $ iorgqr+lorgqropt-1,
494  $ iorglq+lorglqopt-1 )
495  work(1) = lworkopt
496  IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
497  info = -19
498  END IF
499  END IF
500  IF( info .NE. 0 ) THEN
501  CALL xerbla( 'CUNCSD2BY1', -info )
502  RETURN
503  ELSE IF( lquery ) THEN
504  RETURN
505  END IF
506  lorgqr = lwork-iorgqr+1
507  lorglq = lwork-iorglq+1
508 *
509 * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
510 * in which R = MIN(P,M-P,Q,M-Q)
511 *
512  IF( r .EQ. q ) THEN
513 *
514 * Case 1: R = Q
515 *
516 * Simultaneously bidiagonalize X11 and X21
517 *
518  CALL cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
519  $ rwork(iphi), work(itaup1), work(itaup2),
520  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
521 *
522 * Accumulate Householder reflectors
523 *
524  IF( wantu1 .AND. p .GT. 0 ) THEN
525  CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
526  CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
527  $ lorgqr, childinfo )
528  END IF
529  IF( wantu2 .AND. m-p .GT. 0 ) THEN
530  CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
531  CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
532  $ work(iorgqr), lorgqr, childinfo )
533  END IF
534  IF( wantv1t .AND. q .GT. 0 ) THEN
535  v1t(1,1) = one
536  DO j = 2, q
537  v1t(1,j) = zero
538  v1t(j,1) = zero
539  END DO
540  CALL clacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
541  $ ldv1t )
542  CALL cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
543  $ work(iorglq), lorglq, childinfo )
544  END IF
545 *
546 * Simultaneously diagonalize X11 and X21.
547 *
548  CALL cbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
549  $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1,
550  $ rwork(ib11d), rwork(ib11e), rwork(ib12d),
551  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
552  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
553  $ childinfo )
554 *
555 * Permute rows and columns to place zero submatrices in
556 * preferred positions
557 *
558  IF( q .GT. 0 .AND. wantu2 ) THEN
559  DO i = 1, q
560  iwork(i) = m - p - q + i
561  END DO
562  DO i = q + 1, m - p
563  iwork(i) = i - q
564  END DO
565  CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
566  END IF
567  ELSE IF( r .EQ. p ) THEN
568 *
569 * Case 2: R = P
570 *
571 * Simultaneously bidiagonalize X11 and X21
572 *
573  CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
574  $ rwork(iphi), work(itaup1), work(itaup2),
575  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
576 *
577 * Accumulate Householder reflectors
578 *
579  IF( wantu1 .AND. p .GT. 0 ) THEN
580  u1(1,1) = one
581  DO j = 2, p
582  u1(1,j) = zero
583  u1(j,1) = zero
584  END DO
585  CALL clacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
586  CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
587  $ work(iorgqr), lorgqr, childinfo )
588  END IF
589  IF( wantu2 .AND. m-p .GT. 0 ) THEN
590  CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
591  CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
592  $ work(iorgqr), lorgqr, childinfo )
593  END IF
594  IF( wantv1t .AND. q .GT. 0 ) THEN
595  CALL clacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
596  CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
597  $ work(iorglq), lorglq, childinfo )
598  END IF
599 *
600 * Simultaneously diagonalize X11 and X21.
601 *
602  CALL cbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
603  $ rwork(iphi), v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2,
604  $ rwork(ib11d), rwork(ib11e), rwork(ib12d),
605  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
606  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
607  $ childinfo )
608 *
609 * Permute rows and columns to place identity submatrices in
610 * preferred positions
611 *
612  IF( q .GT. 0 .AND. wantu2 ) THEN
613  DO i = 1, q
614  iwork(i) = m - p - q + i
615  END DO
616  DO i = q + 1, m - p
617  iwork(i) = i - q
618  END DO
619  CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
620  END IF
621  ELSE IF( r .EQ. m-p ) THEN
622 *
623 * Case 3: R = M-P
624 *
625 * Simultaneously bidiagonalize X11 and X21
626 *
627  CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
628  $ rwork(iphi), work(itaup1), work(itaup2),
629  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
630 *
631 * Accumulate Householder reflectors
632 *
633  IF( wantu1 .AND. p .GT. 0 ) THEN
634  CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
635  CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
636  $ lorgqr, childinfo )
637  END IF
638  IF( wantu2 .AND. m-p .GT. 0 ) THEN
639  u2(1,1) = one
640  DO j = 2, m-p
641  u2(1,j) = zero
642  u2(j,1) = zero
643  END DO
644  CALL clacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
645  $ ldu2 )
646  CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
647  $ work(itaup2), work(iorgqr), lorgqr, childinfo )
648  END IF
649  IF( wantv1t .AND. q .GT. 0 ) THEN
650  CALL clacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
651  CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
652  $ work(iorglq), lorglq, childinfo )
653  END IF
654 *
655 * Simultaneously diagonalize X11 and X21.
656 *
657  CALL cbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
658  $ theta, rwork(iphi), 0, 1, v1t, ldv1t, u2, ldu2,
659  $ u1, ldu1, rwork(ib11d), rwork(ib11e),
660  $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
661  $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
662  $ rwork(ibbcsd), lbbcsd, childinfo )
663 *
664 * Permute rows and columns to place identity submatrices in
665 * preferred positions
666 *
667  IF( q .GT. r ) THEN
668  DO i = 1, r
669  iwork(i) = q - r + i
670  END DO
671  DO i = r + 1, q
672  iwork(i) = i - r
673  END DO
674  IF( wantu1 ) THEN
675  CALL clapmt( .false., p, q, u1, ldu1, iwork )
676  END IF
677  IF( wantv1t ) THEN
678  CALL clapmr( .false., q, q, v1t, ldv1t, iwork )
679  END IF
680  END IF
681  ELSE
682 *
683 * Case 4: R = M-Q
684 *
685 * Simultaneously bidiagonalize X11 and X21
686 *
687  CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
688  $ rwork(iphi), work(itaup1), work(itaup2),
689  $ work(itauq1), work(iorbdb), work(iorbdb+m),
690  $ lorbdb-m, childinfo )
691 *
692 * Accumulate Householder reflectors
693 *
694  IF( wantu1 .AND. p .GT. 0 ) THEN
695  CALL ccopy( p, work(iorbdb), 1, u1, 1 )
696  DO j = 2, p
697  u1(1,j) = zero
698  END DO
699  CALL clacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
700  $ ldu1 )
701  CALL cungqr( p, p, m-q, u1, ldu1, work(itaup1),
702  $ work(iorgqr), lorgqr, childinfo )
703  END IF
704  IF( wantu2 .AND. m-p .GT. 0 ) THEN
705  CALL ccopy( m-p, work(iorbdb+p), 1, u2, 1 )
706  DO j = 2, m-p
707  u2(1,j) = zero
708  END DO
709  CALL clacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
710  $ ldu2 )
711  CALL cungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
712  $ work(iorgqr), lorgqr, childinfo )
713  END IF
714  IF( wantv1t .AND. q .GT. 0 ) THEN
715  CALL clacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
716  CALL clacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
717  $ v1t(m-q+1,m-q+1), ldv1t )
718  CALL clacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
719  $ v1t(p+1,p+1), ldv1t )
720  CALL cunglq( q, q, q, v1t, ldv1t, work(itauq1),
721  $ work(iorglq), lorglq, childinfo )
722  END IF
723 *
724 * Simultaneously diagonalize X11 and X21.
725 *
726  CALL cbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
727  $ theta, rwork(iphi), u2, ldu2, u1, ldu1, 0, 1, v1t,
728  $ ldv1t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
729  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
730  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
731  $ childinfo )
732 *
733 * Permute rows and columns to place identity submatrices in
734 * preferred positions
735 *
736  IF( p .GT. r ) THEN
737  DO i = 1, r
738  iwork(i) = p - r + i
739  END DO
740  DO i = r + 1, p
741  iwork(i) = i - r
742  END DO
743  IF( wantu1 ) THEN
744  CALL clapmt( .false., p, p, u1, ldu1, iwork )
745  END IF
746  IF( wantv1t ) THEN
747  CALL clapmr( .false., p, q, v1t, ldv1t, iwork )
748  END IF
749  END IF
750  END IF
751 *
752  RETURN
753 *
754 * End of CUNCSD2BY1
755 *
756  END
757