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: LQUERY = LWORK .EQ. -1
289: *
290: IF( M .LT. 0 ) THEN
291: INFO = -4
292: ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
293: INFO = -5
294: ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
295: INFO = -6
296: ELSE IF( LDX11 .LT. MAX( 1, P ) ) THEN
297: INFO = -8
298: ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN
299: INFO = -10
300: ELSE IF( WANTU1 .AND. LDU1 .LT. MAX( 1, P ) ) THEN
301: INFO = -13
302: ELSE IF( WANTU2 .AND. LDU2 .LT. MAX( 1, M - P ) ) THEN
303: INFO = -15
304: ELSE IF( WANTV1T .AND. LDV1T .LT. MAX( 1, Q ) ) THEN
305: INFO = -17
306: END IF
307: *
308: R = MIN( P, M-P, Q, M-Q )
309: *
310: * Compute workspace
311: *
312: * WORK layout:
313: * |-------------------------------------------------------|
314: * | LWORKOPT (1) |
315: * |-------------------------------------------------------|
316: * | PHI (MAX(1,R-1)) |
317: * |-------------------------------------------------------|
318: * | TAUP1 (MAX(1,P)) | B11D (R) |
319: * | TAUP2 (MAX(1,M-P)) | B11E (R-1) |
320: * | TAUQ1 (MAX(1,Q)) | B12D (R) |
321: * |-----------------------------------------| B12E (R-1) |
322: * | DORBDB WORK | DORGQR WORK | DORGLQ WORK | B21D (R) |
323: * | | | | B21E (R-1) |
324: * | | | | B22D (R) |
325: * | | | | B22E (R-1) |
326: * | | | | DBBCSD WORK |
327: * |-------------------------------------------------------|
328: *
329: IF( INFO .EQ. 0 ) THEN
330: IPHI = 2
331: IB11D = IPHI + MAX( 1, R-1 )
332: IB11E = IB11D + MAX( 1, R )
333: IB12D = IB11E + MAX( 1, R - 1 )
334: IB12E = IB12D + MAX( 1, R )
335: IB21D = IB12E + MAX( 1, R - 1 )
336: IB21E = IB21D + MAX( 1, R )
337: IB22D = IB21E + MAX( 1, R - 1 )
338: IB22E = IB22D + MAX( 1, R )
339: IBBCSD = IB22E + MAX( 1, R - 1 )
340: ITAUP1 = IPHI + MAX( 1, R-1 )
341: ITAUP2 = ITAUP1 + MAX( 1, P )
342: ITAUQ1 = ITAUP2 + MAX( 1, M-P )
343: IORBDB = ITAUQ1 + MAX( 1, Q )
344: IORGQR = ITAUQ1 + MAX( 1, Q )
345: IORGLQ = ITAUQ1 + MAX( 1, Q )
346: LORGQRMIN = 1
347: LORGQROPT = 1
348: LORGLQMIN = 1
349: LORGLQOPT = 1
350: IF( R .EQ. Q ) THEN
351: CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA,
352: $ DUM1, DUM1, DUM1, DUM1, WORK,
353: $ -1, CHILDINFO )
354: LORBDB = INT( WORK(1) )
355: IF( WANTU1 .AND. P .GT. 0 ) THEN
356: CALL DORGQR( P, P, Q, U1, LDU1, DUM1, WORK(1), -1,
357: $ CHILDINFO )
358: LORGQRMIN = MAX( LORGQRMIN, P )
359: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
360: ENDIF
361: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
362: CALL DORGQR( M-P, M-P, Q, U2, LDU2, DUM1, WORK(1),
363: $ -1, CHILDINFO )
364: LORGQRMIN = MAX( LORGQRMIN, M-P )
365: LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
366: END IF
367: IF( WANTV1T .AND. Q .GT. 0 ) THEN
368: CALL DORGLQ( Q-1, Q-1, Q-1, V1T, LDV1T,
369: $ DUM1, WORK(1), -1, CHILDINFO )
370: LORGLQMIN = MAX( LORGLQMIN, Q-1 )
371: LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
372: END IF
373: CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, 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: ELSE IF( R .EQ. 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: IF( WANTU1 .AND. P .GT. 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: IF( WANTU2 .AND. M-P .GT. 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: IF( WANTV1T .AND. Q .GT. 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: ELSE IF( R .EQ. 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: IF( WANTU1 .AND. P .GT. 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: IF( WANTU2 .AND. M-P .GT. 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: IF( WANTV1T .AND. Q .GT. 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: IF( WANTU1 .AND. P .GT. 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: IF( WANTU2 .AND. M-P .GT. 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: IF( WANTV1T .AND. Q .GT. 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: 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( '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: IF( R .EQ. 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: IF( WANTU1 .AND. P .GT. 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: IF( WANTU2 .AND. M-P .GT. 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: 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 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: IF( Q .GT. 0 .AND. 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: ELSE IF( R .EQ. 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: IF( WANTU1 .AND. P .GT. 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: IF( WANTU2 .AND. M-P .GT. 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: IF( WANTV1T .AND. Q .GT. 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, 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: 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), 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: IF( Q .GT. 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: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
675: CALL DCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 )
676: END IF
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: 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: IF( WANTV1T .AND. Q .GT. 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, DUM1,
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: IF( P .GT. 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:
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