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: *> \date July 2012
228: *
229: *> \ingroup doubleOTHERcomputational
230: *
231: * =====================================================================
232: SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
233: $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
234: $ LDV1T, WORK, LWORK, IWORK, INFO )
235: *
236: * -- LAPACK computational routine (3.5.0) --
237: * -- LAPACK is a software package provided by Univ. of Tennessee, --
238: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
239: * July 2012
240: *
241: * .. Scalar Arguments ..
242: CHARACTER JOBU1, JOBU2, JOBV1T
243: INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
244: $ M, P, Q
245: * ..
246: * .. Array Arguments ..
247: DOUBLE PRECISION THETA(*)
248: DOUBLE PRECISION U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
249: $ X11(LDX11,*), X21(LDX21,*)
250: INTEGER IWORK(*)
251: * ..
252: *
253: * =====================================================================
254: *
255: * .. Parameters ..
256: DOUBLE PRECISION ONE, ZERO
257: PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 )
258: * ..
259: * .. Local Scalars ..
260: INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
261: $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
262: $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
263: $ J, LBBCSD, LORBDB, LORGLQ, LORGLQMIN,
264: $ LORGLQOPT, LORGQR, LORGQRMIN, LORGQROPT,
265: $ LWORKMIN, LWORKOPT, R
266: LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
267: * ..
268: * .. Local Arrays ..
269: DOUBLE PRECISION DUM1(1), DUM2(1,1)
270: * ..
271: * .. External Subroutines ..
272: EXTERNAL DBBCSD, DCOPY, DLACPY, DLAPMR, DLAPMT, DORBDB1,
273: $ DORBDB2, DORBDB3, DORBDB4, DORGLQ, DORGQR,
274: $ XERBLA
275: * ..
276: * .. External Functions ..
277: LOGICAL LSAME
278: EXTERNAL LSAME
279: * ..
280: * .. Intrinsic Function ..
281: INTRINSIC INT, MAX, MIN
282: * ..
283: * .. Executable Statements ..
284: *
285: * Test input arguments
286: *
287: INFO = 0
288: WANTU1 = LSAME( JOBU1, 'Y' )
289: WANTU2 = LSAME( JOBU2, 'Y' )
290: WANTV1T = LSAME( JOBV1T, 'Y' )
291: LQUERY = LWORK .EQ. -1
292: *
293: IF( M .LT. 0 ) THEN
294: INFO = -4
295: ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
296: INFO = -5
297: ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
298: INFO = -6
299: ELSE IF( LDX11 .LT. MAX( 1, P ) ) THEN
300: INFO = -8
301: ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN
302: INFO = -10
303: ELSE IF( WANTU1 .AND. LDU1 .LT. MAX( 1, P ) ) THEN
304: INFO = -13
305: ELSE IF( WANTU2 .AND. LDU2 .LT. MAX( 1, M - P ) ) THEN
306: INFO = -15
307: ELSE IF( WANTV1T .AND. LDV1T .LT. MAX( 1, Q ) ) THEN
308: INFO = -17
309: END IF
310: *
311: R = MIN( P, M-P, Q, M-Q )
312: *
313: * Compute workspace
314: *
315: * WORK layout:
316: * |-------------------------------------------------------|
317: * | LWORKOPT (1) |
318: * |-------------------------------------------------------|
319: * | PHI (MAX(1,R-1)) |
320: * |-------------------------------------------------------|
321: * | TAUP1 (MAX(1,P)) | B11D (R) |
322: * | TAUP2 (MAX(1,M-P)) | B11E (R-1) |
323: * | TAUQ1 (MAX(1,Q)) | B12D (R) |
324: * |-----------------------------------------| B12E (R-1) |
325: * | DORBDB WORK | DORGQR WORK | DORGLQ WORK | B21D (R) |
326: * | | | | B21E (R-1) |
327: * | | | | B22D (R) |
328: * | | | | B22E (R-1) |
329: * | | | | DBBCSD WORK |
330: * |-------------------------------------------------------|
331: *
332: IF( INFO .EQ. 0 ) THEN
333: IPHI = 2
334: IB11D = IPHI + MAX( 1, R-1 )
335: IB11E = IB11D + MAX( 1, R )
336: IB12D = IB11E + MAX( 1, R - 1 )
337: IB12E = IB12D + MAX( 1, R )
338: IB21D = IB12E + MAX( 1, R - 1 )
339: IB21E = IB21D + MAX( 1, R )
340: IB22D = IB21E + MAX( 1, R - 1 )
341: IB22E = IB22D + MAX( 1, R )
342: IBBCSD = IB22E + MAX( 1, R - 1 )
343: ITAUP1 = IPHI + MAX( 1, R-1 )
344: ITAUP2 = ITAUP1 + MAX( 1, P )
345: ITAUQ1 = ITAUP2 + MAX( 1, M-P )
346: IORBDB = ITAUQ1 + MAX( 1, Q )
347: IORGQR = ITAUQ1 + MAX( 1, Q )
348: IORGLQ = ITAUQ1 + MAX( 1, Q )
349: LORGQRMIN = 1
350: LORGQROPT = 1
351: LORGLQMIN = 1
352: LORGLQOPT = 1
353: IF( R .EQ. Q ) THEN
354: CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA,
355: $ DUM1, DUM1, DUM1, DUM1, WORK,
356: $ -1, CHILDINFO )
357: LORBDB = INT( WORK(1) )
358: IF( WANTU1 .AND. P .GT. 0 ) THEN
359: CALL DORGQR( P, P, Q, U1, LDU1, DUM1, WORK(1), -1,
360: $ CHILDINFO )
361: LORGQRMIN = MAX( LORGQRMIN, P )
362: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
363: ENDIF
364: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
365: CALL DORGQR( M-P, M-P, Q, U2, LDU2, DUM1, WORK(1),
366: $ -1, CHILDINFO )
367: LORGQRMIN = MAX( LORGQRMIN, M-P )
368: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
369: END IF
370: IF( WANTV1T .AND. Q .GT. 0 ) THEN
371: CALL DORGLQ( Q-1, Q-1, Q-1, V1T, LDV1T,
372: $ DUM1, WORK(1), -1, CHILDINFO )
373: LORGLQMIN = MAX( LORGLQMIN, Q-1 )
374: LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
375: END IF
376: CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
377: $ DUM1, U1, LDU1, U2, LDU2, V1T, LDV1T,
378: $ DUM2, 1, DUM1, DUM1, DUM1,
379: $ DUM1, DUM1, DUM1, DUM1,
380: $ DUM1, WORK(1), -1, CHILDINFO )
381: LBBCSD = INT( WORK(1) )
382: ELSE IF( R .EQ. P ) THEN
383: CALL DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
384: $ DUM1, DUM1, DUM1, DUM1,
385: $ WORK(1), -1, CHILDINFO )
386: LORBDB = INT( WORK(1) )
387: IF( WANTU1 .AND. P .GT. 0 ) THEN
388: CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, DUM1,
389: $ WORK(1), -1, CHILDINFO )
390: LORGQRMIN = MAX( LORGQRMIN, P-1 )
391: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
392: END IF
393: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
394: CALL DORGQR( M-P, M-P, Q, U2, LDU2, DUM1, WORK(1),
395: $ -1, CHILDINFO )
396: LORGQRMIN = MAX( LORGQRMIN, M-P )
397: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
398: END IF
399: IF( WANTV1T .AND. Q .GT. 0 ) THEN
400: CALL DORGLQ( Q, Q, R, V1T, LDV1T, DUM1, WORK(1), -1,
401: $ CHILDINFO )
402: LORGLQMIN = MAX( LORGLQMIN, Q )
403: LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
404: END IF
405: CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
406: $ DUM1, V1T, LDV1T, DUM2, 1, U1, LDU1,
407: $ U2, LDU2, DUM1, DUM1, DUM1,
408: $ DUM1, DUM1, DUM1, DUM1,
409: $ DUM1, WORK(1), -1, CHILDINFO )
410: LBBCSD = INT( WORK(1) )
411: ELSE IF( R .EQ. M-P ) THEN
412: CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
413: $ DUM1, DUM1, DUM1, DUM1,
414: $ WORK(1), -1, CHILDINFO )
415: LORBDB = INT( WORK(1) )
416: IF( WANTU1 .AND. P .GT. 0 ) THEN
417: CALL DORGQR( P, P, Q, U1, LDU1, DUM1, WORK(1), -1,
418: $ CHILDINFO )
419: LORGQRMIN = MAX( LORGQRMIN, P )
420: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
421: END IF
422: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
423: CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2,
424: $ DUM1, WORK(1), -1, CHILDINFO )
425: LORGQRMIN = MAX( LORGQRMIN, M-P-1 )
426: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
427: END IF
428: IF( WANTV1T .AND. Q .GT. 0 ) THEN
429: CALL DORGLQ( Q, Q, R, V1T, LDV1T, DUM1, WORK(1), -1,
430: $ CHILDINFO )
431: LORGLQMIN = MAX( LORGLQMIN, Q )
432: LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
433: END IF
434: CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P,
435: $ THETA, DUM1, DUM2, 1, V1T, LDV1T, U2,
436: $ LDU2, U1, LDU1, DUM1, DUM1, DUM1,
437: $ DUM1, DUM1, DUM1, DUM1,
438: $ DUM1, WORK(1), -1, CHILDINFO )
439: LBBCSD = INT( WORK(1) )
440: ELSE
441: CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
442: $ DUM1, DUM1, DUM1, DUM1,
443: $ DUM1, WORK(1), -1, CHILDINFO )
444: LORBDB = M + INT( WORK(1) )
445: IF( WANTU1 .AND. P .GT. 0 ) THEN
446: CALL DORGQR( P, P, M-Q, U1, LDU1, DUM1, WORK(1), -1,
447: $ CHILDINFO )
448: LORGQRMIN = MAX( LORGQRMIN, P )
449: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
450: END IF
451: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
452: CALL DORGQR( M-P, M-P, M-Q, U2, LDU2, DUM1, WORK(1),
453: $ -1, CHILDINFO )
454: LORGQRMIN = MAX( LORGQRMIN, M-P )
455: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
456: END IF
457: IF( WANTV1T .AND. Q .GT. 0 ) THEN
458: CALL DORGLQ( Q, Q, Q, V1T, LDV1T, DUM1, WORK(1), -1,
459: $ CHILDINFO )
460: LORGLQMIN = MAX( LORGLQMIN, Q )
461: LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
462: END IF
463: CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q,
464: $ THETA, DUM1, U2, LDU2, U1, LDU1, DUM2,
465: $ 1, V1T, LDV1T, DUM1, DUM1, DUM1,
466: $ DUM1, DUM1, DUM1, DUM1,
467: $ DUM1, WORK(1), -1, CHILDINFO )
468: LBBCSD = INT( WORK(1) )
469: END IF
470: LWORKMIN = MAX( IORBDB+LORBDB-1,
471: $ IORGQR+LORGQRMIN-1,
472: $ IORGLQ+LORGLQMIN-1,
473: $ IBBCSD+LBBCSD-1 )
474: LWORKOPT = MAX( IORBDB+LORBDB-1,
475: $ IORGQR+LORGQROPT-1,
476: $ IORGLQ+LORGLQOPT-1,
477: $ IBBCSD+LBBCSD-1 )
478: WORK(1) = LWORKOPT
479: IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN
480: INFO = -19
481: END IF
482: END IF
483: IF( INFO .NE. 0 ) THEN
484: CALL XERBLA( 'DORCSD2BY1', -INFO )
485: RETURN
486: ELSE IF( LQUERY ) THEN
487: RETURN
488: END IF
489: LORGQR = LWORK-IORGQR+1
490: LORGLQ = LWORK-IORGLQ+1
491: *
492: * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
493: * in which R = MIN(P,M-P,Q,M-Q)
494: *
495: IF( R .EQ. Q ) THEN
496: *
497: * Case 1: R = Q
498: *
499: * Simultaneously bidiagonalize X11 and X21
500: *
501: CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA,
502: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
503: $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
504: *
505: * Accumulate Householder reflectors
506: *
507: IF( WANTU1 .AND. P .GT. 0 ) THEN
508: CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
509: CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), 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, WORK(ITAUQ1),
526: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
527: END IF
528: *
529: * Simultaneously diagonalize X11 and X21.
530: *
531: CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
532: $ WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T,
533: $ DUM2, 1, WORK(IB11D), WORK(IB11E),
534: $ WORK(IB12D), WORK(IB12E), WORK(IB21D),
535: $ WORK(IB21E), WORK(IB22D), WORK(IB22E),
536: $ WORK(IBBCSD), LBBCSD, CHILDINFO )
537: *
538: * Permute rows and columns to place zero submatrices in
539: * preferred positions
540: *
541: IF( Q .GT. 0 .AND. WANTU2 ) THEN
542: DO I = 1, Q
543: IWORK(I) = M - P - Q + I
544: END DO
545: DO I = Q + 1, M - P
546: IWORK(I) = I - Q
547: END DO
548: CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
549: END IF
550: ELSE IF( R .EQ. P ) THEN
551: *
552: * Case 2: R = P
553: *
554: * Simultaneously bidiagonalize X11 and X21
555: *
556: CALL DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
557: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
558: $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
559: *
560: * Accumulate Householder reflectors
561: *
562: IF( WANTU1 .AND. P .GT. 0 ) THEN
563: U1(1,1) = ONE
564: DO J = 2, P
565: U1(1,J) = ZERO
566: U1(J,1) = ZERO
567: END DO
568: CALL DLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 )
569: CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1),
570: $ WORK(IORGQR), LORGQR, CHILDINFO )
571: END IF
572: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
573: CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
574: CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
575: $ WORK(IORGQR), LORGQR, CHILDINFO )
576: END IF
577: IF( WANTV1T .AND. Q .GT. 0 ) THEN
578: CALL DLACPY( 'U', P, Q, X11, LDX11, V1T, LDV1T )
579: CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
580: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
581: END IF
582: *
583: * Simultaneously diagonalize X11 and X21.
584: *
585: CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
586: $ WORK(IPHI), V1T, LDV1T, DUM2, 1, U1, LDU1, U2,
587: $ LDU2, WORK(IB11D), WORK(IB11E), WORK(IB12D),
588: $ WORK(IB12E), WORK(IB21D), WORK(IB21E),
589: $ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD,
590: $ CHILDINFO )
591: *
592: * Permute rows and columns to place identity submatrices in
593: * preferred positions
594: *
595: IF( Q .GT. 0 .AND. WANTU2 ) THEN
596: DO I = 1, Q
597: IWORK(I) = M - P - Q + I
598: END DO
599: DO I = Q + 1, M - P
600: IWORK(I) = I - Q
601: END DO
602: CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
603: END IF
604: ELSE IF( R .EQ. M-P ) THEN
605: *
606: * Case 3: R = M-P
607: *
608: * Simultaneously bidiagonalize X11 and X21
609: *
610: CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
611: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
612: $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
613: *
614: * Accumulate Householder reflectors
615: *
616: IF( WANTU1 .AND. P .GT. 0 ) THEN
617: CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
618: CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
619: $ LORGQR, CHILDINFO )
620: END IF
621: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
622: U2(1,1) = ONE
623: DO J = 2, M-P
624: U2(1,J) = ZERO
625: U2(J,1) = ZERO
626: END DO
627: CALL DLACPY( 'L', M-P-1, M-P-1, X21(2,1), LDX21, U2(2,2),
628: $ LDU2 )
629: CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2,
630: $ WORK(ITAUP2), WORK(IORGQR), LORGQR, CHILDINFO )
631: END IF
632: IF( WANTV1T .AND. Q .GT. 0 ) THEN
633: CALL DLACPY( 'U', M-P, Q, X21, LDX21, V1T, LDV1T )
634: CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
635: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
636: END IF
637: *
638: * Simultaneously diagonalize X11 and X21.
639: *
640: CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P,
641: $ THETA, WORK(IPHI), DUM2, 1, V1T, LDV1T, U2,
642: $ LDU2, U1, LDU1, WORK(IB11D), WORK(IB11E),
643: $ WORK(IB12D), WORK(IB12E), WORK(IB21D),
644: $ WORK(IB21E), WORK(IB22D), WORK(IB22E),
645: $ WORK(IBBCSD), LBBCSD, CHILDINFO )
646: *
647: * Permute rows and columns to place identity submatrices in
648: * preferred positions
649: *
650: IF( Q .GT. R ) THEN
651: DO I = 1, R
652: IWORK(I) = Q - R + I
653: END DO
654: DO I = R + 1, Q
655: IWORK(I) = I - R
656: END DO
657: IF( WANTU1 ) THEN
658: CALL DLAPMT( .FALSE., P, Q, U1, LDU1, IWORK )
659: END IF
660: IF( WANTV1T ) THEN
661: CALL DLAPMR( .FALSE., Q, Q, V1T, LDV1T, IWORK )
662: END IF
663: END IF
664: ELSE
665: *
666: * Case 4: R = M-Q
667: *
668: * Simultaneously bidiagonalize X11 and X21
669: *
670: CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
671: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
672: $ WORK(ITAUQ1), WORK(IORBDB), WORK(IORBDB+M),
673: $ LORBDB-M, CHILDINFO )
674: *
675: * Accumulate Householder reflectors
676: *
677: IF( WANTU1 .AND. P .GT. 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: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
688: CALL DCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 )
689: DO J = 2, M-P
690: U2(1,J) = ZERO
691: END DO
692: CALL DLACPY( 'L', M-P-1, M-Q-1, X21(2,1), LDX21, U2(2,2),
693: $ LDU2 )
694: CALL DORGQR( M-P, M-P, M-Q, U2, LDU2, WORK(ITAUP2),
695: $ WORK(IORGQR), LORGQR, CHILDINFO )
696: END IF
697: IF( WANTV1T .AND. Q .GT. 0 ) THEN
698: CALL DLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T )
699: CALL DLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11,
700: $ V1T(M-Q+1,M-Q+1), LDV1T )
701: CALL DLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21,
702: $ V1T(P+1,P+1), LDV1T )
703: CALL DORGLQ( Q, Q, Q, V1T, LDV1T, WORK(ITAUQ1),
704: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
705: END IF
706: *
707: * Simultaneously diagonalize X11 and X21.
708: *
709: CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q,
710: $ THETA, WORK(IPHI), U2, LDU2, U1, LDU1, DUM2,
711: $ 1, V1T, LDV1T, WORK(IB11D), WORK(IB11E),
712: $ WORK(IB12D), WORK(IB12E), WORK(IB21D),
713: $ WORK(IB21E), WORK(IB22D), WORK(IB22E),
714: $ WORK(IBBCSD), LBBCSD, CHILDINFO )
715: *
716: * Permute rows and columns to place identity submatrices in
717: * preferred positions
718: *
719: IF( P .GT. R ) THEN
720: DO I = 1, R
721: IWORK(I) = P - R + I
722: END DO
723: DO I = R + 1, P
724: IWORK(I) = I - R
725: END DO
726: IF( WANTU1 ) THEN
727: CALL DLAPMT( .FALSE., P, P, U1, LDU1, IWORK )
728: END IF
729: IF( WANTV1T ) THEN
730: CALL DLAPMR( .FALSE., P, Q, V1T, LDV1T, IWORK )
731: END IF
732: END IF
733: END IF
734: *
735: RETURN
736: *
737: * End of DORCSD2BY1
738: *
739: END
740:
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