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: *> Purpose:
43: *> ========
44: *>
45: *> DORCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with
46: *> orthonormal columns that has been partitioned into a 2-by-1 block
47: *> structure:
48: *>
49: *> [ I 0 0 ]
50: *> [ 0 C 0 ]
51: *> [ X11 ] [ U1 | ] [ 0 0 0 ]
52: *> X = [-----] = [---------] [----------] V1**T .
53: *> [ X21 ] [ | U2 ] [ 0 0 0 ]
54: *> [ 0 S 0 ]
55: *> [ 0 0 I ]
56: *>
57: *> X11 is P-by-Q. The orthogonal matrices U1, U2, and V1 are P-by-P,
58: *> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R
59: *> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which
60: *> R = MIN(P,M-P,Q,M-Q).
61: *> \endverbatim
62: *
63: * Arguments:
64: * ==========
65: *
66: *> \param[in] JOBU1
67: *> \verbatim
68: *> JOBU1 is CHARACTER
69: *> = 'Y': U1 is computed;
70: *> otherwise: U1 is not computed.
71: *> \endverbatim
72: *>
73: *> \param[in] JOBU2
74: *> \verbatim
75: *> JOBU2 is CHARACTER
76: *> = 'Y': U2 is computed;
77: *> otherwise: U2 is not computed.
78: *> \endverbatim
79: *>
80: *> \param[in] JOBV1T
81: *> \verbatim
82: *> JOBV1T is CHARACTER
83: *> = 'Y': V1T is computed;
84: *> otherwise: V1T is not computed.
85: *> \endverbatim
86: *>
87: *> \param[in] M
88: *> \verbatim
89: *> M is INTEGER
90: *> The number of rows in X.
91: *> \endverbatim
92: *>
93: *> \param[in] P
94: *> \verbatim
95: *> P is INTEGER
96: *> The number of rows in X11. 0 <= P <= M.
97: *> \endverbatim
98: *>
99: *> \param[in] Q
100: *> \verbatim
101: *> Q is INTEGER
102: *> The number of columns in X11 and X21. 0 <= Q <= M.
103: *> \endverbatim
104: *>
105: *> \param[in,out] X11
106: *> \verbatim
107: *> X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
108: *> On entry, part of the orthogonal matrix whose CSD is desired.
109: *> \endverbatim
110: *>
111: *> \param[in] LDX11
112: *> \verbatim
113: *> LDX11 is INTEGER
114: *> The leading dimension of X11. LDX11 >= MAX(1,P).
115: *> \endverbatim
116: *>
117: *> \param[in,out] X21
118: *> \verbatim
119: *> X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
120: *> On entry, part of the orthogonal matrix whose CSD is desired.
121: *> \endverbatim
122: *>
123: *> \param[in] LDX21
124: *> \verbatim
125: *> LDX21 is INTEGER
126: *> The leading dimension of X21. LDX21 >= MAX(1,M-P).
127: *> \endverbatim
128: *>
129: *> \param[out] THETA
130: *> \verbatim
131: *> THETA is DOUBLE PRECISION array, dimension (R), in which R =
132: *> MIN(P,M-P,Q,M-Q).
133: *> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
134: *> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
135: *> \endverbatim
136: *>
137: *> \param[out] U1
138: *> \verbatim
139: *> U1 is DOUBLE PRECISION array, dimension (P)
140: *> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
141: *> \endverbatim
142: *>
143: *> \param[in] LDU1
144: *> \verbatim
145: *> LDU1 is INTEGER
146: *> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
147: *> MAX(1,P).
148: *> \endverbatim
149: *>
150: *> \param[out] U2
151: *> \verbatim
152: *> U2 is DOUBLE PRECISION array, dimension (M-P)
153: *> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
154: *> matrix U2.
155: *> \endverbatim
156: *>
157: *> \param[in] LDU2
158: *> \verbatim
159: *> LDU2 is INTEGER
160: *> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
161: *> MAX(1,M-P).
162: *> \endverbatim
163: *>
164: *> \param[out] V1T
165: *> \verbatim
166: *> V1T is DOUBLE PRECISION array, dimension (Q)
167: *> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
168: *> matrix V1**T.
169: *> \endverbatim
170: *>
171: *> \param[in] LDV1T
172: *> \verbatim
173: *> LDV1T is INTEGER
174: *> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
175: *> MAX(1,Q).
176: *> \endverbatim
177: *>
178: *> \param[out] WORK
179: *> \verbatim
180: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
181: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182: *> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
183: *> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
184: *> define the matrix in intermediate bidiagonal-block form
185: *> remaining after nonconvergence. INFO specifies the number
186: *> of nonzero PHI's.
187: *> \endverbatim
188: *>
189: *> \param[in] LWORK
190: *> \verbatim
191: *> LWORK is INTEGER
192: *> The dimension of the array WORK.
193: *>
194: *> If LWORK = -1, then a workspace query is assumed; the routine
195: *> only calculates the optimal size of the WORK array, returns
196: *> this value as the first entry of the work array, and no error
197: *> message related to LWORK is issued by XERBLA.
198: *> \endverbatim
199: *>
200: *> \param[out] IWORK
201: *> \verbatim
202: *> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
203: *> \endverbatim
204: *>
205: *> \param[out] INFO
206: *> \verbatim
207: *> INFO is INTEGER
208: *> = 0: successful exit.
209: *> < 0: if INFO = -i, the i-th argument had an illegal value.
210: *> > 0: DBBCSD did not converge. See the description of WORK
211: *> above for details.
212: *> \endverbatim
213: *
214: *> \par References:
215: * ================
216: *>
217: *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
218: *> Algorithms, 50(1):33-65, 2009.
219: *
220: * Authors:
221: * ========
222: *
223: *> \author Univ. of Tennessee
224: *> \author Univ. of California Berkeley
225: *> \author Univ. of Colorado Denver
226: *> \author NAG Ltd.
227: *
228: *> \date July 2012
229: *
230: *> \ingroup doubleOTHERcomputational
231: *
232: * =====================================================================
233: SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
234: $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
235: $ LDV1T, WORK, LWORK, IWORK, INFO )
236: *
237: * -- LAPACK computational routine (3.5.0) --
238: * -- LAPACK is a software package provided by Univ. of Tennessee, --
239: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
240: * July 2012
241: *
242: * .. Scalar Arguments ..
243: CHARACTER JOBU1, JOBU2, JOBV1T
244: INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
245: $ M, P, Q
246: * ..
247: * .. Array Arguments ..
248: DOUBLE PRECISION THETA(*)
249: DOUBLE PRECISION U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
250: $ X11(LDX11,*), X21(LDX21,*)
251: INTEGER IWORK(*)
252: * ..
253: *
254: * =====================================================================
255: *
256: * .. Parameters ..
257: DOUBLE PRECISION ONE, ZERO
258: PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 )
259: * ..
260: * .. Local Scalars ..
261: INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
262: $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
263: $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
264: $ J, LBBCSD, LORBDB, LORGLQ, LORGLQMIN,
265: $ LORGLQOPT, LORGQR, LORGQRMIN, LORGQROPT,
266: $ LWORKMIN, LWORKOPT, R
267: LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
268: * ..
269: * .. Local Arrays ..
270: DOUBLE PRECISION DUM1(1), DUM2(1,1)
271: * ..
272: * .. External Subroutines ..
273: EXTERNAL DBBCSD, DCOPY, DLACPY, DLAPMR, DLAPMT, DORBDB1,
274: $ DORBDB2, DORBDB3, DORBDB4, DORGLQ, DORGQR,
275: $ XERBLA
276: * ..
277: * .. External Functions ..
278: LOGICAL LSAME
279: EXTERNAL LSAME
280: * ..
281: * .. Intrinsic Function ..
282: INTRINSIC INT, MAX, MIN
283: * ..
284: * .. Executable Statements ..
285: *
286: * Test input arguments
287: *
288: INFO = 0
289: WANTU1 = LSAME( JOBU1, 'Y' )
290: WANTU2 = LSAME( JOBU2, 'Y' )
291: WANTV1T = LSAME( JOBV1T, 'Y' )
292: LQUERY = LWORK .EQ. -1
293: *
294: IF( M .LT. 0 ) THEN
295: INFO = -4
296: ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
297: INFO = -5
298: ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
299: INFO = -6
300: ELSE IF( LDX11 .LT. MAX( 1, P ) ) THEN
301: INFO = -8
302: ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN
303: INFO = -10
304: ELSE IF( WANTU1 .AND. LDU1 .LT. MAX( 1, P ) ) THEN
305: INFO = -13
306: ELSE IF( WANTU2 .AND. LDU2 .LT. MAX( 1, M - P ) ) THEN
307: INFO = -15
308: ELSE IF( WANTV1T .AND. LDV1T .LT. MAX( 1, Q ) ) THEN
309: INFO = -17
310: END IF
311: *
312: R = MIN( P, M-P, Q, M-Q )
313: *
314: * Compute workspace
315: *
316: * WORK layout:
317: * |-------------------------------------------------------|
318: * | LWORKOPT (1) |
319: * |-------------------------------------------------------|
320: * | PHI (MAX(1,R-1)) |
321: * |-------------------------------------------------------|
322: * | TAUP1 (MAX(1,P)) | B11D (R) |
323: * | TAUP2 (MAX(1,M-P)) | B11E (R-1) |
324: * | TAUQ1 (MAX(1,Q)) | B12D (R) |
325: * |-----------------------------------------| B12E (R-1) |
326: * | DORBDB WORK | DORGQR WORK | DORGLQ WORK | B21D (R) |
327: * | | | | B21E (R-1) |
328: * | | | | B22D (R) |
329: * | | | | B22E (R-1) |
330: * | | | | DBBCSD WORK |
331: * |-------------------------------------------------------|
332: *
333: IF( INFO .EQ. 0 ) THEN
334: IPHI = 2
335: IB11D = IPHI + MAX( 1, R-1 )
336: IB11E = IB11D + MAX( 1, R )
337: IB12D = IB11E + MAX( 1, R - 1 )
338: IB12E = IB12D + MAX( 1, R )
339: IB21D = IB12E + MAX( 1, R - 1 )
340: IB21E = IB21D + MAX( 1, R )
341: IB22D = IB21E + MAX( 1, R - 1 )
342: IB22E = IB22D + MAX( 1, R )
343: IBBCSD = IB22E + MAX( 1, R - 1 )
344: ITAUP1 = IPHI + MAX( 1, R-1 )
345: ITAUP2 = ITAUP1 + MAX( 1, P )
346: ITAUQ1 = ITAUP2 + MAX( 1, M-P )
347: IORBDB = ITAUQ1 + MAX( 1, Q )
348: IORGQR = ITAUQ1 + MAX( 1, Q )
349: IORGLQ = ITAUQ1 + MAX( 1, Q )
350: LORGQRMIN = 1
351: LORGQROPT = 1
352: LORGLQMIN = 1
353: LORGLQOPT = 1
354: IF( R .EQ. Q ) THEN
355: CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA,
356: $ DUM1, DUM1, DUM1, DUM1, WORK,
357: $ -1, CHILDINFO )
358: LORBDB = INT( WORK(1) )
359: IF( WANTU1 .AND. P .GT. 0 ) THEN
360: CALL DORGQR( P, P, Q, U1, LDU1, DUM1, WORK(1), -1,
361: $ CHILDINFO )
362: LORGQRMIN = MAX( LORGQRMIN, P )
363: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
364: ENDIF
365: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
366: CALL DORGQR( M-P, M-P, Q, U2, LDU2, DUM1, WORK(1),
367: $ -1, CHILDINFO )
368: LORGQRMIN = MAX( LORGQRMIN, M-P )
369: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
370: END IF
371: IF( WANTV1T .AND. Q .GT. 0 ) THEN
372: CALL DORGLQ( Q-1, Q-1, Q-1, V1T, LDV1T,
373: $ DUM1, WORK(1), -1, CHILDINFO )
374: LORGLQMIN = MAX( LORGLQMIN, Q-1 )
375: LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
376: END IF
377: CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
378: $ DUM1, U1, LDU1, U2, LDU2, V1T, LDV1T,
379: $ DUM2, 1, DUM1, DUM1, DUM1,
380: $ DUM1, DUM1, DUM1, DUM1,
381: $ DUM1, WORK(1), -1, CHILDINFO )
382: LBBCSD = INT( WORK(1) )
383: ELSE IF( R .EQ. P ) THEN
384: CALL DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
385: $ DUM1, DUM1, DUM1, DUM1,
386: $ WORK(1), -1, CHILDINFO )
387: LORBDB = INT( WORK(1) )
388: IF( WANTU1 .AND. P .GT. 0 ) THEN
389: CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, DUM1,
390: $ WORK(1), -1, CHILDINFO )
391: LORGQRMIN = MAX( LORGQRMIN, P-1 )
392: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
393: END IF
394: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
395: CALL DORGQR( M-P, M-P, Q, U2, LDU2, DUM1, WORK(1),
396: $ -1, CHILDINFO )
397: LORGQRMIN = MAX( LORGQRMIN, M-P )
398: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
399: END IF
400: IF( WANTV1T .AND. Q .GT. 0 ) THEN
401: CALL DORGLQ( Q, Q, R, V1T, LDV1T, DUM1, WORK(1), -1,
402: $ CHILDINFO )
403: LORGLQMIN = MAX( LORGLQMIN, Q )
404: LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
405: END IF
406: CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
407: $ DUM1, V1T, LDV1T, DUM2, 1, U1, LDU1,
408: $ U2, LDU2, DUM1, DUM1, DUM1,
409: $ DUM1, DUM1, DUM1, DUM1,
410: $ DUM1, WORK(1), -1, CHILDINFO )
411: LBBCSD = INT( WORK(1) )
412: ELSE IF( R .EQ. M-P ) THEN
413: CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
414: $ DUM1, DUM1, DUM1, DUM1,
415: $ WORK(1), -1, CHILDINFO )
416: LORBDB = INT( WORK(1) )
417: IF( WANTU1 .AND. P .GT. 0 ) THEN
418: CALL DORGQR( P, P, Q, U1, LDU1, DUM1, WORK(1), -1,
419: $ CHILDINFO )
420: LORGQRMIN = MAX( LORGQRMIN, P )
421: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
422: END IF
423: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
424: CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2,
425: $ DUM1, WORK(1), -1, CHILDINFO )
426: LORGQRMIN = MAX( LORGQRMIN, M-P-1 )
427: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
428: END IF
429: IF( WANTV1T .AND. Q .GT. 0 ) THEN
430: CALL DORGLQ( Q, Q, R, V1T, LDV1T, DUM1, WORK(1), -1,
431: $ CHILDINFO )
432: LORGLQMIN = MAX( LORGLQMIN, Q )
433: LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
434: END IF
435: CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P,
436: $ THETA, DUM1, DUM2, 1, V1T, LDV1T, U2,
437: $ LDU2, U1, LDU1, DUM1, DUM1, DUM1,
438: $ DUM1, DUM1, DUM1, DUM1,
439: $ DUM1, WORK(1), -1, CHILDINFO )
440: LBBCSD = INT( WORK(1) )
441: ELSE
442: CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
443: $ DUM1, DUM1, DUM1, DUM1,
444: $ DUM1, WORK(1), -1, CHILDINFO )
445: LORBDB = M + INT( WORK(1) )
446: IF( WANTU1 .AND. P .GT. 0 ) THEN
447: CALL DORGQR( P, P, M-Q, U1, LDU1, DUM1, WORK(1), -1,
448: $ CHILDINFO )
449: LORGQRMIN = MAX( LORGQRMIN, P )
450: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
451: END IF
452: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
453: CALL DORGQR( M-P, M-P, M-Q, U2, LDU2, DUM1, WORK(1),
454: $ -1, CHILDINFO )
455: LORGQRMIN = MAX( LORGQRMIN, M-P )
456: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
457: END IF
458: IF( WANTV1T .AND. Q .GT. 0 ) THEN
459: CALL DORGLQ( Q, Q, Q, V1T, LDV1T, DUM1, WORK(1), -1,
460: $ CHILDINFO )
461: LORGLQMIN = MAX( LORGLQMIN, Q )
462: LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
463: END IF
464: CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q,
465: $ THETA, DUM1, U2, LDU2, U1, LDU1, DUM2,
466: $ 1, V1T, LDV1T, DUM1, DUM1, DUM1,
467: $ DUM1, DUM1, DUM1, DUM1,
468: $ DUM1, WORK(1), -1, CHILDINFO )
469: LBBCSD = INT( WORK(1) )
470: END IF
471: LWORKMIN = MAX( IORBDB+LORBDB-1,
472: $ IORGQR+LORGQRMIN-1,
473: $ IORGLQ+LORGLQMIN-1,
474: $ IBBCSD+LBBCSD-1 )
475: LWORKOPT = MAX( IORBDB+LORBDB-1,
476: $ IORGQR+LORGQROPT-1,
477: $ IORGLQ+LORGLQOPT-1,
478: $ IBBCSD+LBBCSD-1 )
479: WORK(1) = LWORKOPT
480: IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN
481: INFO = -19
482: END IF
483: END IF
484: IF( INFO .NE. 0 ) THEN
485: CALL XERBLA( 'DORCSD2BY1', -INFO )
486: RETURN
487: ELSE IF( LQUERY ) THEN
488: RETURN
489: END IF
490: LORGQR = LWORK-IORGQR+1
491: LORGLQ = LWORK-IORGLQ+1
492: *
493: * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
494: * in which R = MIN(P,M-P,Q,M-Q)
495: *
496: IF( R .EQ. Q ) THEN
497: *
498: * Case 1: R = Q
499: *
500: * Simultaneously bidiagonalize X11 and X21
501: *
502: CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA,
503: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
504: $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
505: *
506: * Accumulate Householder reflectors
507: *
508: IF( WANTU1 .AND. P .GT. 0 ) THEN
509: CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
510: CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
511: $ LORGQR, CHILDINFO )
512: END IF
513: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
514: CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
515: CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
516: $ WORK(IORGQR), LORGQR, CHILDINFO )
517: END IF
518: IF( WANTV1T .AND. Q .GT. 0 ) THEN
519: V1T(1,1) = ONE
520: DO J = 2, Q
521: V1T(1,J) = ZERO
522: V1T(J,1) = ZERO
523: END DO
524: CALL DLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2),
525: $ LDV1T )
526: CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, 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), LDU1 )
570: CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1),
571: $ WORK(IORGQR), LORGQR, CHILDINFO )
572: END IF
573: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
574: CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
575: CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
576: $ WORK(IORGQR), LORGQR, CHILDINFO )
577: END IF
578: IF( WANTV1T .AND. Q .GT. 0 ) THEN
579: CALL DLACPY( 'U', P, Q, X11, LDX11, V1T, LDV1T )
580: CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
581: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
582: END IF
583: *
584: * Simultaneously diagonalize X11 and X21.
585: *
586: CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
587: $ WORK(IPHI), V1T, LDV1T, DUM2, 1, U1, LDU1, U2,
588: $ LDU2, WORK(IB11D), WORK(IB11E), WORK(IB12D),
589: $ WORK(IB12E), WORK(IB21D), WORK(IB21E),
590: $ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD,
591: $ CHILDINFO )
592: *
593: * Permute rows and columns to place identity submatrices in
594: * preferred positions
595: *
596: IF( Q .GT. 0 .AND. WANTU2 ) THEN
597: DO I = 1, Q
598: IWORK(I) = M - P - Q + I
599: END DO
600: DO I = Q + 1, M - P
601: IWORK(I) = I - Q
602: END DO
603: CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
604: END IF
605: ELSE IF( R .EQ. M-P ) THEN
606: *
607: * Case 3: R = M-P
608: *
609: * Simultaneously bidiagonalize X11 and X21
610: *
611: CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
612: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
613: $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
614: *
615: * Accumulate Householder reflectors
616: *
617: IF( WANTU1 .AND. P .GT. 0 ) THEN
618: CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
619: CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
620: $ LORGQR, CHILDINFO )
621: END IF
622: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
623: U2(1,1) = ONE
624: DO J = 2, M-P
625: U2(1,J) = ZERO
626: U2(J,1) = ZERO
627: END DO
628: CALL DLACPY( 'L', M-P-1, M-P-1, X21(2,1), LDX21, U2(2,2),
629: $ LDU2 )
630: CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2,
631: $ WORK(ITAUP2), WORK(IORGQR), LORGQR, CHILDINFO )
632: END IF
633: IF( WANTV1T .AND. Q .GT. 0 ) THEN
634: CALL DLACPY( 'U', M-P, Q, X21, LDX21, V1T, LDV1T )
635: CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
636: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
637: END IF
638: *
639: * Simultaneously diagonalize X11 and X21.
640: *
641: CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P,
642: $ THETA, WORK(IPHI), DUM2, 1, V1T, LDV1T, U2,
643: $ LDU2, U1, LDU1, WORK(IB11D), WORK(IB11E),
644: $ WORK(IB12D), WORK(IB12E), WORK(IB21D),
645: $ WORK(IB21E), WORK(IB22D), WORK(IB22E),
646: $ WORK(IBBCSD), LBBCSD, CHILDINFO )
647: *
648: * Permute rows and columns to place identity submatrices in
649: * preferred positions
650: *
651: IF( Q .GT. R ) THEN
652: DO I = 1, R
653: IWORK(I) = Q - R + I
654: END DO
655: DO I = R + 1, Q
656: IWORK(I) = I - R
657: END DO
658: IF( WANTU1 ) THEN
659: CALL DLAPMT( .FALSE., P, Q, U1, LDU1, IWORK )
660: END IF
661: IF( WANTV1T ) THEN
662: CALL DLAPMR( .FALSE., Q, Q, V1T, LDV1T, IWORK )
663: END IF
664: END IF
665: ELSE
666: *
667: * Case 4: R = M-Q
668: *
669: * Simultaneously bidiagonalize X11 and X21
670: *
671: CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
672: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
673: $ WORK(ITAUQ1), WORK(IORBDB), WORK(IORBDB+M),
674: $ LORBDB-M, CHILDINFO )
675: *
676: * Accumulate Householder reflectors
677: *
678: IF( WANTU1 .AND. P .GT. 0 ) THEN
679: CALL DCOPY( P, WORK(IORBDB), 1, U1, 1 )
680: DO J = 2, P
681: U1(1,J) = ZERO
682: END DO
683: CALL DLACPY( 'L', P-1, M-Q-1, X11(2,1), LDX11, U1(2,2),
684: $ LDU1 )
685: CALL DORGQR( P, P, M-Q, U1, LDU1, WORK(ITAUP1),
686: $ WORK(IORGQR), LORGQR, CHILDINFO )
687: END IF
688: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
689: CALL DCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 )
690: DO J = 2, M-P
691: U2(1,J) = ZERO
692: END DO
693: CALL DLACPY( 'L', M-P-1, M-Q-1, X21(2,1), LDX21, U2(2,2),
694: $ LDU2 )
695: CALL DORGQR( M-P, M-P, M-Q, U2, LDU2, WORK(ITAUP2),
696: $ WORK(IORGQR), LORGQR, CHILDINFO )
697: END IF
698: IF( WANTV1T .AND. Q .GT. 0 ) THEN
699: CALL DLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T )
700: CALL DLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11,
701: $ V1T(M-Q+1,M-Q+1), LDV1T )
702: CALL DLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21,
703: $ V1T(P+1,P+1), LDV1T )
704: CALL DORGLQ( Q, Q, Q, V1T, LDV1T, WORK(ITAUQ1),
705: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
706: END IF
707: *
708: * Simultaneously diagonalize X11 and X21.
709: *
710: CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q,
711: $ THETA, WORK(IPHI), U2, LDU2, U1, LDU1, DUM2,
712: $ 1, V1T, LDV1T, WORK(IB11D), WORK(IB11E),
713: $ WORK(IB12D), WORK(IB12E), WORK(IB21D),
714: $ WORK(IB21E), WORK(IB22D), WORK(IB22E),
715: $ WORK(IBBCSD), LBBCSD, CHILDINFO )
716: *
717: * Permute rows and columns to place identity submatrices in
718: * preferred positions
719: *
720: IF( P .GT. R ) THEN
721: DO I = 1, R
722: IWORK(I) = P - R + I
723: END DO
724: DO I = R + 1, P
725: IWORK(I) = I - R
726: END DO
727: IF( WANTU1 ) THEN
728: CALL DLAPMT( .FALSE., P, P, U1, LDU1, IWORK )
729: END IF
730: IF( WANTV1T ) THEN
731: CALL DLAPMR( .FALSE., P, Q, V1T, LDV1T, IWORK )
732: END IF
733: END IF
734: END IF
735: *
736: RETURN
737: *
738: * End of DORCSD2BY1
739: *
740: END
741:
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dorcsd2by1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack