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