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, V1, and V2 are P-by-P,
58: *> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
59: *> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
60: *> which R = MIN(P,M-P,Q,M-Q).
61: *>
62: *>\endverbatim
63: *
64: * Arguments:
65: * ==========
66: *
67: *> \param[in] JOBU1
68: *> \verbatim
69: *> JOBU1 is CHARACTER
70: *> = 'Y': U1 is computed;
71: *> otherwise: U1 is not computed.
72: *> \endverbatim
73: *>
74: *> \param[in] JOBU2
75: *> \verbatim
76: *> JOBU2 is CHARACTER
77: *> = 'Y': U2 is computed;
78: *> otherwise: U2 is not computed.
79: *> \endverbatim
80: *>
81: *> \param[in] JOBV1T
82: *> \verbatim
83: *> JOBV1T is CHARACTER
84: *> = 'Y': V1T is computed;
85: *> otherwise: V1T is not computed.
86: *> \endverbatim
87: *>
88: *> \param[in] M
89: *> \verbatim
90: *> M is INTEGER
91: *> The number of rows and columns in X.
92: *> \endverbatim
93: *>
94: *> \param[in] P
95: *> \verbatim
96: *> P is INTEGER
97: *> The number of rows in X11 and X12. 0 <= P <= M.
98: *> \endverbatim
99: *>
100: *> \param[in] Q
101: *> \verbatim
102: *> Q is INTEGER
103: *> The number of columns in X11 and X21. 0 <= Q <= M.
104: *> \endverbatim
105: *>
106: *> \param[in,out] X11
107: *> \verbatim
108: *> X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
109: *> On entry, part of the orthogonal matrix whose CSD is
110: *> desired.
111: *> \endverbatim
112: *>
113: *> \param[in] LDX11
114: *> \verbatim
115: *> LDX11 is INTEGER
116: *> The leading dimension of X11. LDX11 >= MAX(1,P).
117: *> \endverbatim
118: *>
119: *> \param[in,out] X21
120: *> \verbatim
121: *> X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
122: *> On entry, part of the orthogonal matrix whose CSD is
123: *> desired.
124: *> \endverbatim
125: *>
126: *> \param[in] LDX21
127: *> \verbatim
128: *> LDX21 is INTEGER
129: *> The leading dimension of X21. LDX21 >= MAX(1,M-P).
130: *> \endverbatim
131: *>
132: *> \param[out] THETA
133: *> \verbatim
134: *> THETA is DOUBLE PRECISION array, dimension (R), in which R =
135: *> MIN(P,M-P,Q,M-Q).
136: *> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
137: *> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
138: *> \endverbatim
139: *>
140: *> \param[out] U1
141: *> \verbatim
142: *> U1 is DOUBLE PRECISION array, dimension (P)
143: *> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
144: *> \endverbatim
145: *>
146: *> \param[in] LDU1
147: *> \verbatim
148: *> LDU1 is INTEGER
149: *> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
150: *> MAX(1,P).
151: *> \endverbatim
152: *>
153: *> \param[out] U2
154: *> \verbatim
155: *> U2 is DOUBLE PRECISION array, dimension (M-P)
156: *> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
157: *> matrix U2.
158: *> \endverbatim
159: *>
160: *> \param[in] LDU2
161: *> \verbatim
162: *> LDU2 is INTEGER
163: *> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
164: *> MAX(1,M-P).
165: *> \endverbatim
166: *>
167: *> \param[out] V1T
168: *> \verbatim
169: *> V1T is DOUBLE PRECISION array, dimension (Q)
170: *> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
171: *> matrix V1**T.
172: *> \endverbatim
173: *>
174: *> \param[in] LDV1T
175: *> \verbatim
176: *> LDV1T is INTEGER
177: *> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
178: *> MAX(1,Q).
179: *> \endverbatim
180: *>
181: *> \param[out] WORK
182: *> \verbatim
183: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
184: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
185: *> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
186: *> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
187: *> define the matrix in intermediate bidiagonal-block form
188: *> remaining after nonconvergence. INFO specifies the number
189: *> of nonzero PHI's.
190: *> \endverbatim
191: *>
192: *> \param[in] LWORK
193: *> \verbatim
194: *> LWORK is INTEGER
195: *> The dimension of the array WORK.
196: *> \endverbatim
197: *>
198: *> If LWORK = -1, then a workspace query is assumed; the routine
199: *> only calculates the optimal size of the WORK array, returns
200: *> this value as the first entry of the work array, and no error
201: *> message related to LWORK is issued by XERBLA.
202: *> \param[out] IWORK
203: *> \verbatim
204: *> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
205: *> \endverbatim
206: *>
207: *> \param[out] INFO
208: *> \verbatim
209: *> INFO is INTEGER
210: *> = 0: successful exit.
211: *> < 0: if INFO = -i, the i-th argument had an illegal value.
212: *> > 0: DBBCSD did not converge. See the description of WORK
213: *> above for details.
214: *> \endverbatim
215: *
216: * Authors:
217: * ========
218: *
219: *> \author Univ. of Tennessee
220: *> \author Univ. of California Berkeley
221: *> \author Univ. of Colorado Denver
222: *> \author NAG Ltd.
223: *
224: *> \date July 2012
225: *
226: *> \ingroup doubleOTHERcomputational
227: *
228: *> \par References:
229: * ================
230: *>
231: *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
232: *> Algorithms, 50(1):33-65, 2009.
233: *> \endverbatim
234: *>
235: * =====================================================================
236: SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
237: $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
238: $ LDV1T, WORK, LWORK, IWORK, INFO )
239: *
240: * -- LAPACK computational routine (3.5.0) --
241: * -- LAPACK is a software package provided by Univ. of Tennessee, --
242: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
243: * July 2012
244: *
245: * .. Scalar Arguments ..
246: CHARACTER JOBU1, JOBU2, JOBV1T
247: INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
248: $ M, P, Q
249: * ..
250: * .. Array Arguments ..
251: DOUBLE PRECISION THETA(*)
252: DOUBLE PRECISION U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
253: $ X11(LDX11,*), X21(LDX21,*)
254: INTEGER IWORK(*)
255: * ..
256: *
257: * =====================================================================
258: *
259: * .. Parameters ..
260: DOUBLE PRECISION ONE, ZERO
261: PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 )
262: * ..
263: * .. Local Scalars ..
264: INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
265: $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
266: $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
267: $ J, LBBCSD, LORBDB, LORGLQ, LORGLQMIN,
268: $ LORGLQOPT, LORGQR, LORGQRMIN, LORGQROPT,
269: $ LWORKMIN, LWORKOPT, R
270: LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
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. P ) THEN
305: INFO = -13
306: ELSE IF( WANTU2 .AND. LDU2 .LT. M - P ) THEN
307: INFO = -15
308: ELSE IF( WANTV1T .AND. LDV1T .LT. 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: IF( R .EQ. Q ) THEN
351: CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0,
352: $ 0, 0, WORK, -1, CHILDINFO )
353: LORBDB = INT( WORK(1) )
354: IF( P .GE. M-P ) THEN
355: CALL DORGQR( P, P, Q, U1, LDU1, 0, WORK(1), -1,
356: $ CHILDINFO )
357: LORGQRMIN = MAX( 1, P )
358: LORGQROPT = INT( WORK(1) )
359: ELSE
360: CALL DORGQR( M-P, M-P, Q, U2, LDU2, 0, WORK(1), -1,
361: $ CHILDINFO )
362: LORGQRMIN = MAX( 1, M-P )
363: LORGQROPT = INT( WORK(1) )
364: END IF
365: CALL DORGLQ( MAX(0,Q-1), MAX(0,Q-1), MAX(0,Q-1), V1T, LDV1T,
366: $ 0, WORK(1), -1, CHILDINFO )
367: LORGLQMIN = MAX( 1, Q-1 )
368: LORGLQOPT = INT( WORK(1) )
369: CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
370: $ 0, U1, LDU1, U2, LDU2, V1T, LDV1T, 0, 1, 0, 0,
371: $ 0, 0, 0, 0, 0, 0, WORK(1), -1, CHILDINFO )
372: LBBCSD = INT( WORK(1) )
373: ELSE IF( R .EQ. P ) THEN
374: CALL DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0,
375: $ 0, 0, WORK(1), -1, CHILDINFO )
376: LORBDB = INT( WORK(1) )
377: IF( P-1 .GE. M-P ) THEN
378: CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, 0, WORK(1),
379: $ -1, CHILDINFO )
380: LORGQRMIN = MAX( 1, P-1 )
381: LORGQROPT = INT( WORK(1) )
382: ELSE
383: CALL DORGQR( M-P, M-P, Q, U2, LDU2, 0, WORK(1), -1,
384: $ CHILDINFO )
385: LORGQRMIN = MAX( 1, M-P )
386: LORGQROPT = INT( WORK(1) )
387: END IF
388: CALL DORGLQ( Q, Q, R, V1T, LDV1T, 0, WORK(1), -1,
389: $ CHILDINFO )
390: LORGLQMIN = MAX( 1, Q )
391: LORGLQOPT = INT( WORK(1) )
392: CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
393: $ 0, V1T, LDV1T, 0, 1, U1, LDU1, U2, LDU2, 0, 0,
394: $ 0, 0, 0, 0, 0, 0, WORK(1), -1, CHILDINFO )
395: LBBCSD = INT( WORK(1) )
396: ELSE IF( R .EQ. M-P ) THEN
397: CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0,
398: $ 0, 0, WORK(1), -1, CHILDINFO )
399: LORBDB = INT( WORK(1) )
400: IF( P .GE. M-P-1 ) THEN
401: CALL DORGQR( P, P, Q, U1, LDU1, 0, WORK(1), -1,
402: $ CHILDINFO )
403: LORGQRMIN = MAX( 1, P )
404: LORGQROPT = INT( WORK(1) )
405: ELSE
406: CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2, 0,
407: $ WORK(1), -1, CHILDINFO )
408: LORGQRMIN = MAX( 1, M-P-1 )
409: LORGQROPT = INT( WORK(1) )
410: END IF
411: CALL DORGLQ( Q, Q, R, V1T, LDV1T, 0, WORK(1), -1,
412: $ CHILDINFO )
413: LORGLQMIN = MAX( 1, Q )
414: LORGLQOPT = INT( WORK(1) )
415: CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P,
416: $ THETA, 0, 0, 1, V1T, LDV1T, U2, LDU2, U1, LDU1,
417: $ 0, 0, 0, 0, 0, 0, 0, 0, WORK(1), -1,
418: $ CHILDINFO )
419: LBBCSD = INT( WORK(1) )
420: ELSE
421: CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0,
422: $ 0, 0, 0, WORK(1), -1, CHILDINFO )
423: LORBDB = M + INT( WORK(1) )
424: IF( P .GE. M-P ) THEN
425: CALL DORGQR( P, P, M-Q, U1, LDU1, 0, WORK(1), -1,
426: $ CHILDINFO )
427: LORGQRMIN = MAX( 1, P )
428: LORGQROPT = INT( WORK(1) )
429: ELSE
430: CALL DORGQR( M-P, M-P, M-Q, U2, LDU2, 0, WORK(1), -1,
431: $ CHILDINFO )
432: LORGQRMIN = MAX( 1, M-P )
433: LORGQROPT = INT( WORK(1) )
434: END IF
435: CALL DORGLQ( Q, Q, Q, V1T, LDV1T, 0, WORK(1), -1,
436: $ CHILDINFO )
437: LORGLQMIN = MAX( 1, Q )
438: LORGLQOPT = INT( WORK(1) )
439: CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q,
440: $ THETA, 0, U2, LDU2, U1, LDU1, 0, 1, V1T, LDV1T,
441: $ 0, 0, 0, 0, 0, 0, 0, 0, WORK(1), -1,
442: $ CHILDINFO )
443: LBBCSD = INT( WORK(1) )
444: END IF
445: LWORKMIN = MAX( IORBDB+LORBDB-1,
446: $ IORGQR+LORGQRMIN-1,
447: $ IORGLQ+LORGLQMIN-1,
448: $ IBBCSD+LBBCSD-1 )
449: LWORKOPT = MAX( IORBDB+LORBDB-1,
450: $ IORGQR+LORGQROPT-1,
451: $ IORGLQ+LORGLQOPT-1,
452: $ IBBCSD+LBBCSD-1 )
453: WORK(1) = LWORKOPT
454: IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN
455: INFO = -19
456: END IF
457: END IF
458: IF( INFO .NE. 0 ) THEN
459: CALL XERBLA( 'DORCSD2BY1', -INFO )
460: RETURN
461: ELSE IF( LQUERY ) THEN
462: RETURN
463: END IF
464: LORGQR = LWORK-IORGQR+1
465: LORGLQ = LWORK-IORGLQ+1
466: *
467: * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
468: * in which R = MIN(P,M-P,Q,M-Q)
469: *
470: IF( R .EQ. Q ) THEN
471: *
472: * Case 1: R = Q
473: *
474: * Simultaneously bidiagonalize X11 and X21
475: *
476: CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA,
477: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
478: $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
479: *
480: * Accumulate Householder reflectors
481: *
482: IF( WANTU1 .AND. P .GT. 0 ) THEN
483: CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
484: CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
485: $ LORGQR, CHILDINFO )
486: END IF
487: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
488: CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
489: CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
490: $ WORK(IORGQR), LORGQR, CHILDINFO )
491: END IF
492: IF( WANTV1T .AND. Q .GT. 0 ) THEN
493: V1T(1,1) = ONE
494: DO J = 2, Q
495: V1T(1,J) = ZERO
496: V1T(J,1) = ZERO
497: END DO
498: CALL DLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2),
499: $ LDV1T )
500: CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
501: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
502: END IF
503: *
504: * Simultaneously diagonalize X11 and X21.
505: *
506: CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
507: $ WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, 0, 1,
508: $ WORK(IB11D), WORK(IB11E), WORK(IB12D),
509: $ WORK(IB12E), WORK(IB21D), WORK(IB21E),
510: $ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD,
511: $ CHILDINFO )
512: *
513: * Permute rows and columns to place zero submatrices in
514: * preferred positions
515: *
516: IF( Q .GT. 0 .AND. WANTU2 ) THEN
517: DO I = 1, Q
518: IWORK(I) = M - P - Q + I
519: END DO
520: DO I = Q + 1, M - P
521: IWORK(I) = I - Q
522: END DO
523: CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
524: END IF
525: ELSE IF( R .EQ. P ) THEN
526: *
527: * Case 2: R = P
528: *
529: * Simultaneously bidiagonalize X11 and X21
530: *
531: CALL DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
532: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
533: $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
534: *
535: * Accumulate Householder reflectors
536: *
537: IF( WANTU1 .AND. P .GT. 0 ) THEN
538: U1(1,1) = ONE
539: DO J = 2, P
540: U1(1,J) = ZERO
541: U1(J,1) = ZERO
542: END DO
543: CALL DLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 )
544: CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1),
545: $ WORK(IORGQR), LORGQR, CHILDINFO )
546: END IF
547: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
548: CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
549: CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
550: $ WORK(IORGQR), LORGQR, CHILDINFO )
551: END IF
552: IF( WANTV1T .AND. Q .GT. 0 ) THEN
553: CALL DLACPY( 'U', P, Q, X11, LDX11, V1T, LDV1T )
554: CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
555: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
556: END IF
557: *
558: * Simultaneously diagonalize X11 and X21.
559: *
560: CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
561: $ WORK(IPHI), V1T, LDV1T, 0, 1, U1, LDU1, U2, LDU2,
562: $ WORK(IB11D), WORK(IB11E), WORK(IB12D),
563: $ WORK(IB12E), WORK(IB21D), WORK(IB21E),
564: $ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD,
565: $ CHILDINFO )
566: *
567: * Permute rows and columns to place identity submatrices in
568: * preferred positions
569: *
570: IF( Q .GT. 0 .AND. WANTU2 ) THEN
571: DO I = 1, Q
572: IWORK(I) = M - P - Q + I
573: END DO
574: DO I = Q + 1, M - P
575: IWORK(I) = I - Q
576: END DO
577: CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
578: END IF
579: ELSE IF( R .EQ. M-P ) THEN
580: *
581: * Case 3: R = M-P
582: *
583: * Simultaneously bidiagonalize X11 and X21
584: *
585: CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
586: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
587: $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
588: *
589: * Accumulate Householder reflectors
590: *
591: IF( WANTU1 .AND. P .GT. 0 ) THEN
592: CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
593: CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
594: $ LORGQR, CHILDINFO )
595: END IF
596: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
597: U2(1,1) = ONE
598: DO J = 2, M-P
599: U2(1,J) = ZERO
600: U2(J,1) = ZERO
601: END DO
602: CALL DLACPY( 'L', M-P-1, M-P-1, X21(2,1), LDX21, U2(2,2),
603: $ LDU2 )
604: CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2,
605: $ WORK(ITAUP2), WORK(IORGQR), LORGQR, CHILDINFO )
606: END IF
607: IF( WANTV1T .AND. Q .GT. 0 ) THEN
608: CALL DLACPY( 'U', M-P, Q, X21, LDX21, V1T, LDV1T )
609: CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
610: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
611: END IF
612: *
613: * Simultaneously diagonalize X11 and X21.
614: *
615: CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P,
616: $ THETA, WORK(IPHI), 0, 1, V1T, LDV1T, U2, LDU2, U1,
617: $ LDU1, WORK(IB11D), WORK(IB11E), WORK(IB12D),
618: $ WORK(IB12E), WORK(IB21D), WORK(IB21E),
619: $ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD,
620: $ CHILDINFO )
621: *
622: * Permute rows and columns to place identity submatrices in
623: * preferred positions
624: *
625: IF( Q .GT. R ) THEN
626: DO I = 1, R
627: IWORK(I) = Q - R + I
628: END DO
629: DO I = R + 1, Q
630: IWORK(I) = I - R
631: END DO
632: IF( WANTU1 ) THEN
633: CALL DLAPMT( .FALSE., P, Q, U1, LDU1, IWORK )
634: END IF
635: IF( WANTV1T ) THEN
636: CALL DLAPMR( .FALSE., Q, Q, V1T, LDV1T, IWORK )
637: END IF
638: END IF
639: ELSE
640: *
641: * Case 4: R = M-Q
642: *
643: * Simultaneously bidiagonalize X11 and X21
644: *
645: CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
646: $ WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
647: $ WORK(ITAUQ1), WORK(IORBDB), WORK(IORBDB+M),
648: $ LORBDB-M, CHILDINFO )
649: *
650: * Accumulate Householder reflectors
651: *
652: IF( WANTU1 .AND. P .GT. 0 ) THEN
653: CALL DCOPY( P, WORK(IORBDB), 1, U1, 1 )
654: DO J = 2, P
655: U1(1,J) = ZERO
656: END DO
657: CALL DLACPY( 'L', P-1, M-Q-1, X11(2,1), LDX11, U1(2,2),
658: $ LDU1 )
659: CALL DORGQR( P, P, M-Q, U1, LDU1, WORK(ITAUP1),
660: $ WORK(IORGQR), LORGQR, CHILDINFO )
661: END IF
662: IF( WANTU2 .AND. M-P .GT. 0 ) THEN
663: CALL DCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 )
664: DO J = 2, M-P
665: U2(1,J) = ZERO
666: END DO
667: CALL DLACPY( 'L', M-P-1, M-Q-1, X21(2,1), LDX21, U2(2,2),
668: $ LDU2 )
669: CALL DORGQR( M-P, M-P, M-Q, U2, LDU2, WORK(ITAUP2),
670: $ WORK(IORGQR), LORGQR, CHILDINFO )
671: END IF
672: IF( WANTV1T .AND. Q .GT. 0 ) THEN
673: CALL DLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T )
674: CALL DLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11,
675: $ V1T(M-Q+1,M-Q+1), LDV1T )
676: CALL DLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21,
677: $ V1T(P+1,P+1), LDV1T )
678: CALL DORGLQ( Q, Q, Q, V1T, LDV1T, WORK(ITAUQ1),
679: $ WORK(IORGLQ), LORGLQ, CHILDINFO )
680: END IF
681: *
682: * Simultaneously diagonalize X11 and X21.
683: *
684: CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q,
685: $ THETA, WORK(IPHI), U2, LDU2, U1, LDU1, 0, 1, V1T,
686: $ LDV1T, WORK(IB11D), WORK(IB11E), WORK(IB12D),
687: $ WORK(IB12E), WORK(IB21D), WORK(IB21E),
688: $ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD,
689: $ CHILDINFO )
690: *
691: * Permute rows and columns to place identity submatrices in
692: * preferred positions
693: *
694: IF( P .GT. R ) THEN
695: DO I = 1, R
696: IWORK(I) = P - R + I
697: END DO
698: DO I = R + 1, P
699: IWORK(I) = I - R
700: END DO
701: IF( WANTU1 ) THEN
702: CALL DLAPMT( .FALSE., P, P, U1, LDU1, IWORK )
703: END IF
704: IF( WANTV1T ) THEN
705: CALL DLAPMR( .FALSE., P, Q, V1T, LDV1T, IWORK )
706: END IF
707: END IF
708: END IF
709: *
710: RETURN
711: *
712: * End of DORCSD2BY1
713: *
714: END
715:
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