LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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