LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Files Functions 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 
subroutine cuncsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
CUNCSD2BY1
Definition: cuncsd2by1.f:260
subroutine clapmr(FORWRD, M, N, X, LDX, K)
CLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: clapmr.f:105
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
subroutine cunbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB1
Definition: cunbdb1.f:202
subroutine cbbcsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E, B22D, B22E, RWORK, LRWORK, INFO)
CBBCSD
Definition: cbbcsd.f:330
subroutine cunbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
CUNBDB4
Definition: cunbdb4.f:212
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:104
subroutine clapmt(FORWRD, M, N, X, LDX, K)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: clapmt.f:105
subroutine cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
Definition: cunglq.f:128
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:51
subroutine cunbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB2
Definition: cunbdb2.f:202
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
Definition: cungqr.f:129
subroutine cunbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB3
Definition: cunbdb3.f:202