LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
sorcsd.f
Go to the documentation of this file.
1 *> \brief \b SORCSD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SORCSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorcsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorcsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorcsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * RECURSIVE SUBROUTINE SORCSD( 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 * REAL THETA( * )
35 * REAL 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 *> SORCSD 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (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 REAL array, dimension (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 REAL array, dimension (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 REAL array, dimension (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 REAL 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: SBBCSD 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 *> \date November 2013
293 *
294 *> \ingroup realOTHERcomputational
295 *
296 * =====================================================================
297  RECURSIVE SUBROUTINE sorcsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
298  $ signs, m, p, q, x11, ldx11, x12,
299  $ ldx12, x21, ldx21, x22, ldx22, theta,
300  $ u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
301  $ ldv2t, work, lwork, iwork, info )
302 *
303 * -- LAPACK computational routine (version 3.5.0) --
304 * -- LAPACK is a software package provided by Univ. of Tennessee, --
305 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
306 * November 2013
307 *
308 * .. Scalar Arguments ..
309  CHARACTER jobu1, jobu2, jobv1t, jobv2t, signs, trans
310  INTEGER info, ldu1, ldu2, ldv1t, ldv2t, ldx11, ldx12,
311  $ ldx21, ldx22, lwork, m, p, q
312 * ..
313 * .. Array Arguments ..
314  INTEGER iwork( * )
315  REAL theta( * )
316  REAL u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
317  $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
318  $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
319  $ * )
320 * ..
321 *
322 * ===================================================================
323 *
324 * .. Parameters ..
325  REAL one, zero
326  parameter( one = 1.0e+0,
327  $ zero = 0.0e+0 )
328 * ..
329 * .. Local Arrays ..
330  REAL dummy(1)
331 * ..
332 * .. Local Scalars ..
333  CHARACTER transt, signst
334  INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
335  $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
336  $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
337  $ itauq2, j, lbbcsdwork, lbbcsdworkmin,
338  $ lbbcsdworkopt, lorbdbwork, lorbdbworkmin,
339  $ lorbdbworkopt, lorglqwork, lorglqworkmin,
340  $ lorglqworkopt, lorgqrwork, lorgqrworkmin,
341  $ lorgqrworkopt, lworkmin, lworkopt
342  LOGICAL colmajor, defaultsigns, lquery, wantu1, wantu2,
343  $ wantv1t, wantv2t
344 * ..
345 * .. External Subroutines ..
346  EXTERNAL sbbcsd, slacpy, slapmr, slapmt, slascl, slaset,
348 * ..
349 * .. External Functions ..
350  LOGICAL lsame
351  EXTERNAL lsame
352 * ..
353 * .. Intrinsic Functions
354  INTRINSIC int, max, min
355 * ..
356 * .. Executable Statements ..
357 *
358 * Test input arguments
359 *
360  info = 0
361  wantu1 = lsame( jobu1, 'Y' )
362  wantu2 = lsame( jobu2, 'Y' )
363  wantv1t = lsame( jobv1t, 'Y' )
364  wantv2t = lsame( jobv2t, 'Y' )
365  colmajor = .NOT. lsame( trans, 'T' )
366  defaultsigns = .NOT. lsame( signs, 'O' )
367  lquery = lwork .EQ. -1
368  IF( m .LT. 0 ) THEN
369  info = -7
370  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
371  info = -8
372  ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
373  info = -9
374  ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
375  info = -11
376  ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
377  info = -11
378  ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
379  info = -13
380  ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
381  info = -13
382  ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
383  info = -15
384  ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
385  info = -15
386  ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
387  info = -17
388  ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
389  info = -17
390  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
391  info = -20
392  ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
393  info = -22
394  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
395  info = -24
396  ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
397  info = -26
398  END IF
399 *
400 * Work with transpose if convenient
401 *
402  IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) ) THEN
403  IF( colmajor ) THEN
404  transt = 'T'
405  ELSE
406  transt = 'N'
407  END IF
408  IF( defaultsigns ) THEN
409  signst = 'O'
410  ELSE
411  signst = 'D'
412  END IF
413  CALL sorcsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
414  $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
415  $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
416  $ u2, ldu2, work, lwork, iwork, info )
417  RETURN
418  END IF
419 *
420 * Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
421 * convenient
422 *
423  IF( info .EQ. 0 .AND. m-q .LT. q ) THEN
424  IF( defaultsigns ) THEN
425  signst = 'O'
426  ELSE
427  signst = 'D'
428  END IF
429  CALL sorcsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
430  $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
431  $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
432  $ ldv1t, work, lwork, iwork, info )
433  RETURN
434  END IF
435 *
436 * Compute workspace
437 *
438  IF( info .EQ. 0 ) THEN
439 *
440  iphi = 2
441  itaup1 = iphi + max( 1, q - 1 )
442  itaup2 = itaup1 + max( 1, p )
443  itauq1 = itaup2 + max( 1, m - p )
444  itauq2 = itauq1 + max( 1, q )
445  iorgqr = itauq2 + max( 1, m - q )
446  CALL sorgqr( m-q, m-q, m-q, dummy, max(1,m-q), dummy, work, -1,
447  $ childinfo )
448  lorgqrworkopt = int( work(1) )
449  lorgqrworkmin = max( 1, m - q )
450  iorglq = itauq2 + max( 1, m - q )
451  CALL sorglq( m-q, m-q, m-q, dummy, max(1,m-q), dummy, work, -1,
452  $ childinfo )
453  lorglqworkopt = int( work(1) )
454  lorglqworkmin = max( 1, m - q )
455  iorbdb = itauq2 + max( 1, m - q )
456  CALL sorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
457  $ x21, ldx21, x22, ldx22, dummy, dummy, dummy, dummy, dummy,
458  $ dummy,work,-1,childinfo )
459  lorbdbworkopt = int( work(1) )
460  lorbdbworkmin = lorbdbworkopt
461  ib11d = itauq2 + max( 1, m - q )
462  ib11e = ib11d + max( 1, q )
463  ib12d = ib11e + max( 1, q - 1 )
464  ib12e = ib12d + max( 1, q )
465  ib21d = ib12e + max( 1, q - 1 )
466  ib21e = ib21d + max( 1, q )
467  ib22d = ib21e + max( 1, q - 1 )
468  ib22e = ib22d + max( 1, q )
469  ibbcsd = ib22e + max( 1, q - 1 )
470  CALL sbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
471  $ dummy, dummy, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
472  $ ldv2t, dummy, dummy, dummy, dummy, dummy, dummy,
473  $ dummy, dummy, work, -1, childinfo )
474  lbbcsdworkopt = int( work(1) )
475  lbbcsdworkmin = lbbcsdworkopt
476  lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
477  $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkopt ) - 1
478  lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
479  $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkmin ) - 1
480  work(1) = max(lworkopt,lworkmin)
481 *
482  IF( lwork .LT. lworkmin .AND. .NOT. lquery ) THEN
483  info = -22
484  ELSE
485  lorgqrwork = lwork - iorgqr + 1
486  lorglqwork = lwork - iorglq + 1
487  lorbdbwork = lwork - iorbdb + 1
488  lbbcsdwork = lwork - ibbcsd + 1
489  END IF
490  END IF
491 *
492 * Abort if any illegal arguments
493 *
494  IF( info .NE. 0 ) THEN
495  CALL xerbla( 'SORCSD', -info )
496  RETURN
497  ELSE IF( lquery ) THEN
498  RETURN
499  END IF
500 *
501 * Transform to bidiagonal block form
502 *
503  CALL sorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
504  $ ldx21, x22, ldx22, theta, work(iphi), work(itaup1),
505  $ work(itaup2), work(itauq1), work(itauq2),
506  $ work(iorbdb), lorbdbwork, childinfo )
507 *
508 * Accumulate Householder reflectors
509 *
510  IF( colmajor ) THEN
511  IF( wantu1 .AND. p .GT. 0 ) THEN
512  CALL slacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
513  CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
514  $ lorgqrwork, info)
515  END IF
516  IF( wantu2 .AND. m-p .GT. 0 ) THEN
517  CALL slacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
518  CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
519  $ work(iorgqr), lorgqrwork, info )
520  END IF
521  IF( wantv1t .AND. q .GT. 0 ) THEN
522  CALL slacpy( 'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
523  $ ldv1t )
524  v1t(1, 1) = one
525  DO j = 2, q
526  v1t(1,j) = zero
527  v1t(j,1) = zero
528  END DO
529  CALL sorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
530  $ work(iorglq), lorglqwork, info )
531  END IF
532  IF( wantv2t .AND. m-q .GT. 0 ) THEN
533  CALL slacpy( 'U', p, m-q, x12, ldx12, v2t, ldv2t )
534  CALL slacpy( 'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
535  $ v2t(p+1,p+1), ldv2t )
536  CALL sorglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
537  $ work(iorglq), lorglqwork, info )
538  END IF
539  ELSE
540  IF( wantu1 .AND. p .GT. 0 ) THEN
541  CALL slacpy( 'U', q, p, x11, ldx11, u1, ldu1 )
542  CALL sorglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
543  $ lorglqwork, info)
544  END IF
545  IF( wantu2 .AND. m-p .GT. 0 ) THEN
546  CALL slacpy( 'U', q, m-p, x21, ldx21, u2, ldu2 )
547  CALL sorglq( m-p, m-p, q, u2, ldu2, work(itaup2),
548  $ work(iorglq), lorglqwork, info )
549  END IF
550  IF( wantv1t .AND. q .GT. 0 ) THEN
551  CALL slacpy( 'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
552  $ ldv1t )
553  v1t(1, 1) = one
554  DO j = 2, q
555  v1t(1,j) = zero
556  v1t(j,1) = zero
557  END DO
558  CALL sorgqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
559  $ work(iorgqr), lorgqrwork, info )
560  END IF
561  IF( wantv2t .AND. m-q .GT. 0 ) THEN
562  CALL slacpy( 'L', m-q, p, x12, ldx12, v2t, ldv2t )
563  CALL slacpy( 'L', m-p-q, m-p-q, x22(p+1,q+1), ldx22,
564  $ v2t(p+1,p+1), ldv2t )
565  CALL sorgqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
566  $ work(iorgqr), lorgqrwork, info )
567  END IF
568  END IF
569 *
570 * Compute the CSD of the matrix in bidiagonal-block form
571 *
572  CALL sbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
573  $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
574  $ ldv2t, work(ib11d), work(ib11e), work(ib12d),
575  $ work(ib12e), work(ib21d), work(ib21e), work(ib22d),
576  $ work(ib22e), work(ibbcsd), lbbcsdwork, info )
577 *
578 * Permute rows and columns to place identity submatrices in top-
579 * left corner of (1,1)-block and/or bottom-right corner of (1,2)-
580 * block and/or bottom-right corner of (2,1)-block and/or top-left
581 * corner of (2,2)-block
582 *
583  IF( q .GT. 0 .AND. wantu2 ) THEN
584  DO i = 1, q
585  iwork(i) = m - p - q + i
586  END DO
587  DO i = q + 1, m - p
588  iwork(i) = i - q
589  END DO
590  IF( colmajor ) THEN
591  CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
592  ELSE
593  CALL slapmr( .false., m-p, m-p, u2, ldu2, iwork )
594  END IF
595  END IF
596  IF( m .GT. 0 .AND. wantv2t ) THEN
597  DO i = 1, p
598  iwork(i) = m - p - q + i
599  END DO
600  DO i = p + 1, m - q
601  iwork(i) = i - p
602  END DO
603  IF( .NOT. colmajor ) THEN
604  CALL slapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
605  ELSE
606  CALL slapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
607  END IF
608  END IF
609 *
610  RETURN
611 *
612 * End SORCSD
613 *
614  END
615