OpenRadioss 2025.1.11
OpenRadioss project
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*> \htmlonly
9*> Download DORCSD2BY1 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd2by1.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd2by1.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd2by1.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DORCSD2BY1( 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* DOUBLE PRECISION THETA(*)
32* DOUBLE PRECISION 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*> DORCSD2BY1 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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: DBBCSD 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 doubleOTHERcomputational
228*
229* =====================================================================
230 SUBROUTINE dorcsd2by1( 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 (3.5.0) --
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 DOUBLE PRECISION THETA(*)
245 DOUBLE PRECISION U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
246 $ x11(ldx11,*), x21(ldx21,*)
247 INTEGER IWORK(*)
248* ..
249*
250* =====================================================================
251*
252* .. Parameters ..
253 DOUBLE PRECISION ONE, ZERO
254 PARAMETER ( ONE = 1.0d0, zero = 0.0d0 )
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 DOUBLE PRECISION DUM1(1), DUM2(1,1)
267* ..
268* .. External Subroutines ..
269 EXTERNAL dbbcsd, dcopy, dlacpy, dlapmr, dlapmt, 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.EQ. LQUERY = LWORK -1
289*
290.LT. IF( M 0 ) THEN
291 INFO = -4
292.LT..OR..GT. ELSE IF( P 0 P M ) THEN
293 INFO = -5
294.LT..OR..GT. ELSE IF( Q 0 Q M ) THEN
295 INFO = -6
296.LT. ELSE IF( LDX11 MAX( 1, P ) ) THEN
297 INFO = -8
298.LT. ELSE IF( LDX21 MAX( 1, M-P ) ) THEN
299 INFO = -10
300.AND..LT. ELSE IF( WANTU1 LDU1 MAX( 1, P ) ) THEN
301 INFO = -13
302.AND..LT. ELSE IF( WANTU2 LDU2 MAX( 1, M - P ) ) THEN
303 INFO = -15
304.AND..LT. ELSE IF( WANTV1T LDV1T 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.EQ. IF( INFO 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.EQ. IF( R 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.AND..GT. IF( WANTU1 P 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.AND..GT. IF( WANTU2 M-P 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.AND..GT. IF( WANTV1T Q 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, THETA,
374 $ DUM1, U1, LDU1, U2, LDU2, V1T, LDV1T,
375 $ DUM2, 1, DUM1, DUM1, DUM1,
376 $ DUM1, DUM1, DUM1, DUM1,
377 $ DUM1, WORK(1), -1, CHILDINFO )
378 LBBCSD = INT( WORK(1) )
379.EQ. ELSE IF( R P ) THEN
380 CALL DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
381 $ DUM1, DUM1, DUM1, DUM1,
382 $ WORK(1), -1, CHILDINFO )
383 LORBDB = INT( WORK(1) )
384.AND..GT. IF( WANTU1 P 0 ) THEN
385 CALL DORGQR( 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.AND..GT. IF( WANTU2 M-P 0 ) THEN
391 CALL DORGQR( M-P, M-P, Q, U2, LDU2, DUM1, WORK(1),
392 $ -1, CHILDINFO )
393 LORGQRMIN = MAX( LORGQRMIN, M-P )
394 LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
395 END IF
396.AND..GT. IF( WANTV1T Q 0 ) THEN
397 CALL DORGLQ( 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 DBBCSD( JOBV1T, 'n', JOBU1, JOBU2, 't', M, Q, P, THETA,
403 $ DUM1, V1T, LDV1T, DUM2, 1, U1, LDU1,
404 $ U2, LDU2, DUM1, DUM1, DUM1,
405 $ DUM1, DUM1, DUM1, DUM1,
406 $ DUM1, WORK(1), -1, CHILDINFO )
407 LBBCSD = INT( WORK(1) )
408.EQ. ELSE IF( R M-P ) THEN
409 CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
410 $ DUM1, DUM1, DUM1, DUM1,
411 $ WORK(1), -1, CHILDINFO )
412 LORBDB = INT( WORK(1) )
413.AND..GT. IF( WANTU1 P 0 ) THEN
414 CALL DORGQR( 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.AND..GT. IF( WANTU2 M-P 0 ) THEN
420 CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2,
421 $ DUM1, WORK(1), -1, CHILDINFO )
422 LORGQRMIN = MAX( LORGQRMIN, M-P-1 )
423 LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
424 END IF
425.AND..GT. IF( WANTV1T Q 0 ) THEN
426 CALL DORGLQ( 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 DBBCSD( 'n', JOBV1T, JOBU2, JOBU1, 't', M, M-Q, M-P,
432 $ THETA, DUM1, DUM2, 1, V1T, LDV1T, U2,
433 $ LDU2, U1, LDU1, DUM1, DUM1, DUM1,
434 $ DUM1, DUM1, DUM1, DUM1,
435 $ DUM1, WORK(1), -1, CHILDINFO )
436 LBBCSD = INT( WORK(1) )
437 ELSE
438 CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
439 $ DUM1, DUM1, DUM1, DUM1,
440 $ DUM1, WORK(1), -1, CHILDINFO )
441 LORBDB = M + INT( WORK(1) )
442.AND..GT. IF( WANTU1 P 0 ) THEN
443 CALL DORGQR( 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.AND..GT. IF( WANTU2 M-P 0 ) THEN
449 CALL DORGQR( 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.AND..GT. IF( WANTV1T Q 0 ) THEN
455 CALL DORGLQ( 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 DBBCSD( JOBU2, JOBU1, 'n', JOBV1T, 'n', M, M-P, M-Q,
461 $ THETA, DUM1, U2, LDU2, U1, LDU1, DUM2,
462 $ 1, V1T, LDV1T, DUM1, DUM1, DUM1,
463 $ DUM1, DUM1, DUM1, DUM1,
464 $ DUM1, WORK(1), -1, 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.LT..AND..NOT. IF( LWORK LWORKMIN LQUERY ) THEN
477 INFO = -19
478 END IF
479 END IF
480.NE. IF( INFO 0 ) THEN
481 CALL XERBLA( 'dorcsd2by1', -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.EQ. IF( R Q ) THEN
493*
494* Case 1: R = Q
495*
496* Simultaneously bidiagonalize X11 and X21
497*
498 CALL DORBDB1( 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.AND..GT. IF( WANTU1 P 0 ) THEN
505 CALL DLACPY( 'l', P, Q, X11, LDX11, U1, LDU1 )
506 CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
507 $ LORGQR, CHILDINFO )
508 END IF
509.AND..GT. IF( WANTU2 M-P 0 ) THEN
510 CALL DLACPY( 'l', M-P, Q, X21, LDX21, U2, LDU2 )
511 CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
512 $ WORK(IORGQR), LORGQR, CHILDINFO )
513 END IF
514.AND..GT. IF( WANTV1T Q 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 DLACPY( 'u', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2),
521 $ LDV1T )
522 CALL DORGLQ( 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 DBBCSD( 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),
531 $ WORK(IB12D), WORK(IB12E), WORK(IB21D),
532 $ WORK(IB21E), WORK(IB22D), WORK(IB22E),
533 $ WORK(IBBCSD), LBBCSD, 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 DLAPMT( .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 DORBDB2( 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 DLACPY( 'l', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 )
566 CALL DORGQR( 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 DLACPY( 'l', M-P, Q, X21, LDX21, U2, LDU2 )
571 CALL DORGQR( 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 DLACPY( 'u', p, q, x11, ldx11, v1t, ldv1t )
576 CALL dorglq( 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 dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
583 $ work(iphi), v1t, ldv1t, dum2, 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 IF( q .GT. 0 .AND. 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 dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
600 END IF
601 ELSE IF( r .EQ. m-p ) THEN
602*
603* Case 3: R = M-P
604*
605* Simultaneously bidiagonalize X11 and X21
606*
607 CALL dorbdb3( 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 IF( wantu1 .AND. p .GT. 0 ) THEN
614 CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
615 CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
616 $ lorgqr, childinfo )
617 END IF
618 IF( wantu2 .AND. m-p .GT. 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 dlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
625 $ ldu2 )
626 CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
627 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
628 END IF
629 IF( wantv1t .AND. q .GT. 0 ) THEN
630 CALL dlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
631 CALL dorglq( 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 dbbcsd( 'n', JOBV1T, JOBU2, JOBU1, 't', M, M-Q, M-P,
638 $ THETA, WORK(IPHI), DUM2, 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 DLAPMT( .FALSE., P, Q, U1, LDU1, IWORK )
656 END IF
657 IF( WANTV1T ) THEN
658 CALL DLAPMR( .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 DORBDB4( 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 DCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 )
676 END IF
677.AND..GT. IF( WANTU1 P 0 ) THEN
678 CALL DCOPY( P, WORK(IORBDB), 1, U1, 1 )
679 DO J = 2, P
680 U1(1,J) = ZERO
681 END DO
682 CALL DLACPY( 'l', P-1, M-Q-1, X11(2,1), LDX11, U1(2,2),
683 $ LDU1 )
684 CALL DORGQR( 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 DLACPY( 'l', M-P-1, M-Q-1, X21(2,1), LDX21, U2(2,2),
692 $ LDU2 )
693 CALL DORGQR( 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 DLACPY( 'u', M-Q, Q, X21, LDX21, V1T, LDV1T )
698 CALL DLACPY( '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 DLACPY( 'u', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21,
701 $ V1T(P+1,P+1), LDV1T )
702 CALL DORGLQ( 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 DBBCSD( JOBU2, JOBU1, 'n', JOBV1T, 'n', M, M-P, M-Q,
709 $ THETA, WORK(IPHI), U2, LDU2, U1, LDU1, DUM2,
710 $ 1, V1T, LDV1T, WORK(IB11D), WORK(IB11E),
711 $ WORK(IB12D), WORK(IB12E), WORK(IB21D),
712 $ WORK(IB21E), WORK(IB22D), WORK(IB22E),
713 $ WORK(IBBCSD), LBBCSD, 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 DLAPMT( .FALSE., P, P, U1, LDU1, IWORK )
727 END IF
728 IF( WANTV1T ) THEN
729 CALL DLAPMR( .FALSE., P, Q, V1T, LDV1T, IWORK )
730 END IF
731 END IF
732 END IF
733*
734 RETURN
735*
736* End of DORCSD2BY1
737*
738 END
739
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:103
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
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:104
dlapmr
subroutine dlapmr(forwrd, m, n, x, ldx, k)
DLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition dlapmr.f:104
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:332
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:233
dorbdb3
subroutine dorbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB3
Definition dorbdb3.f:201
dorbdb2
subroutine dorbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB2
Definition dorbdb2.f:202
dorbdb1
subroutine dorbdb1(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB1
Definition dorbdb1.f:203
dorbdb4
subroutine dorbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
DORBDB4
Definition dorbdb4.f:213
dorgqr
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:128
dorglq
subroutine dorglq(m, n, k, a, lda, tau, work, lwork, info)
DORGLQ
Definition dorglq.f:127
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21