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