LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
cuncsd.f
Go to the documentation of this file.
1 *> \brief \b CUNCSD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CUNCSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cuncsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cuncsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cuncsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * RECURSIVE SUBROUTINE CUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
22 * SIGNS, M, P, Q, X11, LDX11, X12,
23 * LDX12, X21, LDX21, X22, LDX22, THETA,
24 * U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
25 * LDV2T, WORK, LWORK, RWORK, LRWORK,
26 * IWORK, INFO )
27 *
28 * .. Scalar Arguments ..
29 * CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
30 * INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
31 * $ LDX21, LDX22, LRWORK, LWORK, M, P, Q
32 * ..
33 * .. Array Arguments ..
34 * INTEGER IWORK( * )
35 * REAL THETA( * )
36 * REAL RWORK( * )
37 * COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
38 * $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
39 * $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
40 * $ * )
41 * ..
42 *
43 *
44 *> \par Purpose:
45 * =============
46 *>
47 *> \verbatim
48 *>
49 *> CUNCSD computes the CS decomposition of an M-by-M partitioned
50 *> unitary matrix X:
51 *>
52 *> [ I 0 0 | 0 0 0 ]
53 *> [ 0 C 0 | 0 -S 0 ]
54 *> [ X11 | X12 ] [ U1 | ] [ 0 0 0 | 0 0 -I ] [ V1 | ]**H
55 *> X = [-----------] = [---------] [---------------------] [---------] .
56 *> [ X21 | X22 ] [ | U2 ] [ 0 0 0 | I 0 0 ] [ | V2 ]
57 *> [ 0 S 0 | 0 C 0 ]
58 *> [ 0 0 I | 0 0 0 ]
59 *>
60 *> X11 is P-by-Q. The unitary matrices U1, U2, V1, and V2 are P-by-P,
61 *> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
62 *> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
63 *> which R = MIN(P,M-P,Q,M-Q).
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] JOBV2T
91 *> \verbatim
92 *> JOBV2T is CHARACTER
93 *> = 'Y': V2T is computed;
94 *> otherwise: V2T is not computed.
95 *> \endverbatim
96 *>
97 *> \param[in] TRANS
98 *> \verbatim
99 *> TRANS is CHARACTER
100 *> = 'T': X, U1, U2, V1T, and V2T are stored in row-major
101 *> order;
102 *> otherwise: X, U1, U2, V1T, and V2T are stored in column-
103 *> major order.
104 *> \endverbatim
105 *>
106 *> \param[in] SIGNS
107 *> \verbatim
108 *> SIGNS is CHARACTER
109 *> = 'O': The lower-left block is made nonpositive (the
110 *> "other" convention);
111 *> otherwise: The upper-right block is made nonpositive (the
112 *> "default" convention).
113 *> \endverbatim
114 *>
115 *> \param[in] M
116 *> \verbatim
117 *> M is INTEGER
118 *> The number of rows and columns in X.
119 *> \endverbatim
120 *>
121 *> \param[in] P
122 *> \verbatim
123 *> P is INTEGER
124 *> The number of rows in X11 and X12. 0 <= P <= M.
125 *> \endverbatim
126 *>
127 *> \param[in] Q
128 *> \verbatim
129 *> Q is INTEGER
130 *> The number of columns in X11 and X21. 0 <= Q <= M.
131 *> \endverbatim
132 *>
133 *> \param[in,out] X11
134 *> \verbatim
135 *> X11 is COMPLEX array, dimension (LDX11,Q)
136 *> On entry, part of the unitary matrix whose CSD is desired.
137 *> \endverbatim
138 *>
139 *> \param[in] LDX11
140 *> \verbatim
141 *> LDX11 is INTEGER
142 *> The leading dimension of X11. LDX11 >= MAX(1,P).
143 *> \endverbatim
144 *>
145 *> \param[in,out] X12
146 *> \verbatim
147 *> X12 is COMPLEX array, dimension (LDX12,M-Q)
148 *> On entry, part of the unitary matrix whose CSD is desired.
149 *> \endverbatim
150 *>
151 *> \param[in] LDX12
152 *> \verbatim
153 *> LDX12 is INTEGER
154 *> The leading dimension of X12. LDX12 >= MAX(1,P).
155 *> \endverbatim
156 *>
157 *> \param[in,out] X21
158 *> \verbatim
159 *> X21 is COMPLEX array, dimension (LDX21,Q)
160 *> On entry, part of the unitary matrix whose CSD is desired.
161 *> \endverbatim
162 *>
163 *> \param[in] LDX21
164 *> \verbatim
165 *> LDX21 is INTEGER
166 *> The leading dimension of X11. LDX21 >= MAX(1,M-P).
167 *> \endverbatim
168 *>
169 *> \param[in,out] X22
170 *> \verbatim
171 *> X22 is COMPLEX array, dimension (LDX22,M-Q)
172 *> On entry, part of the unitary matrix whose CSD is desired.
173 *> \endverbatim
174 *>
175 *> \param[in] LDX22
176 *> \verbatim
177 *> LDX22 is INTEGER
178 *> The leading dimension of X11. LDX22 >= MAX(1,M-P).
179 *> \endverbatim
180 *>
181 *> \param[out] THETA
182 *> \verbatim
183 *> THETA is REAL array, dimension (R), in which R =
184 *> MIN(P,M-P,Q,M-Q).
185 *> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
186 *> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
187 *> \endverbatim
188 *>
189 *> \param[out] U1
190 *> \verbatim
191 *> U1 is COMPLEX array, dimension (LDU1,P)
192 *> If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1.
193 *> \endverbatim
194 *>
195 *> \param[in] LDU1
196 *> \verbatim
197 *> LDU1 is INTEGER
198 *> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
199 *> MAX(1,P).
200 *> \endverbatim
201 *>
202 *> \param[out] U2
203 *> \verbatim
204 *> U2 is COMPLEX array, dimension (LDU2,M-P)
205 *> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
206 *> matrix U2.
207 *> \endverbatim
208 *>
209 *> \param[in] LDU2
210 *> \verbatim
211 *> LDU2 is INTEGER
212 *> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
213 *> MAX(1,M-P).
214 *> \endverbatim
215 *>
216 *> \param[out] V1T
217 *> \verbatim
218 *> V1T is COMPLEX array, dimension (LDV1T,Q)
219 *> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
220 *> matrix V1**H.
221 *> \endverbatim
222 *>
223 *> \param[in] LDV1T
224 *> \verbatim
225 *> LDV1T is INTEGER
226 *> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
227 *> MAX(1,Q).
228 *> \endverbatim
229 *>
230 *> \param[out] V2T
231 *> \verbatim
232 *> V2T is COMPLEX array, dimension (LDV2T,M-Q)
233 *> If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) unitary
234 *> matrix V2**H.
235 *> \endverbatim
236 *>
237 *> \param[in] LDV2T
238 *> \verbatim
239 *> LDV2T is INTEGER
240 *> The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
241 *> MAX(1,M-Q).
242 *> \endverbatim
243 *>
244 *> \param[out] WORK
245 *> \verbatim
246 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
247 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
248 *> \endverbatim
249 *>
250 *> \param[in] LWORK
251 *> \verbatim
252 *> LWORK is INTEGER
253 *> The dimension of the array WORK.
254 *>
255 *> If LWORK = -1, then a workspace query is assumed; the routine
256 *> only calculates the optimal size of the WORK array, returns
257 *> this value as the first entry of the work array, and no error
258 *> message related to LWORK is issued by XERBLA.
259 *> \endverbatim
260 *>
261 *> \param[out] RWORK
262 *> \verbatim
263 *> RWORK is REAL array, dimension MAX(1,LRWORK)
264 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
265 *> If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
266 *> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
267 *> define the matrix in intermediate bidiagonal-block form
268 *> remaining after nonconvergence. INFO specifies the number
269 *> of nonzero PHI's.
270 *> \endverbatim
271 *>
272 *> \param[in] LRWORK
273 *> \verbatim
274 *> LRWORK is INTEGER
275 *> The dimension of the array RWORK.
276 *>
277 *> If LRWORK = -1, then a workspace query is assumed; the routine
278 *> only calculates the optimal size of the RWORK array, returns
279 *> this value as the first entry of the work array, and no error
280 *> message related to LRWORK is issued by XERBLA.
281 *> \endverbatim
282 *>
283 *> \param[out] IWORK
284 *> \verbatim
285 *> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
286 *> \endverbatim
287 *>
288 *> \param[out] INFO
289 *> \verbatim
290 *> INFO is INTEGER
291 *> = 0: successful exit.
292 *> < 0: if INFO = -i, the i-th argument had an illegal value.
293 *> > 0: CBBCSD did not converge. See the description of RWORK
294 *> above for details.
295 *> \endverbatim
296 *
297 *> \par References:
298 * ================
299 *>
300 *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
301 *> Algorithms, 50(1):33-65, 2009.
302 *
303 * Authors:
304 * ========
305 *
306 *> \author Univ. of Tennessee
307 *> \author Univ. of California Berkeley
308 *> \author Univ. of Colorado Denver
309 *> \author NAG Ltd.
310 *
311 *> \ingroup complexOTHERcomputational
312 *
313 * =====================================================================
314  RECURSIVE SUBROUTINE cuncsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
315  $ SIGNS, M, P, Q, X11, LDX11, X12,
316  $ LDX12, X21, LDX21, X22, LDX22, THETA,
317  $ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
318  $ LDV2T, WORK, LWORK, RWORK, LRWORK,
319  $ IWORK, INFO )
320 *
321 * -- LAPACK computational routine --
322 * -- LAPACK is a software package provided by Univ. of Tennessee, --
323 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
324 *
325 * .. Scalar Arguments ..
326  CHARACTER jobu1, jobu2, jobv1t, jobv2t, signs, trans
327  INTEGER info, ldu1, ldu2, ldv1t, ldv2t, ldx11, ldx12,
328  $ ldx21, ldx22, lrwork, lwork, m, p, q
329 * ..
330 * .. Array Arguments ..
331  INTEGER iwork( * )
332  REAL theta( * )
333  REAL rwork( * )
334  COMPLEX u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
335  $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
336  $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
337  $ * )
338 * ..
339 *
340 * ===================================================================
341 *
342 * .. Parameters ..
343  COMPLEX one, zero
344  PARAMETER ( one = (1.0e0,0.0e0),
345  $ zero = (0.0e0,0.0e0) )
346 * ..
347 * .. Local Scalars ..
348  CHARACTER transt, signst
349  INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
350  $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
351  $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
352  $ itauq2, j, lbbcsdwork, lbbcsdworkmin,
353  $ lbbcsdworkopt, lorbdbwork, lorbdbworkmin,
354  $ lorbdbworkopt, lorglqwork, lorglqworkmin,
355  $ lorglqworkopt, lorgqrwork, lorgqrworkmin,
356  $ lorgqrworkopt, lworkmin, lworkopt, p1, q1
357  LOGICAL colmajor, defaultsigns, lquery, wantu1, wantu2,
358  $ wantv1t, wantv2t
359  INTEGER lrworkmin, lrworkopt
360  LOGICAL lrquery
361 * ..
362 * .. External Subroutines ..
363  EXTERNAL xerbla, cbbcsd, clacpy, clapmr, clapmt,
364  $ cunbdb, cunglq, cungqr
365 * ..
366 * .. External Functions ..
367  LOGICAL lsame
368  EXTERNAL lsame
369 * ..
370 * .. Intrinsic Functions
371  INTRINSIC int, max, min
372 * ..
373 * .. Executable Statements ..
374 *
375 * Test input arguments
376 *
377  info = 0
378  wantu1 = lsame( jobu1, 'Y' )
379  wantu2 = lsame( jobu2, 'Y' )
380  wantv1t = lsame( jobv1t, 'Y' )
381  wantv2t = lsame( jobv2t, 'Y' )
382  colmajor = .NOT. lsame( trans, 'T' )
383  defaultsigns = .NOT. lsame( signs, 'O' )
384  lquery = lwork .EQ. -1
385  lrquery = lrwork .EQ. -1
386  IF( m .LT. 0 ) THEN
387  info = -7
388  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
389  info = -8
390  ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
391  info = -9
392  ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
393  info = -11
394  ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
395  info = -11
396  ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
397  info = -13
398  ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
399  info = -13
400  ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
401  info = -15
402  ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
403  info = -15
404  ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
405  info = -17
406  ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
407  info = -17
408  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
409  info = -20
410  ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
411  info = -22
412  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
413  info = -24
414  ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
415  info = -26
416  END IF
417 *
418 * Work with transpose if convenient
419 *
420  IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) ) THEN
421  IF( colmajor ) THEN
422  transt = 'T'
423  ELSE
424  transt = 'N'
425  END IF
426  IF( defaultsigns ) THEN
427  signst = 'O'
428  ELSE
429  signst = 'D'
430  END IF
431  CALL cuncsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
432  $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
433  $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
434  $ u2, ldu2, work, lwork, rwork, lrwork, iwork,
435  $ info )
436  RETURN
437  END IF
438 *
439 * Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
440 * convenient
441 *
442  IF( info .EQ. 0 .AND. m-q .LT. q ) THEN
443  IF( defaultsigns ) THEN
444  signst = 'O'
445  ELSE
446  signst = 'D'
447  END IF
448  CALL cuncsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
449  $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
450  $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
451  $ ldv1t, work, lwork, rwork, lrwork, iwork, info )
452  RETURN
453  END IF
454 *
455 * Compute workspace
456 *
457  IF( info .EQ. 0 ) THEN
458 *
459 * Real workspace
460 *
461  iphi = 2
462  ib11d = iphi + max( 1, q - 1 )
463  ib11e = ib11d + max( 1, q )
464  ib12d = ib11e + max( 1, q - 1 )
465  ib12e = ib12d + max( 1, q )
466  ib21d = ib12e + max( 1, q - 1 )
467  ib21e = ib21d + max( 1, q )
468  ib22d = ib21e + max( 1, q - 1 )
469  ib22e = ib22d + max( 1, q )
470  ibbcsd = ib22e + max( 1, q - 1 )
471  CALL cbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
472  $ theta, theta, u1, ldu1, u2, ldu2, v1t, ldv1t,
473  $ v2t, ldv2t, theta, theta, theta, theta, theta,
474  $ theta, theta, theta, rwork, -1, childinfo )
475  lbbcsdworkopt = int( rwork(1) )
476  lbbcsdworkmin = lbbcsdworkopt
477  lrworkopt = ibbcsd + lbbcsdworkopt - 1
478  lrworkmin = ibbcsd + lbbcsdworkmin - 1
479  rwork(1) = lrworkopt
480 *
481 * Complex workspace
482 *
483  itaup1 = 2
484  itaup2 = itaup1 + max( 1, p )
485  itauq1 = itaup2 + max( 1, m - p )
486  itauq2 = itauq1 + max( 1, q )
487  iorgqr = itauq2 + max( 1, m - q )
488  CALL cungqr( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
489  $ childinfo )
490  lorgqrworkopt = int( work(1) )
491  lorgqrworkmin = max( 1, m - q )
492  iorglq = itauq2 + max( 1, m - q )
493  CALL cunglq( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
494  $ childinfo )
495  lorglqworkopt = int( work(1) )
496  lorglqworkmin = max( 1, m - q )
497  iorbdb = itauq2 + max( 1, m - q )
498  CALL cunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
499  $ x21, ldx21, x22, ldx22, theta, theta, u1, u2,
500  $ v1t, v2t, work, -1, childinfo )
501  lorbdbworkopt = int( work(1) )
502  lorbdbworkmin = lorbdbworkopt
503  lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
504  $ iorbdb + lorbdbworkopt ) - 1
505  lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
506  $ iorbdb + lorbdbworkmin ) - 1
507  work(1) = max(lworkopt,lworkmin)
508 *
509  IF( lwork .LT. lworkmin
510  $ .AND. .NOT. ( lquery .OR. lrquery ) ) THEN
511  info = -22
512  ELSE IF( lrwork .LT. lrworkmin
513  $ .AND. .NOT. ( lquery .OR. lrquery ) ) THEN
514  info = -24
515  ELSE
516  lorgqrwork = lwork - iorgqr + 1
517  lorglqwork = lwork - iorglq + 1
518  lorbdbwork = lwork - iorbdb + 1
519  lbbcsdwork = lrwork - ibbcsd + 1
520  END IF
521  END IF
522 *
523 * Abort if any illegal arguments
524 *
525  IF( info .NE. 0 ) THEN
526  CALL xerbla( 'CUNCSD', -info )
527  RETURN
528  ELSE IF( lquery .OR. lrquery ) THEN
529  RETURN
530  END IF
531 *
532 * Transform to bidiagonal block form
533 *
534  CALL cunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
535  $ ldx21, x22, ldx22, theta, rwork(iphi), work(itaup1),
536  $ work(itaup2), work(itauq1), work(itauq2),
537  $ work(iorbdb), lorbdbwork, childinfo )
538 *
539 * Accumulate Householder reflectors
540 *
541  IF( colmajor ) THEN
542  IF( wantu1 .AND. p .GT. 0 ) THEN
543  CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
544  CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
545  $ lorgqrwork, info)
546  END IF
547  IF( wantu2 .AND. m-p .GT. 0 ) THEN
548  CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
549  CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
550  $ work(iorgqr), lorgqrwork, info )
551  END IF
552  IF( wantv1t .AND. q .GT. 0 ) THEN
553  CALL clacpy( 'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
554  $ ldv1t )
555  v1t(1, 1) = one
556  DO j = 2, q
557  v1t(1,j) = zero
558  v1t(j,1) = zero
559  END DO
560  CALL cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
561  $ work(iorglq), lorglqwork, info )
562  END IF
563  IF( wantv2t .AND. m-q .GT. 0 ) THEN
564  CALL clacpy( 'U', p, m-q, x12, ldx12, v2t, ldv2t )
565  IF( m-p .GT. q ) THEN
566  CALL clacpy( 'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
567  $ v2t(p+1,p+1), ldv2t )
568  END IF
569  IF( m .GT. q ) THEN
570  CALL cunglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
571  $ work(iorglq), lorglqwork, info )
572  END IF
573  END IF
574  ELSE
575  IF( wantu1 .AND. p .GT. 0 ) THEN
576  CALL clacpy( 'U', q, p, x11, ldx11, u1, ldu1 )
577  CALL cunglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
578  $ lorglqwork, info)
579  END IF
580  IF( wantu2 .AND. m-p .GT. 0 ) THEN
581  CALL clacpy( 'U', q, m-p, x21, ldx21, u2, ldu2 )
582  CALL cunglq( m-p, m-p, q, u2, ldu2, work(itaup2),
583  $ work(iorglq), lorglqwork, info )
584  END IF
585  IF( wantv1t .AND. q .GT. 0 ) THEN
586  CALL clacpy( 'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
587  $ ldv1t )
588  v1t(1, 1) = one
589  DO j = 2, q
590  v1t(1,j) = zero
591  v1t(j,1) = zero
592  END DO
593  CALL cungqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
594  $ work(iorgqr), lorgqrwork, info )
595  END IF
596  IF( wantv2t .AND. m-q .GT. 0 ) THEN
597  p1 = min( p+1, m )
598  q1 = min( q+1, m )
599  CALL clacpy( 'L', m-q, p, x12, ldx12, v2t, ldv2t )
600  IF ( m .GT. p+q ) THEN
601  CALL clacpy( 'L', m-p-q, m-p-q, x22(p1,q1), ldx22,
602  $ v2t(p+1,p+1), ldv2t )
603  END IF
604  CALL cungqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
605  $ work(iorgqr), lorgqrwork, info )
606  END IF
607  END IF
608 *
609 * Compute the CSD of the matrix in bidiagonal-block form
610 *
611  CALL cbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
612  $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
613  $ ldv2t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
614  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
615  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
616  $ lbbcsdwork, info )
617 *
618 * Permute rows and columns to place identity submatrices in top-
619 * left corner of (1,1)-block and/or bottom-right corner of (1,2)-
620 * block and/or bottom-right corner of (2,1)-block and/or top-left
621 * corner of (2,2)-block
622 *
623  IF( q .GT. 0 .AND. wantu2 ) THEN
624  DO i = 1, q
625  iwork(i) = m - p - q + i
626  END DO
627  DO i = q + 1, m - p
628  iwork(i) = i - q
629  END DO
630  IF( colmajor ) THEN
631  CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
632  ELSE
633  CALL clapmr( .false., m-p, m-p, u2, ldu2, iwork )
634  END IF
635  END IF
636  IF( m .GT. 0 .AND. wantv2t ) THEN
637  DO i = 1, p
638  iwork(i) = m - p - q + i
639  END DO
640  DO i = p + 1, m - q
641  iwork(i) = i - p
642  END DO
643  IF( .NOT. colmajor ) THEN
644  CALL clapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
645  ELSE
646  CALL clapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
647  END IF
648  END IF
649 *
650  RETURN
651 *
652 * End CUNCSD
653 *
654  END
655 
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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 cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
Definition: cunglq.f:127
recursive subroutine cuncsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
CUNCSD
Definition: cuncsd.f:320
subroutine cunbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
CUNBDB
Definition: cunbdb.f:287
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 cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
Definition: cungqr.f:128