LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
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 (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 REAL 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 REAL 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 REAL 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 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 June 2017
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.7.1) --
304 * -- LAPACK is a software package provided by Univ. of Tennessee, --
305 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
306 * June 2017
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,
347  $ sorbdb, sorglq, sorgqr, xerbla
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 
sorglq
subroutine sorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGLQ
Definition: sorglq.f:129
sorbdb
subroutine sorbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
SORBDB
Definition: sorbdb.f:289
sorcsd
recursive subroutine sorcsd(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)
SORCSD
Definition: sorcsd.f:302
slacpy
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
sorgqr
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
Definition: sorgqr.f:130
slapmr
subroutine slapmr(FORWRD, M, N, X, LDX, K)
SLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: slapmr.f:106
slapmt
subroutine slapmt(FORWRD, M, N, X, LDX, K)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: slapmt.f:106
sbbcsd
subroutine sbbcsd(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)
SBBCSD
Definition: sbbcsd.f:334