LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dorcsd2by1.f
Go to the documentation of this file.
1*> \brief \b DORCSD2BY1
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DORCSD2BY1 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd2by1.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd2by1.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd2by1.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
20* X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
21* LDV1T, WORK, LWORK, IWORK, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER JOBU1, JOBU2, JOBV1T
25* INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
26* $ M, P, Q
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION THETA(*)
30* DOUBLE PRECISION U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
31* $ X11(LDX11,*), X21(LDX21,*)
32* INTEGER IWORK(*)
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*>\verbatim
40*>
41*> DORCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with
42*> orthonormal columns that has been partitioned into a 2-by-1 block
43*> structure:
44*>
45*> [ I1 0 0 ]
46*> [ 0 C 0 ]
47*> [ X11 ] [ U1 | ] [ 0 0 0 ]
48*> X = [-----] = [---------] [----------] V1**T .
49*> [ X21 ] [ | U2 ] [ 0 0 0 ]
50*> [ 0 S 0 ]
51*> [ 0 0 I2]
52*>
53*> X11 is P-by-Q. The orthogonal matrices U1, U2, and V1 are P-by-P,
54*> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R
55*> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which
56*> R = MIN(P,M-P,Q,M-Q). I1 is a K1-by-K1 identity matrix and I2 is a
57*> K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0).
58*> \endverbatim
59*
60* Arguments:
61* ==========
62*
63*> \param[in] JOBU1
64*> \verbatim
65*> JOBU1 is CHARACTER
66*> = 'Y': U1 is computed;
67*> otherwise: U1 is not computed.
68*> \endverbatim
69*>
70*> \param[in] JOBU2
71*> \verbatim
72*> JOBU2 is CHARACTER
73*> = 'Y': U2 is computed;
74*> otherwise: U2 is not computed.
75*> \endverbatim
76*>
77*> \param[in] JOBV1T
78*> \verbatim
79*> JOBV1T is CHARACTER
80*> = 'Y': V1T is computed;
81*> otherwise: V1T is not computed.
82*> \endverbatim
83*>
84*> \param[in] M
85*> \verbatim
86*> M is INTEGER
87*> The number of rows in X.
88*> \endverbatim
89*>
90*> \param[in] P
91*> \verbatim
92*> P is INTEGER
93*> The number of rows in X11. 0 <= P <= M.
94*> \endverbatim
95*>
96*> \param[in] Q
97*> \verbatim
98*> Q is INTEGER
99*> The number of columns in X11 and X21. 0 <= Q <= M.
100*> \endverbatim
101*>
102*> \param[in,out] X11
103*> \verbatim
104*> X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
105*> On entry, part of the orthogonal matrix whose CSD is desired.
106*> \endverbatim
107*>
108*> \param[in] LDX11
109*> \verbatim
110*> LDX11 is INTEGER
111*> The leading dimension of X11. LDX11 >= MAX(1,P).
112*> \endverbatim
113*>
114*> \param[in,out] X21
115*> \verbatim
116*> X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
117*> On entry, part of the orthogonal matrix whose CSD is desired.
118*> \endverbatim
119*>
120*> \param[in] LDX21
121*> \verbatim
122*> LDX21 is INTEGER
123*> The leading dimension of X21. LDX21 >= MAX(1,M-P).
124*> \endverbatim
125*>
126*> \param[out] THETA
127*> \verbatim
128*> THETA is DOUBLE PRECISION array, dimension (R), in which R =
129*> MIN(P,M-P,Q,M-Q).
130*> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
131*> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
132*> \endverbatim
133*>
134*> \param[out] U1
135*> \verbatim
136*> U1 is DOUBLE PRECISION array, dimension (P)
137*> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
138*> \endverbatim
139*>
140*> \param[in] LDU1
141*> \verbatim
142*> LDU1 is INTEGER
143*> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
144*> MAX(1,P).
145*> \endverbatim
146*>
147*> \param[out] U2
148*> \verbatim
149*> U2 is DOUBLE PRECISION array, dimension (M-P)
150*> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
151*> matrix U2.
152*> \endverbatim
153*>
154*> \param[in] LDU2
155*> \verbatim
156*> LDU2 is INTEGER
157*> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
158*> MAX(1,M-P).
159*> \endverbatim
160*>
161*> \param[out] V1T
162*> \verbatim
163*> V1T is DOUBLE PRECISION array, dimension (Q)
164*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
165*> matrix V1**T.
166*> \endverbatim
167*>
168*> \param[in] LDV1T
169*> \verbatim
170*> LDV1T is INTEGER
171*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
172*> MAX(1,Q).
173*> \endverbatim
174*>
175*> \param[out] WORK
176*> \verbatim
177*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
178*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
179*> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
180*> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
181*> define the matrix in intermediate bidiagonal-block form
182*> remaining after nonconvergence. INFO specifies the number
183*> of nonzero PHI's.
184*> \endverbatim
185*>
186*> \param[in] LWORK
187*> \verbatim
188*> LWORK is INTEGER
189*> The dimension of the array WORK.
190*>
191*> If LWORK = -1, then a workspace query is assumed; the routine
192*> only calculates the optimal size of the WORK array, returns
193*> this value as the first entry of the work array, and no error
194*> message related to LWORK is issued by XERBLA.
195*> \endverbatim
196*>
197*> \param[out] IWORK
198*> \verbatim
199*> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
200*> \endverbatim
201*>
202*> \param[out] INFO
203*> \verbatim
204*> INFO is INTEGER
205*> = 0: successful exit.
206*> < 0: if INFO = -i, the i-th argument had an illegal value.
207*> > 0: DBBCSD did not converge. See the description of WORK
208*> above for details.
209*> \endverbatim
210*
211*> \par References:
212* ================
213*>
214*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
215*> Algorithms, 50(1):33-65, 2009.
216*
217* Authors:
218* ========
219*
220*> \author Univ. of Tennessee
221*> \author Univ. of California Berkeley
222*> \author Univ. of Colorado Denver
223*> \author NAG Ltd.
224*
225*> \ingroup uncsd2by1
226*
227* =====================================================================
228 SUBROUTINE dorcsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11,
229 $ LDX11,
230 $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
231 $ LDV1T, WORK, LWORK, IWORK, INFO )
232*
233* -- LAPACK computational routine (3.5.0) --
234* -- LAPACK is a software package provided by Univ. of Tennessee, --
235* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
236*
237* .. Scalar Arguments ..
238 CHARACTER JOBU1, JOBU2, JOBV1T
239 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
240 $ M, P, Q
241* ..
242* .. Array Arguments ..
243 DOUBLE PRECISION THETA(*)
244 DOUBLE PRECISION U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
245 $ X11(LDX11,*), X21(LDX21,*)
246 INTEGER IWORK(*)
247* ..
248*
249* =====================================================================
250*
251* .. Parameters ..
252 DOUBLE PRECISION ONE, ZERO
253 PARAMETER ( ONE = 1.0d0, zero = 0.0d0 )
254* ..
255* .. Local Scalars ..
256 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
257 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
258 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
259 $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
260 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
261 $ lworkmin, lworkopt, r
262 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
263* ..
264* .. Local Arrays ..
265 DOUBLE PRECISION DUM1(1), DUM2(1,1)
266* ..
267* .. External Subroutines ..
268 EXTERNAL dbbcsd, dcopy, dlacpy, dlapmr, dlapmt,
269 $ dorbdb1,
270 $ dorbdb2, dorbdb3, dorbdb4, dorglq, dorgqr,
271 $ xerbla
272* ..
273* .. External Functions ..
274 LOGICAL LSAME
275 EXTERNAL LSAME
276* ..
277* .. Intrinsic Function ..
278 INTRINSIC int, max, min
279* ..
280* .. Executable Statements ..
281*
282* Test input arguments
283*
284 info = 0
285 wantu1 = lsame( jobu1, 'Y' )
286 wantu2 = lsame( jobu2, 'Y' )
287 wantv1t = lsame( jobv1t, 'Y' )
288 lquery = lwork .EQ. -1
289*
290 IF( m .LT. 0 ) THEN
291 info = -4
292 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
293 info = -5
294 ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
295 info = -6
296 ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
297 info = -8
298 ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
299 info = -10
300 ELSE IF( wantu1 .AND. ldu1 .LT. max( 1, p ) ) THEN
301 info = -13
302 ELSE IF( wantu2 .AND. ldu2 .LT. max( 1, m - p ) ) THEN
303 info = -15
304 ELSE IF( wantv1t .AND. ldv1t .LT. max( 1, q ) ) THEN
305 info = -17
306 END IF
307*
308 r = min( p, m-p, q, m-q )
309*
310* Compute workspace
311*
312* WORK layout:
313* |-------------------------------------------------------|
314* | LWORKOPT (1) |
315* |-------------------------------------------------------|
316* | PHI (MAX(1,R-1)) |
317* |-------------------------------------------------------|
318* | TAUP1 (MAX(1,P)) | B11D (R) |
319* | TAUP2 (MAX(1,M-P)) | B11E (R-1) |
320* | TAUQ1 (MAX(1,Q)) | B12D (R) |
321* |-----------------------------------------| B12E (R-1) |
322* | DORBDB WORK | DORGQR WORK | DORGLQ WORK | B21D (R) |
323* | | | | B21E (R-1) |
324* | | | | B22D (R) |
325* | | | | B22E (R-1) |
326* | | | | DBBCSD WORK |
327* |-------------------------------------------------------|
328*
329 IF( info .EQ. 0 ) THEN
330 iphi = 2
331 ib11d = iphi + max( 1, r-1 )
332 ib11e = ib11d + max( 1, r )
333 ib12d = ib11e + max( 1, r - 1 )
334 ib12e = ib12d + max( 1, r )
335 ib21d = ib12e + max( 1, r - 1 )
336 ib21e = ib21d + max( 1, r )
337 ib22d = ib21e + max( 1, r - 1 )
338 ib22e = ib22d + max( 1, r )
339 ibbcsd = ib22e + max( 1, r - 1 )
340 itaup1 = iphi + max( 1, r-1 )
341 itaup2 = itaup1 + max( 1, p )
342 itauq1 = itaup2 + max( 1, m-p )
343 iorbdb = itauq1 + max( 1, q )
344 iorgqr = itauq1 + max( 1, q )
345 iorglq = itauq1 + max( 1, q )
346 lorgqrmin = 1
347 lorgqropt = 1
348 lorglqmin = 1
349 lorglqopt = 1
350 IF( r .EQ. q ) THEN
351 CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
352 $ dum1, dum1, dum1, dum1, work,
353 $ -1, childinfo )
354 lorbdb = int( work(1) )
355 IF( wantu1 .AND. p .GT. 0 ) THEN
356 CALL dorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
357 $ childinfo )
358 lorgqrmin = max( lorgqrmin, p )
359 lorgqropt = max( lorgqropt, int( work(1) ) )
360 ENDIF
361 IF( wantu2 .AND. m-p .GT. 0 ) THEN
362 CALL dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
363 $ -1, childinfo )
364 lorgqrmin = max( lorgqrmin, m-p )
365 lorgqropt = max( lorgqropt, int( work(1) ) )
366 END IF
367 IF( wantv1t .AND. q .GT. 0 ) THEN
368 CALL dorglq( q-1, q-1, q-1, v1t, ldv1t,
369 $ dum1, work(1), -1, childinfo )
370 lorglqmin = max( lorglqmin, q-1 )
371 lorglqopt = max( lorglqopt, int( work(1) ) )
372 END IF
373 CALL dbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q,
374 $ theta,
375 $ dum1, u1, ldu1, u2, ldu2, v1t, ldv1t,
376 $ dum2, 1, dum1, dum1, dum1,
377 $ dum1, dum1, dum1, dum1,
378 $ dum1, work(1), -1, childinfo )
379 lbbcsd = int( work(1) )
380 ELSE IF( r .EQ. p ) THEN
381 CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
382 $ dum1, dum1, dum1, dum1,
383 $ work(1), -1, childinfo )
384 lorbdb = int( work(1) )
385 IF( wantu1 .AND. p .GT. 0 ) THEN
386 CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, dum1,
387 $ work(1), -1, childinfo )
388 lorgqrmin = max( lorgqrmin, p-1 )
389 lorgqropt = max( lorgqropt, int( work(1) ) )
390 END IF
391 IF( wantu2 .AND. m-p .GT. 0 ) THEN
392 CALL dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
393 $ -1, childinfo )
394 lorgqrmin = max( lorgqrmin, m-p )
395 lorgqropt = max( lorgqropt, int( work(1) ) )
396 END IF
397 IF( wantv1t .AND. q .GT. 0 ) THEN
398 CALL dorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
399 $ childinfo )
400 lorglqmin = max( lorglqmin, q )
401 lorglqopt = max( lorglqopt, int( work(1) ) )
402 END IF
403 CALL dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p,
404 $ theta,
405 $ dum1, v1t, ldv1t, dum2, 1, u1, ldu1,
406 $ u2, ldu2, dum1, dum1, dum1,
407 $ dum1, dum1, dum1, dum1,
408 $ dum1, work(1), -1, childinfo )
409 lbbcsd = int( work(1) )
410 ELSE IF( r .EQ. m-p ) THEN
411 CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
412 $ dum1, dum1, dum1, dum1,
413 $ work(1), -1, childinfo )
414 lorbdb = int( work(1) )
415 IF( wantu1 .AND. p .GT. 0 ) THEN
416 CALL dorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
417 $ childinfo )
418 lorgqrmin = max( lorgqrmin, p )
419 lorgqropt = max( lorgqropt, int( work(1) ) )
420 END IF
421 IF( wantu2 .AND. m-p .GT. 0 ) THEN
422 CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
423 $ dum1, work(1), -1, childinfo )
424 lorgqrmin = max( lorgqrmin, m-p-1 )
425 lorgqropt = max( lorgqropt, int( work(1) ) )
426 END IF
427 IF( wantv1t .AND. q .GT. 0 ) THEN
428 CALL dorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
429 $ childinfo )
430 lorglqmin = max( lorglqmin, q )
431 lorglqopt = max( lorglqopt, int( work(1) ) )
432 END IF
433 CALL dbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
434 $ theta, dum1, dum2, 1, v1t, ldv1t, u2,
435 $ ldu2, u1, ldu1, dum1, dum1, dum1,
436 $ dum1, dum1, dum1, dum1,
437 $ dum1, work(1), -1, childinfo )
438 lbbcsd = int( work(1) )
439 ELSE
440 CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
441 $ dum1, dum1, dum1, dum1,
442 $ dum1, work(1), -1, childinfo )
443 lorbdb = m + int( work(1) )
444 IF( wantu1 .AND. p .GT. 0 ) THEN
445 CALL dorgqr( p, p, m-q, u1, ldu1, dum1, work(1), -1,
446 $ childinfo )
447 lorgqrmin = max( lorgqrmin, p )
448 lorgqropt = max( lorgqropt, int( work(1) ) )
449 END IF
450 IF( wantu2 .AND. m-p .GT. 0 ) THEN
451 CALL dorgqr( m-p, m-p, m-q, u2, ldu2, dum1, work(1),
452 $ -1, childinfo )
453 lorgqrmin = max( lorgqrmin, m-p )
454 lorgqropt = max( lorgqropt, int( work(1) ) )
455 END IF
456 IF( wantv1t .AND. q .GT. 0 ) THEN
457 CALL dorglq( q, q, q, v1t, ldv1t, dum1, work(1), -1,
458 $ childinfo )
459 lorglqmin = max( lorglqmin, q )
460 lorglqopt = max( lorglqopt, int( work(1) ) )
461 END IF
462 CALL dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
463 $ theta, dum1, u2, ldu2, u1, ldu1, dum2,
464 $ 1, v1t, ldv1t, dum1, dum1, dum1,
465 $ dum1, dum1, dum1, dum1,
466 $ dum1, work(1), -1, childinfo )
467 lbbcsd = int( work(1) )
468 END IF
469 lworkmin = max( iorbdb+lorbdb-1,
470 $ iorgqr+lorgqrmin-1,
471 $ iorglq+lorglqmin-1,
472 $ ibbcsd+lbbcsd-1 )
473 lworkopt = max( iorbdb+lorbdb-1,
474 $ iorgqr+lorgqropt-1,
475 $ iorglq+lorglqopt-1,
476 $ ibbcsd+lbbcsd-1 )
477 work(1) = lworkopt
478 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
479 info = -19
480 END IF
481 END IF
482 IF( info .NE. 0 ) THEN
483 CALL xerbla( 'DORCSD2BY1', -info )
484 RETURN
485 ELSE IF( lquery ) THEN
486 RETURN
487 END IF
488 lorgqr = lwork-iorgqr+1
489 lorglq = lwork-iorglq+1
490*
491* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
492* in which R = MIN(P,M-P,Q,M-Q)
493*
494 IF( r .EQ. q ) THEN
495*
496* Case 1: R = Q
497*
498* Simultaneously bidiagonalize X11 and X21
499*
500 CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
501 $ work(iphi), work(itaup1), work(itaup2),
502 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
503*
504* Accumulate Householder reflectors
505*
506 IF( wantu1 .AND. p .GT. 0 ) THEN
507 CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
508 CALL dorgqr( p, p, q, u1, ldu1, work(itaup1),
509 $ work(iorgqr),
510 $ lorgqr, childinfo )
511 END IF
512 IF( wantu2 .AND. m-p .GT. 0 ) THEN
513 CALL dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
514 CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
515 $ work(iorgqr), lorgqr, childinfo )
516 END IF
517 IF( wantv1t .AND. q .GT. 0 ) THEN
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 dlacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
524 $ ldv1t )
525 CALL dorglq( q-1, q-1, q-1, v1t(2,2), ldv1t,
526 $ work(itauq1),
527 $ work(iorglq), lorglq, childinfo )
528 END IF
529*
530* Simultaneously diagonalize X11 and X21.
531*
532 CALL dbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
533 $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t,
534 $ dum2, 1, work(ib11d), work(ib11e),
535 $ work(ib12d), work(ib12e), work(ib21d),
536 $ work(ib21e), work(ib22d), work(ib22e),
537 $ work(ibbcsd), lbbcsd, childinfo )
538*
539* Permute rows and columns to place zero submatrices in
540* preferred positions
541*
542 IF( q .GT. 0 .AND. wantu2 ) THEN
543 DO i = 1, q
544 iwork(i) = m - p - q + i
545 END DO
546 DO i = q + 1, m - p
547 iwork(i) = i - q
548 END DO
549 CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
550 END IF
551 ELSE IF( r .EQ. p ) THEN
552*
553* Case 2: R = P
554*
555* Simultaneously bidiagonalize X11 and X21
556*
557 CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
558 $ work(iphi), work(itaup1), work(itaup2),
559 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
560*
561* Accumulate Householder reflectors
562*
563 IF( wantu1 .AND. p .GT. 0 ) THEN
564 u1(1,1) = one
565 DO j = 2, p
566 u1(1,j) = zero
567 u1(j,1) = zero
568 END DO
569 CALL dlacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2),
570 $ ldu1 )
571 CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
572 $ work(iorgqr), lorgqr, childinfo )
573 END IF
574 IF( wantu2 .AND. m-p .GT. 0 ) THEN
575 CALL dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
576 CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
577 $ work(iorgqr), lorgqr, childinfo )
578 END IF
579 IF( wantv1t .AND. q .GT. 0 ) THEN
580 CALL dlacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
581 CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
582 $ work(iorglq), lorglq, childinfo )
583 END IF
584*
585* Simultaneously diagonalize X11 and X21.
586*
587 CALL dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
588 $ work(iphi), v1t, ldv1t, dum1, 1, u1, ldu1, u2,
589 $ ldu2, work(ib11d), work(ib11e), work(ib12d),
590 $ work(ib12e), work(ib21d), work(ib21e),
591 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
592 $ childinfo )
593*
594* Permute rows and columns to place identity submatrices in
595* preferred positions
596*
597 IF( q .GT. 0 .AND. wantu2 ) THEN
598 DO i = 1, q
599 iwork(i) = m - p - q + i
600 END DO
601 DO i = q + 1, m - p
602 iwork(i) = i - q
603 END DO
604 CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
605 END IF
606 ELSE IF( r .EQ. m-p ) THEN
607*
608* Case 3: R = M-P
609*
610* Simultaneously bidiagonalize X11 and X21
611*
612 CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
613 $ work(iphi), work(itaup1), work(itaup2),
614 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
615*
616* Accumulate Householder reflectors
617*
618 IF( wantu1 .AND. p .GT. 0 ) THEN
619 CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
620 CALL dorgqr( p, p, q, u1, ldu1, work(itaup1),
621 $ work(iorgqr),
622 $ lorgqr, childinfo )
623 END IF
624 IF( wantu2 .AND. m-p .GT. 0 ) THEN
625 u2(1,1) = one
626 DO j = 2, m-p
627 u2(1,j) = zero
628 u2(j,1) = zero
629 END DO
630 CALL dlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
631 $ ldu2 )
632 CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
633 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
634 END IF
635 IF( wantv1t .AND. q .GT. 0 ) THEN
636 CALL dlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
637 CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
638 $ work(iorglq), lorglq, childinfo )
639 END IF
640*
641* Simultaneously diagonalize X11 and X21.
642*
643 CALL dbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
644 $ theta, work(iphi), dum1, 1, v1t, ldv1t, u2,
645 $ ldu2, u1, ldu1, work(ib11d), work(ib11e),
646 $ work(ib12d), work(ib12e), work(ib21d),
647 $ work(ib21e), work(ib22d), work(ib22e),
648 $ work(ibbcsd), lbbcsd, childinfo )
649*
650* Permute rows and columns to place identity submatrices in
651* preferred positions
652*
653 IF( q .GT. r ) THEN
654 DO i = 1, r
655 iwork(i) = q - r + i
656 END DO
657 DO i = r + 1, q
658 iwork(i) = i - r
659 END DO
660 IF( wantu1 ) THEN
661 CALL dlapmt( .false., p, q, u1, ldu1, iwork )
662 END IF
663 IF( wantv1t ) THEN
664 CALL dlapmr( .false., q, q, v1t, ldv1t, iwork )
665 END IF
666 END IF
667 ELSE
668*
669* Case 4: R = M-Q
670*
671* Simultaneously bidiagonalize X11 and X21
672*
673 CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
674 $ work(iphi), work(itaup1), work(itaup2),
675 $ work(itauq1), work(iorbdb), work(iorbdb+m),
676 $ lorbdb-m, childinfo )
677*
678* Accumulate Householder reflectors
679*
680 IF( wantu2 .AND. m-p .GT. 0 ) THEN
681 CALL dcopy( m-p, work(iorbdb+p), 1, u2, 1 )
682 END IF
683 IF( wantu1 .AND. p .GT. 0 ) THEN
684 CALL dcopy( p, work(iorbdb), 1, u1, 1 )
685 DO j = 2, p
686 u1(1,j) = zero
687 END DO
688 CALL dlacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
689 $ ldu1 )
690 CALL dorgqr( p, p, m-q, u1, ldu1, work(itaup1),
691 $ work(iorgqr), lorgqr, childinfo )
692 END IF
693 IF( wantu2 .AND. m-p .GT. 0 ) THEN
694 DO j = 2, m-p
695 u2(1,j) = zero
696 END DO
697 CALL dlacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
698 $ ldu2 )
699 CALL dorgqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
700 $ work(iorgqr), lorgqr, childinfo )
701 END IF
702 IF( wantv1t .AND. q .GT. 0 ) THEN
703 CALL dlacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
704 CALL dlacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1),
705 $ ldx11,
706 $ v1t(m-q+1,m-q+1), ldv1t )
707 CALL dlacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
708 $ v1t(p+1,p+1), ldv1t )
709 CALL dorglq( q, q, q, v1t, ldv1t, work(itauq1),
710 $ work(iorglq), lorglq, childinfo )
711 END IF
712*
713* Simultaneously diagonalize X11 and X21.
714*
715 CALL dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
716 $ theta, work(iphi), u2, ldu2, u1, ldu1, dum1,
717 $ 1, v1t, ldv1t, work(ib11d), work(ib11e),
718 $ work(ib12d), work(ib12e), work(ib21d),
719 $ work(ib21e), work(ib22d), work(ib22e),
720 $ work(ibbcsd), lbbcsd, childinfo )
721*
722* Permute rows and columns to place identity submatrices in
723* preferred positions
724*
725 IF( p .GT. r ) THEN
726 DO i = 1, r
727 iwork(i) = p - r + i
728 END DO
729 DO i = r + 1, p
730 iwork(i) = i - r
731 END DO
732 IF( wantu1 ) THEN
733 CALL dlapmt( .false., p, p, u1, ldu1, iwork )
734 END IF
735 IF( wantv1t ) THEN
736 CALL dlapmr( .false., p, q, v1t, ldv1t, iwork )
737 END IF
738 END IF
739 END IF
740*
741 RETURN
742*
743* End of DORCSD2BY1
744*
745 END
746
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dbbcsd
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:331
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dlapmr
subroutine dlapmr(forwrd, m, n, x, ldx, k)
DLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition dlapmr.f:102
dlapmt
subroutine dlapmt(forwrd, m, n, x, ldx, k)
DLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition dlapmt.f:102
dorbdb1
subroutine dorbdb1(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB1
Definition dorbdb1.f:202
dorbdb2
subroutine dorbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB2
Definition dorbdb2.f:201
dorbdb3
subroutine dorbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB3
Definition dorbdb3.f:200
dorbdb4
subroutine dorbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
DORBDB4
Definition dorbdb4.f:212
dorcsd2by1
subroutine dorcsd2by1(jobu1, jobu2, jobv1t, m, p, q, x11, ldx11, x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work, lwork, iwork, info)
DORCSD2BY1
Definition dorcsd2by1.f:232
dorglq
subroutine dorglq(m, n, k, a, lda, tau, work, lwork, info)
DORGLQ
Definition dorglq.f:125
dorgqr
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:126