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