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