OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sorcsd2by1.f
Go to the documentation of this file.
1*> \brief \b SORCSD2BY1
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SORCSD2BY1 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorcsd2by1.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorcsd2by1.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorcsd2by1.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
22* X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
23* LDV1T, WORK, LWORK, IWORK, INFO )
24*
25* .. Scalar Arguments ..
26* CHARACTER JOBU1, JOBU2, JOBV1T
27* INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
28* $ M, P, Q
29* ..
30* .. Array Arguments ..
31* REAL THETA(*)
32* REAL U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
33* $ X11(LDX11,*), X21(LDX21,*)
34* INTEGER IWORK(*)
35* ..
36*
37*
38*> \par Purpose:
39* =============
40*>
41*>\verbatim
42*>
43*> SORCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with
44*> orthonormal columns that has been partitioned into a 2-by-1 block
45*> structure:
46*>
47*> [ I1 0 0 ]
48*> [ 0 C 0 ]
49*> [ X11 ] [ U1 | ] [ 0 0 0 ]
50*> X = [-----] = [---------] [----------] V1**T .
51*> [ X21 ] [ | U2 ] [ 0 0 0 ]
52*> [ 0 S 0 ]
53*> [ 0 0 I2]
54*>
55*> X11 is P-by-Q. The orthogonal matrices U1, U2, and V1 are P-by-P,
56*> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R
57*> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which
58*> R = MIN(P,M-P,Q,M-Q). I1 is a K1-by-K1 identity matrix and I2 is a
59*> K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0).
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] M
87*> \verbatim
88*> M is INTEGER
89*> The number of rows in X.
90*> \endverbatim
91*>
92*> \param[in] P
93*> \verbatim
94*> P is INTEGER
95*> The number of rows in X11. 0 <= P <= M.
96*> \endverbatim
97*>
98*> \param[in] Q
99*> \verbatim
100*> Q is INTEGER
101*> The number of columns in X11 and X21. 0 <= Q <= M.
102*> \endverbatim
103*>
104*> \param[in,out] X11
105*> \verbatim
106*> X11 is REAL array, dimension (LDX11,Q)
107*> On entry, part of the orthogonal matrix whose CSD is desired.
108*> \endverbatim
109*>
110*> \param[in] LDX11
111*> \verbatim
112*> LDX11 is INTEGER
113*> The leading dimension of X11. LDX11 >= MAX(1,P).
114*> \endverbatim
115*>
116*> \param[in,out] X21
117*> \verbatim
118*> X21 is REAL array, dimension (LDX21,Q)
119*> On entry, part of the orthogonal matrix whose CSD is desired.
120*> \endverbatim
121*>
122*> \param[in] LDX21
123*> \verbatim
124*> LDX21 is INTEGER
125*> The leading dimension of X21. LDX21 >= MAX(1,M-P).
126*> \endverbatim
127*>
128*> \param[out] THETA
129*> \verbatim
130*> THETA is REAL array, dimension (R), in which R =
131*> MIN(P,M-P,Q,M-Q).
132*> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
133*> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
134*> \endverbatim
135*>
136*> \param[out] U1
137*> \verbatim
138*> U1 is REAL array, dimension (P)
139*> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
140*> \endverbatim
141*>
142*> \param[in] LDU1
143*> \verbatim
144*> LDU1 is INTEGER
145*> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
146*> MAX(1,P).
147*> \endverbatim
148*>
149*> \param[out] U2
150*> \verbatim
151*> U2 is REAL array, dimension (M-P)
152*> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
153*> matrix U2.
154*> \endverbatim
155*>
156*> \param[in] LDU2
157*> \verbatim
158*> LDU2 is INTEGER
159*> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
160*> MAX(1,M-P).
161*> \endverbatim
162*>
163*> \param[out] V1T
164*> \verbatim
165*> V1T is REAL array, dimension (Q)
166*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
167*> matrix V1**T.
168*> \endverbatim
169*>
170*> \param[in] LDV1T
171*> \verbatim
172*> LDV1T is INTEGER
173*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
174*> MAX(1,Q).
175*> \endverbatim
176*>
177*> \param[out] WORK
178*> \verbatim
179*> WORK is REAL array, dimension (MAX(1,LWORK))
180*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
181*> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
182*> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
183*> define the matrix in intermediate bidiagonal-block form
184*> remaining after nonconvergence. INFO specifies the number
185*> of nonzero PHI's.
186*> \endverbatim
187*>
188*> \param[in] LWORK
189*> \verbatim
190*> LWORK is INTEGER
191*> The dimension of the array WORK.
192*>
193*> If LWORK = -1, then a workspace query is assumed; the routine
194*> only calculates the optimal size of the WORK array, returns
195*> this value as the first entry of the work array, and no error
196*> message related to LWORK is issued by XERBLA.
197*> \endverbatim
198*>
199*> \param[out] IWORK
200*> \verbatim
201*> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
202*> \endverbatim
203*>
204*> \param[out] INFO
205*> \verbatim
206*> INFO is INTEGER
207*> = 0: successful exit.
208*> < 0: if INFO = -i, the i-th argument had an illegal value.
209*> > 0: SBBCSD did not converge. See the description of WORK
210*> above for details.
211*> \endverbatim
212*
213*> \par References:
214* ================
215*>
216*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
217*> Algorithms, 50(1):33-65, 2009.
218*
219* Authors:
220* ========
221*
222*> \author Univ. of Tennessee
223*> \author Univ. of California Berkeley
224*> \author Univ. of Colorado Denver
225*> \author NAG Ltd.
226*
227*> \ingroup realOTHERcomputational
228*
229* =====================================================================
230 SUBROUTINE sorcsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
231 $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
232 $ LDV1T, WORK, LWORK, IWORK, INFO )
233*
234* -- LAPACK computational routine --
235* -- LAPACK is a software package provided by Univ. of Tennessee, --
236* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
237*
238* .. Scalar Arguments ..
239 CHARACTER JOBU1, JOBU2, JOBV1T
240 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
241 $ m, p, q
242* ..
243* .. Array Arguments ..
244 REAL THETA(*)
245 REAL U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
246 $ x11(ldx11,*), x21(ldx21,*)
247 INTEGER IWORK(*)
248* ..
249*
250* =====================================================================
251*
252* .. Parameters ..
253 REAL ONE, ZERO
254 PARAMETER ( ONE = 1.0e0, zero = 0.0e0 )
255* ..
256* .. Local Scalars ..
257 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
258 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
259 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
260 $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
261 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
262 $ lworkmin, lworkopt, r
263 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
264* ..
265* .. Local Arrays ..
266 REAL DUM1(1), DUM2(1,1)
267* ..
268* .. External Subroutines ..
269 EXTERNAL sbbcsd, scopy, slacpy, slapmr, slapmt, sorbdb1,
270 $ sorbdb2, sorbdb3, sorbdb4, sorglq, sorgqr,
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* | SORBDB WORK | SORGQR WORK | SORGLQ WORK | B21D (R) |
323* | | | | B21E (R-1) |
324* | | | | B22D (R) |
325* | | | | B22E (R-1) |
326* | | | | SBBCSD 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 sorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
352 $ dum1, dum1, dum1, dum1, work, -1,
353 $ childinfo )
354 lorbdb = int( work(1) )
355 IF( wantu1 .AND. p .GT. 0 ) THEN
356 CALL sorgqr( 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 sorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1), -1,
363 $ 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 sorglq( 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 sbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
374 $ dum1, u1, ldu1, u2, ldu2, v1t, ldv1t, dum2,
375 $ 1, dum1, dum1, dum1, dum1, dum1,
376 $ dum1, dum1, dum1, work(1), -1, childinfo
377 $ )
378 lbbcsd = int( work(1) )
379 ELSE IF( r .EQ. p ) THEN
380 CALL sorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
381 $ dum1, dum1, dum1, dum1, work(1), -1,
382 $ childinfo )
383 lorbdb = int( work(1) )
384 IF( wantu1 .AND. p .GT. 0 ) THEN
385 CALL sorgqr( p-1, p-1, p-1, u1(2,2), ldu1, dum1,
386 $ work(1), -1, childinfo )
387 lorgqrmin = max( lorgqrmin, p-1 )
388 lorgqropt = max( lorgqropt, int( work(1) ) )
389 END IF
390 IF( wantu2 .AND. m-p .GT. 0 ) THEN
391 CALL sorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1), -1,
392 $ childinfo )
393 lorgqrmin = max( lorgqrmin, m-p )
394 lorgqropt = max( lorgqropt, int( work(1) ) )
395 END IF
396 IF( wantv1t .AND. q .GT. 0 ) THEN
397 CALL sorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
398 $ childinfo )
399 lorglqmin = max( lorglqmin, q )
400 lorglqopt = max( lorglqopt, int( work(1) ) )
401 END IF
402 CALL sbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
403 $ dum1, v1t, ldv1t, dum2, 1, u1, ldu1, u2,
404 $ ldu2, dum1, dum1, dum1, dum1, dum1,
405 $ dum1, dum1, dum1, work(1), -1, childinfo
406 $ )
407 lbbcsd = int( work(1) )
408 ELSE IF( r .EQ. m-p ) THEN
409 CALL sorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
410 $ dum1, dum1, dum1, dum1, work(1), -1,
411 $ childinfo )
412 lorbdb = int( work(1) )
413 IF( wantu1 .AND. p .GT. 0 ) THEN
414 CALL sorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
415 $ childinfo )
416 lorgqrmin = max( lorgqrmin, p )
417 lorgqropt = max( lorgqropt, int( work(1) ) )
418 END IF
419 IF( wantu2 .AND. m-p .GT. 0 ) THEN
420 CALL sorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, dum1,
421 $ work(1), -1, childinfo )
422 lorgqrmin = max( lorgqrmin, m-p-1 )
423 lorgqropt = max( lorgqropt, int( work(1) ) )
424 END IF
425 IF( wantv1t .AND. q .GT. 0 ) THEN
426 CALL sorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
427 $ childinfo )
428 lorglqmin = max( lorglqmin, q )
429 lorglqopt = max( lorglqopt, int( work(1) ) )
430 END IF
431 CALL sbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
432 $ theta, dum1, dum2, 1, v1t, ldv1t, u2, ldu2,
433 $ u1, ldu1, dum1, dum1, dum1, dum1,
434 $ dum1, dum1, dum1, dum1, work(1), -1,
435 $ childinfo )
436 lbbcsd = int( work(1) )
437 ELSE
438 CALL sorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
439 $ dum1, dum1, dum1, dum1, dum1,
440 $ work(1), -1, childinfo )
441 lorbdb = m + int( work(1) )
442 IF( wantu1 .AND. p .GT. 0 ) THEN
443 CALL sorgqr( p, p, m-q, u1, ldu1, dum1, work(1), -1,
444 $ childinfo )
445 lorgqrmin = max( lorgqrmin, p )
446 lorgqropt = max( lorgqropt, int( work(1) ) )
447 END IF
448 IF( wantu2 .AND. m-p .GT. 0 ) THEN
449 CALL sorgqr( m-p, m-p, m-q, u2, ldu2, dum1, work(1),
450 $ -1, childinfo )
451 lorgqrmin = max( lorgqrmin, m-p )
452 lorgqropt = max( lorgqropt, int( work(1) ) )
453 END IF
454 IF( wantv1t .AND. q .GT. 0 ) THEN
455 CALL sorglq( q, q, q, v1t, ldv1t, dum1, work(1), -1,
456 $ childinfo )
457 lorglqmin = max( lorglqmin, q )
458 lorglqopt = max( lorglqopt, int( work(1) ) )
459 END IF
460 CALL sbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
461 $ theta, dum1, u2, ldu2, u1, ldu1, dum2, 1,
462 $ v1t, ldv1t, dum1, dum1, dum1, dum1,
463 $ dum1, dum1, dum1, dum1, work(1), -1,
464 $ childinfo )
465 lbbcsd = int( work(1) )
466 END IF
467 lworkmin = max( iorbdb+lorbdb-1,
468 $ iorgqr+lorgqrmin-1,
469 $ iorglq+lorglqmin-1,
470 $ ibbcsd+lbbcsd-1 )
471 lworkopt = max( iorbdb+lorbdb-1,
472 $ iorgqr+lorgqropt-1,
473 $ iorglq+lorglqopt-1,
474 $ ibbcsd+lbbcsd-1 )
475 work(1) = lworkopt
476 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
477 info = -19
478 END IF
479 END IF
480 IF( info .NE. 0 ) THEN
481 CALL xerbla( 'SORCSD2BY1', -info )
482 RETURN
483 ELSE IF( lquery ) THEN
484 RETURN
485 END IF
486 lorgqr = lwork-iorgqr+1
487 lorglq = lwork-iorglq+1
488*
489* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
490* in which R = MIN(P,M-P,Q,M-Q)
491*
492 IF( r .EQ. q ) THEN
493*
494* Case 1: R = Q
495*
496* Simultaneously bidiagonalize X11 and X21
497*
498 CALL sorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
499 $ work(iphi), work(itaup1), work(itaup2),
500 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
501*
502* Accumulate Householder reflectors
503*
504 IF( wantu1 .AND. p .GT. 0 ) THEN
505 CALL slacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
506 CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
507 $ lorgqr, childinfo )
508 END IF
509 IF( wantu2 .AND. m-p .GT. 0 ) THEN
510 CALL slacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
511 CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
512 $ work(iorgqr), lorgqr, childinfo )
513 END IF
514 IF( wantv1t .AND. q .GT. 0 ) THEN
515 v1t(1,1) = one
516 DO j = 2, q
517 v1t(1,j) = zero
518 v1t(j,1) = zero
519 END DO
520 CALL slacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
521 $ ldv1t )
522 CALL sorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
523 $ work(iorglq), lorglq, childinfo )
524 END IF
525*
526* Simultaneously diagonalize X11 and X21.
527*
528 CALL sbbcsd( jobu1, jobu2, jobv1t, 'n', 'n', M, P, Q, THETA,
529 $ WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T,
530 $ DUM2, 1, WORK(IB11D), WORK(IB11E), WORK(IB12D),
531 $ WORK(IB12E), WORK(IB21D), WORK(IB21E),
532 $ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD,
533 $ CHILDINFO )
534*
535* Permute rows and columns to place zero submatrices in
536* preferred positions
537*
538.GT..AND. IF( Q 0 WANTU2 ) THEN
539 DO I = 1, Q
540 IWORK(I) = M - P - Q + I
541 END DO
542 DO I = Q + 1, M - P
543 IWORK(I) = I - Q
544 END DO
545 CALL SLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
546 END IF
547.EQ. ELSE IF( R P ) THEN
548*
549* Case 2: R = P
550*
551* Simultaneously bidiagonalize X11 and X21
552*
553 CALL SORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
554 $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
555 $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
556*
557* Accumulate Householder reflectors
558*
559.AND..GT. IF( WANTU1 P 0 ) THEN
560 U1(1,1) = ONE
561 DO J = 2, P
562 U1(1,J) = ZERO
563 U1(J,1) = ZERO
564 END DO
565 CALL SLACPY( 'l', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 )
566 CALL SORGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1),
567 $ WORK(IORGQR), LORGQR, CHILDINFO )
568 END IF
569.AND..GT. IF( WANTU2 M-P 0 ) THEN
570 CALL SLACPY( 'l', M-P, Q, X21, LDX21, U2, LDU2 )
571 CALL SORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
572 $ WORK(IORGQR), LORGQR, CHILDINFO )
573 END IF
574.AND..GT. IF( WANTV1T Q 0 ) THEN
575 CALL SLACPY( 'u', P, Q, X11, LDX11, V1T, LDV1T )
576 CALL SORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
577 $ WORK(IORGLQ), LORGLQ, CHILDINFO )
578 END IF
579*
580* Simultaneously diagonalize X11 and X21.
581*
582 CALL SBBCSD( JOBV1T, 'n', JOBU1, JOBU2, 't', M, Q, P, THETA,
583 $ WORK(IPHI), V1T, LDV1T, DUM1, 1, U1, LDU1, U2,
584 $ LDU2, WORK(IB11D), WORK(IB11E), WORK(IB12D),
585 $ WORK(IB12E), WORK(IB21D), WORK(IB21E),
586 $ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD,
587 $ CHILDINFO )
588*
589* Permute rows and columns to place identity submatrices in
590* preferred positions
591*
592.GT..AND. IF( Q 0 WANTU2 ) THEN
593 DO I = 1, Q
594 IWORK(I) = M - P - Q + I
595 END DO
596 DO I = Q + 1, M - P
597 IWORK(I) = I - Q
598 END DO
599 CALL SLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
600 END IF
601.EQ. ELSE IF( R M-P ) THEN
602*
603* Case 3: R = M-P
604*
605* Simultaneously bidiagonalize X11 and X21
606*
607 CALL SORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
608 $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
609 $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
610*
611* Accumulate Householder reflectors
612*
613.AND..GT. IF( WANTU1 P 0 ) THEN
614 CALL SLACPY( 'l', P, Q, X11, LDX11, U1, LDU1 )
615 CALL SORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
616 $ LORGQR, CHILDINFO )
617 END IF
618.AND..GT. IF( WANTU2 M-P 0 ) THEN
619 U2(1,1) = ONE
620 DO J = 2, M-P
621 U2(1,J) = ZERO
622 U2(J,1) = ZERO
623 END DO
624 CALL SLACPY( 'l', M-P-1, M-P-1, X21(2,1), LDX21, U2(2,2),
625 $ LDU2 )
626 CALL SORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2,
627 $ WORK(ITAUP2), WORK(IORGQR), LORGQR, CHILDINFO )
628 END IF
629.AND..GT. IF( WANTV1T Q 0 ) THEN
630 CALL SLACPY( 'u', M-P, Q, X21, LDX21, V1T, LDV1T )
631 CALL SORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
632 $ WORK(IORGLQ), LORGLQ, CHILDINFO )
633 END IF
634*
635* Simultaneously diagonalize X11 and X21.
636*
637 CALL SBBCSD( 'n', JOBV1T, JOBU2, JOBU1, 't', M, M-Q, M-P,
638 $ THETA, WORK(IPHI), DUM1, 1, V1T, LDV1T, U2,
639 $ LDU2, U1, LDU1, WORK(IB11D), WORK(IB11E),
640 $ WORK(IB12D), WORK(IB12E), WORK(IB21D),
641 $ WORK(IB21E), WORK(IB22D), WORK(IB22E),
642 $ WORK(IBBCSD), LBBCSD, CHILDINFO )
643*
644* Permute rows and columns to place identity submatrices in
645* preferred positions
646*
647.GT. IF( Q R ) THEN
648 DO I = 1, R
649 IWORK(I) = Q - R + I
650 END DO
651 DO I = R + 1, Q
652 IWORK(I) = I - R
653 END DO
654 IF( WANTU1 ) THEN
655 CALL SLAPMT( .FALSE., P, Q, U1, LDU1, IWORK )
656 END IF
657 IF( WANTV1T ) THEN
658 CALL SLAPMR( .FALSE., Q, Q, V1T, LDV1T, IWORK )
659 END IF
660 END IF
661 ELSE
662*
663* Case 4: R = M-Q
664*
665* Simultaneously bidiagonalize X11 and X21
666*
667 CALL SORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
668 $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
669 $ WORK(ITAUQ1), WORK(IORBDB), WORK(IORBDB+M),
670 $ LORBDB-M, CHILDINFO )
671*
672* Accumulate Householder reflectors
673*
674.AND..GT. IF( WANTU2 M-P 0 ) THEN
675 CALL SCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 )
676 END IF
677.AND..GT. IF( WANTU1 P 0 ) THEN
678 CALL SCOPY( P, WORK(IORBDB), 1, U1, 1 )
679 DO J = 2, P
680 U1(1,J) = ZERO
681 END DO
682 CALL SLACPY( 'l', P-1, M-Q-1, X11(2,1), LDX11, U1(2,2),
683 $ LDU1 )
684 CALL SORGQR( P, P, M-Q, U1, LDU1, WORK(ITAUP1),
685 $ WORK(IORGQR), LORGQR, CHILDINFO )
686 END IF
687.AND..GT. IF( WANTU2 M-P 0 ) THEN
688 DO J = 2, M-P
689 U2(1,J) = ZERO
690 END DO
691 CALL SLACPY( 'l', M-P-1, M-Q-1, X21(2,1), LDX21, U2(2,2),
692 $ LDU2 )
693 CALL SORGQR( M-P, M-P, M-Q, U2, LDU2, WORK(ITAUP2),
694 $ WORK(IORGQR), LORGQR, CHILDINFO )
695 END IF
696.AND..GT. IF( WANTV1T Q 0 ) THEN
697 CALL SLACPY( 'u', M-Q, Q, X21, LDX21, V1T, LDV1T )
698 CALL SLACPY( 'u', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11,
699 $ V1T(M-Q+1,M-Q+1), LDV1T )
700 CALL SLACPY( 'u', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21,
701 $ V1T(P+1,P+1), LDV1T )
702 CALL SORGLQ( Q, Q, Q, V1T, LDV1T, WORK(ITAUQ1),
703 $ WORK(IORGLQ), LORGLQ, CHILDINFO )
704 END IF
705*
706* Simultaneously diagonalize X11 and X21.
707*
708 CALL SBBCSD( JOBU2, JOBU1, 'n', JOBV1T, 'n', M, M-P, M-Q,
709 $ THETA, WORK(IPHI), U2, LDU2, U1, LDU1, DUM1, 1,
710 $ V1T, LDV1T, WORK(IB11D), WORK(IB11E), WORK(IB12D),
711 $ WORK(IB12E), WORK(IB21D), WORK(IB21E),
712 $ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD,
713 $ CHILDINFO )
714*
715* Permute rows and columns to place identity submatrices in
716* preferred positions
717*
718.GT. IF( P R ) THEN
719 DO I = 1, R
720 IWORK(I) = P - R + I
721 END DO
722 DO I = R + 1, P
723 IWORK(I) = I - R
724 END DO
725 IF( WANTU1 ) THEN
726 CALL SLAPMT( .FALSE., P, P, U1, LDU1, IWORK )
727 END IF
728 IF( WANTV1T ) THEN
729 CALL SLAPMR( .FALSE., P, Q, V1T, LDV1T, IWORK )
730 END IF
731 END IF
732 END IF
733*
734 RETURN
735*
736* End of SORCSD2BY1
737*
738 END
739
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:103
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
slapmr
subroutine slapmr(forwrd, m, n, x, ldx, k)
SLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition slapmr.f:104
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:104
sorbdb1
subroutine sorbdb1(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
SORBDB1
Definition sorbdb1.f:203
sorbdb3
subroutine sorbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
SORBDB3
Definition sorbdb3.f:202
sorbdb2
subroutine sorbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
SORBDB2
Definition sorbdb2.f:201
sorcsd2by1
subroutine sorcsd2by1(jobu1, jobu2, jobv1t, m, p, q, x11, ldx11, x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work, lwork, iwork, info)
SORCSD2BY1
Definition sorcsd2by1.f:233
sorglq
subroutine sorglq(m, n, k, a, lda, tau, work, lwork, info)
SORGLQ
Definition sorglq.f:127
sorbdb4
subroutine sorbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
SORBDB4
Definition sorbdb4.f:214
sorgqr
subroutine sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
Definition sorgqr.f:128
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:332
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21