1: SUBROUTINE ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
2: $ THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T,
3: $ V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E,
4: $ B22D, B22E, RWORK, LRWORK, INFO )
5: IMPLICIT NONE
6: *
7: * -- LAPACK routine (version 3.3.0) --
8: *
9: * -- Contributed by Brian Sutton of the Randolph-Macon College --
10: * -- November 2010
11: *
12: * -- LAPACK is a software package provided by Univ. of Tennessee, --
13: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
14: *
15: * .. Scalar Arguments ..
16: CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
17: INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q
18: * ..
19: * .. Array Arguments ..
20: DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
21: $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
22: $ PHI( * ), THETA( * ), RWORK( * )
23: COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
24: $ V2T( LDV2T, * )
25: * ..
26: *
27: * Purpose
28: * =======
29: *
30: * ZBBCSD computes the CS decomposition of a unitary matrix in
31: * bidiagonal-block form,
32: *
33: *
34: * [ B11 | B12 0 0 ]
35: * [ 0 | 0 -I 0 ]
36: * X = [----------------]
37: * [ B21 | B22 0 0 ]
38: * [ 0 | 0 0 I ]
39: *
40: * [ C | -S 0 0 ]
41: * [ U1 | ] [ 0 | 0 -I 0 ] [ V1 | ]**H
42: * = [---------] [---------------] [---------] .
43: * [ | U2 ] [ S | C 0 0 ] [ | V2 ]
44: * [ 0 | 0 0 I ]
45: *
46: * X is M-by-M, its top-left block is P-by-Q, and Q must be no larger
47: * than P, M-P, or M-Q. (If Q is not the smallest index, then X must be
48: * transposed and/or permuted. This can be done in constant time using
49: * the TRANS and SIGNS options. See ZUNCSD for details.)
50: *
51: * The bidiagonal matrices B11, B12, B21, and B22 are represented
52: * implicitly by angles THETA(1:Q) and PHI(1:Q-1).
53: *
54: * The unitary matrices U1, U2, V1T, and V2T are input/output.
55: * The input matrices are pre- or post-multiplied by the appropriate
56: * singular vector matrices.
57: *
58: * Arguments
59: * =========
60: *
61: * JOBU1 (input) CHARACTER
62: * = 'Y': U1 is updated;
63: * otherwise: U1 is not updated.
64: *
65: * JOBU2 (input) CHARACTER
66: * = 'Y': U2 is updated;
67: * otherwise: U2 is not updated.
68: *
69: * JOBV1T (input) CHARACTER
70: * = 'Y': V1T is updated;
71: * otherwise: V1T is not updated.
72: *
73: * JOBV2T (input) CHARACTER
74: * = 'Y': V2T is updated;
75: * otherwise: V2T is not updated.
76: *
77: * TRANS (input) CHARACTER
78: * = 'T': X, U1, U2, V1T, and V2T are stored in row-major
79: * order;
80: * otherwise: X, U1, U2, V1T, and V2T are stored in column-
81: * major order.
82: *
83: * M (input) INTEGER
84: * The number of rows and columns in X, the unitary matrix in
85: * bidiagonal-block form.
86: *
87: * P (input) INTEGER
88: * The number of rows in the top-left block of X. 0 <= P <= M.
89: *
90: * Q (input) INTEGER
91: * The number of columns in the top-left block of X.
92: * 0 <= Q <= MIN(P,M-P,M-Q).
93: *
94: * THETA (input/output) DOUBLE PRECISION array, dimension (Q)
95: * On entry, the angles THETA(1),...,THETA(Q) that, along with
96: * PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block
97: * form. On exit, the angles whose cosines and sines define the
98: * diagonal blocks in the CS decomposition.
99: *
100: * PHI (input/workspace) DOUBLE PRECISION array, dimension (Q-1)
101: * The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),...,
102: * THETA(Q), define the matrix in bidiagonal-block form.
103: *
104: * U1 (input/output) COMPLEX*16 array, dimension (LDU1,P)
105: * On entry, an LDU1-by-P matrix. On exit, U1 is postmultiplied
106: * by the left singular vector matrix common to [ B11 ; 0 ] and
107: * [ B12 0 0 ; 0 -I 0 0 ].
108: *
109: * LDU1 (input) INTEGER
110: * The leading dimension of the array U1.
111: *
112: * U2 (input/output) COMPLEX*16 array, dimension (LDU2,M-P)
113: * On entry, an LDU2-by-(M-P) matrix. On exit, U2 is
114: * postmultiplied by the left singular vector matrix common to
115: * [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ].
116: *
117: * LDU2 (input) INTEGER
118: * The leading dimension of the array U2.
119: *
120: * V1T (input/output) COMPLEX*16 array, dimension (LDV1T,Q)
121: * On entry, a LDV1T-by-Q matrix. On exit, V1T is premultiplied
122: * by the conjugate transpose of the right singular vector
123: * matrix common to [ B11 ; 0 ] and [ B21 ; 0 ].
124: *
125: * LDV1T (input) INTEGER
126: * The leading dimension of the array V1T.
127: *
128: * V2T (input/output) COMPLEX*16 array, dimenison (LDV2T,M-Q)
129: * On entry, a LDV2T-by-(M-Q) matrix. On exit, V2T is
130: * premultiplied by the conjugate transpose of the right
131: * singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and
132: * [ B22 0 0 ; 0 0 I ].
133: *
134: * LDV2T (input) INTEGER
135: * The leading dimension of the array V2T.
136: *
137: * B11D (output) DOUBLE PRECISION array, dimension (Q)
138: * When ZBBCSD converges, B11D contains the cosines of THETA(1),
139: * ..., THETA(Q). If ZBBCSD fails to converge, then B11D
140: * contains the diagonal of the partially reduced top-left
141: * block.
142: *
143: * B11E (output) DOUBLE PRECISION array, dimension (Q-1)
144: * When ZBBCSD converges, B11E contains zeros. If ZBBCSD fails
145: * to converge, then B11E contains the superdiagonal of the
146: * partially reduced top-left block.
147: *
148: * B12D (output) DOUBLE PRECISION array, dimension (Q)
149: * When ZBBCSD converges, B12D contains the negative sines of
150: * THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
151: * B12D contains the diagonal of the partially reduced top-right
152: * block.
153: *
154: * B12E (output) DOUBLE PRECISION array, dimension (Q-1)
155: * When ZBBCSD converges, B12E contains zeros. If ZBBCSD fails
156: * to converge, then B12E contains the subdiagonal of the
157: * partially reduced top-right block.
158: *
159: * RWORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
160: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
161: *
162: * LRWORK (input) INTEGER
163: * The dimension of the array RWORK. LRWORK >= MAX(1,8*Q).
164: *
165: * If LRWORK = -1, then a workspace query is assumed; the
166: * routine only calculates the optimal size of the RWORK array,
167: * returns this value as the first entry of the work array, and
168: * no error message related to LRWORK is issued by XERBLA.
169: *
170: * INFO (output) INTEGER
171: * = 0: successful exit.
172: * < 0: if INFO = -i, the i-th argument had an illegal value.
173: * > 0: if ZBBCSD did not converge, INFO specifies the number
174: * of nonzero entries in PHI, and B11D, B11E, etc.,
175: * contain the partially reduced matrix.
176: *
177: * Reference
178: * =========
179: *
180: * [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
181: * Algorithms, 50(1):33-65, 2009.
182: *
183: * Internal Parameters
184: * ===================
185: *
186: * TOLMUL DOUBLE PRECISION, default = MAX(10,MIN(100,EPS**(-1/8)))
187: * TOLMUL controls the convergence criterion of the QR loop.
188: * Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they
189: * are within TOLMUL*EPS of either bound.
190: *
191: * ===================================================================
192: *
193: * .. Parameters ..
194: INTEGER MAXITR
195: PARAMETER ( MAXITR = 6 )
196: DOUBLE PRECISION HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO
197: PARAMETER ( HUNDRED = 100.0D0, MEIGHTH = -0.125D0,
198: $ ONE = 1.0D0, PIOVER2 = 1.57079632679489662D0,
199: $ TEN = 10.0D0, ZERO = 0.0D0 )
200: COMPLEX*16 NEGONECOMPLEX
201: PARAMETER ( NEGONECOMPLEX = (-1.0D0,0.0D0) )
202: * ..
203: * .. Local Scalars ..
204: LOGICAL COLMAJOR, LQUERY, RESTART11, RESTART12,
205: $ RESTART21, RESTART22, WANTU1, WANTU2, WANTV1T,
206: $ WANTV2T
207: INTEGER I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS,
208: $ IU2SN, IV1TCS, IV1TSN, IV2TCS, IV2TSN, J,
209: $ LRWORKMIN, LRWORKOPT, MAXIT, MINI
210: DOUBLE PRECISION B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY,
211: $ EPS, MU, NU, R, SIGMA11, SIGMA21,
212: $ TEMP, THETAMAX, THETAMIN, THRESH, TOL, TOLMUL,
213: $ UNFL, X1, X2, Y1, Y2
214: *
215: EXTERNAL DLARTGP, DLARTGS, DLAS2, XERBLA, ZLASR, ZSCAL,
216: $ ZSWAP
217: * ..
218: * .. External Functions ..
219: DOUBLE PRECISION DLAMCH
220: LOGICAL LSAME
221: EXTERNAL LSAME, DLAMCH
222: * ..
223: * .. Intrinsic Functions ..
224: INTRINSIC ABS, ATAN2, COS, MAX, MIN, SIN, SQRT
225: * ..
226: * .. Executable Statements ..
227: *
228: * Test input arguments
229: *
230: INFO = 0
231: LQUERY = LRWORK .EQ. -1
232: WANTU1 = LSAME( JOBU1, 'Y' )
233: WANTU2 = LSAME( JOBU2, 'Y' )
234: WANTV1T = LSAME( JOBV1T, 'Y' )
235: WANTV2T = LSAME( JOBV2T, 'Y' )
236: COLMAJOR = .NOT. LSAME( TRANS, 'T' )
237: *
238: IF( M .LT. 0 ) THEN
239: INFO = -6
240: ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
241: INFO = -7
242: ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
243: INFO = -8
244: ELSE IF( Q .GT. P .OR. Q .GT. M-P .OR. Q .GT. M-Q ) THEN
245: INFO = -8
246: ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
247: INFO = -12
248: ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
249: INFO = -14
250: ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
251: INFO = -16
252: ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
253: INFO = -18
254: END IF
255: *
256: * Quick return if Q = 0
257: *
258: IF( INFO .EQ. 0 .AND. Q .EQ. 0 ) THEN
259: LRWORKMIN = 1
260: RWORK(1) = LRWORKMIN
261: RETURN
262: END IF
263: *
264: * Compute workspace
265: *
266: IF( INFO .EQ. 0 ) THEN
267: IU1CS = 1
268: IU1SN = IU1CS + Q
269: IU2CS = IU1SN + Q
270: IU2SN = IU2CS + Q
271: IV1TCS = IU2SN + Q
272: IV1TSN = IV1TCS + Q
273: IV2TCS = IV1TSN + Q
274: IV2TSN = IV2TCS + Q
275: LRWORKOPT = IV2TSN + Q - 1
276: LRWORKMIN = LRWORKOPT
277: RWORK(1) = LRWORKOPT
278: IF( LRWORK .LT. LRWORKMIN .AND. .NOT. LQUERY ) THEN
279: INFO = -28
280: END IF
281: END IF
282: *
283: IF( INFO .NE. 0 ) THEN
284: CALL XERBLA( 'ZBBCSD', -INFO )
285: RETURN
286: ELSE IF( LQUERY ) THEN
287: RETURN
288: END IF
289: *
290: * Get machine constants
291: *
292: EPS = DLAMCH( 'Epsilon' )
293: UNFL = DLAMCH( 'Safe minimum' )
294: TOLMUL = MAX( TEN, MIN( HUNDRED, EPS**MEIGHTH ) )
295: TOL = TOLMUL*EPS
296: THRESH = MAX( TOL, MAXITR*Q*Q*UNFL )
297: *
298: * Test for negligible sines or cosines
299: *
300: DO I = 1, Q
301: IF( THETA(I) .LT. THRESH ) THEN
302: THETA(I) = ZERO
303: ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
304: THETA(I) = PIOVER2
305: END IF
306: END DO
307: DO I = 1, Q-1
308: IF( PHI(I) .LT. THRESH ) THEN
309: PHI(I) = ZERO
310: ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
311: PHI(I) = PIOVER2
312: END IF
313: END DO
314: *
315: * Initial deflation
316: *
317: IMAX = Q
318: DO WHILE( ( IMAX .GT. 1 ) .AND. ( PHI(IMAX-1) .EQ. ZERO ) )
319: IMAX = IMAX - 1
320: END DO
321: IMIN = IMAX - 1
322: IF ( IMIN .GT. 1 ) THEN
323: DO WHILE( PHI(IMIN-1) .NE. ZERO )
324: IMIN = IMIN - 1
325: IF ( IMIN .LE. 1 ) EXIT
326: END DO
327: END IF
328: *
329: * Initialize iteration counter
330: *
331: MAXIT = MAXITR*Q*Q
332: ITER = 0
333: *
334: * Begin main iteration loop
335: *
336: DO WHILE( IMAX .GT. 1 )
337: *
338: * Compute the matrix entries
339: *
340: B11D(IMIN) = COS( THETA(IMIN) )
341: B21D(IMIN) = -SIN( THETA(IMIN) )
342: DO I = IMIN, IMAX - 1
343: B11E(I) = -SIN( THETA(I) ) * SIN( PHI(I) )
344: B11D(I+1) = COS( THETA(I+1) ) * COS( PHI(I) )
345: B12D(I) = SIN( THETA(I) ) * COS( PHI(I) )
346: B12E(I) = COS( THETA(I+1) ) * SIN( PHI(I) )
347: B21E(I) = -COS( THETA(I) ) * SIN( PHI(I) )
348: B21D(I+1) = -SIN( THETA(I+1) ) * COS( PHI(I) )
349: B22D(I) = COS( THETA(I) ) * COS( PHI(I) )
350: B22E(I) = -SIN( THETA(I+1) ) * SIN( PHI(I) )
351: END DO
352: B12D(IMAX) = SIN( THETA(IMAX) )
353: B22D(IMAX) = COS( THETA(IMAX) )
354: *
355: * Abort if not converging; otherwise, increment ITER
356: *
357: IF( ITER .GT. MAXIT ) THEN
358: INFO = 0
359: DO I = 1, Q
360: IF( PHI(I) .NE. ZERO )
361: $ INFO = INFO + 1
362: END DO
363: RETURN
364: END IF
365: *
366: ITER = ITER + IMAX - IMIN
367: *
368: * Compute shifts
369: *
370: THETAMAX = THETA(IMIN)
371: THETAMIN = THETA(IMIN)
372: DO I = IMIN+1, IMAX
373: IF( THETA(I) > THETAMAX )
374: $ THETAMAX = THETA(I)
375: IF( THETA(I) < THETAMIN )
376: $ THETAMIN = THETA(I)
377: END DO
378: *
379: IF( THETAMAX .GT. PIOVER2 - THRESH ) THEN
380: *
381: * Zero on diagonals of B11 and B22; induce deflation with a
382: * zero shift
383: *
384: MU = ZERO
385: NU = ONE
386: *
387: ELSE IF( THETAMIN .LT. THRESH ) THEN
388: *
389: * Zero on diagonals of B12 and B22; induce deflation with a
390: * zero shift
391: *
392: MU = ONE
393: NU = ZERO
394: *
395: ELSE
396: *
397: * Compute shifts for B11 and B21 and use the lesser
398: *
399: CALL DLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11,
400: $ DUMMY )
401: CALL DLAS2( B21D(IMAX-1), B21E(IMAX-1), B21D(IMAX), SIGMA21,
402: $ DUMMY )
403: *
404: IF( SIGMA11 .LE. SIGMA21 ) THEN
405: MU = SIGMA11
406: NU = SQRT( ONE - MU**2 )
407: IF( MU .LT. THRESH ) THEN
408: MU = ZERO
409: NU = ONE
410: END IF
411: ELSE
412: NU = SIGMA21
413: MU = SQRT( 1.0 - NU**2 )
414: IF( NU .LT. THRESH ) THEN
415: MU = ONE
416: NU = ZERO
417: END IF
418: END IF
419: END IF
420: *
421: * Rotate to produce bulges in B11 and B21
422: *
423: IF( MU .LE. NU ) THEN
424: CALL DLARTGS( B11D(IMIN), B11E(IMIN), MU,
425: $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) )
426: ELSE
427: CALL DLARTGS( B21D(IMIN), B21E(IMIN), NU,
428: $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) )
429: END IF
430: *
431: TEMP = RWORK(IV1TCS+IMIN-1)*B11D(IMIN) +
432: $ RWORK(IV1TSN+IMIN-1)*B11E(IMIN)
433: B11E(IMIN) = RWORK(IV1TCS+IMIN-1)*B11E(IMIN) -
434: $ RWORK(IV1TSN+IMIN-1)*B11D(IMIN)
435: B11D(IMIN) = TEMP
436: B11BULGE = RWORK(IV1TSN+IMIN-1)*B11D(IMIN+1)
437: B11D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B11D(IMIN+1)
438: TEMP = RWORK(IV1TCS+IMIN-1)*B21D(IMIN) +
439: $ RWORK(IV1TSN+IMIN-1)*B21E(IMIN)
440: B21E(IMIN) = RWORK(IV1TCS+IMIN-1)*B21E(IMIN) -
441: $ RWORK(IV1TSN+IMIN-1)*B21D(IMIN)
442: B21D(IMIN) = TEMP
443: B21BULGE = RWORK(IV1TSN+IMIN-1)*B21D(IMIN+1)
444: B21D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B21D(IMIN+1)
445: *
446: * Compute THETA(IMIN)
447: *
448: THETA( IMIN ) = ATAN2( SQRT( B21D(IMIN)**2+B21BULGE**2 ),
449: $ SQRT( B11D(IMIN)**2+B11BULGE**2 ) )
450: *
451: * Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
452: *
453: IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN
454: CALL DLARTGP( B11BULGE, B11D(IMIN), RWORK(IU1SN+IMIN-1),
455: $ RWORK(IU1CS+IMIN-1), R )
456: ELSE IF( MU .LE. NU ) THEN
457: CALL DLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU,
458: $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) )
459: ELSE
460: CALL DLARTGS( B12D( IMIN ), B12E( IMIN ), NU,
461: $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) )
462: END IF
463: IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN
464: CALL DLARTGP( B21BULGE, B21D(IMIN), RWORK(IU2SN+IMIN-1),
465: $ RWORK(IU2CS+IMIN-1), R )
466: ELSE IF( NU .LT. MU ) THEN
467: CALL DLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU,
468: $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) )
469: ELSE
470: CALL DLARTGS( B22D(IMIN), B22E(IMIN), MU,
471: $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) )
472: END IF
473: RWORK(IU2CS+IMIN-1) = -RWORK(IU2CS+IMIN-1)
474: RWORK(IU2SN+IMIN-1) = -RWORK(IU2SN+IMIN-1)
475: *
476: TEMP = RWORK(IU1CS+IMIN-1)*B11E(IMIN) +
477: $ RWORK(IU1SN+IMIN-1)*B11D(IMIN+1)
478: B11D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11D(IMIN+1) -
479: $ RWORK(IU1SN+IMIN-1)*B11E(IMIN)
480: B11E(IMIN) = TEMP
481: IF( IMAX .GT. IMIN+1 ) THEN
482: B11BULGE = RWORK(IU1SN+IMIN-1)*B11E(IMIN+1)
483: B11E(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11E(IMIN+1)
484: END IF
485: TEMP = RWORK(IU1CS+IMIN-1)*B12D(IMIN) +
486: $ RWORK(IU1SN+IMIN-1)*B12E(IMIN)
487: B12E(IMIN) = RWORK(IU1CS+IMIN-1)*B12E(IMIN) -
488: $ RWORK(IU1SN+IMIN-1)*B12D(IMIN)
489: B12D(IMIN) = TEMP
490: B12BULGE = RWORK(IU1SN+IMIN-1)*B12D(IMIN+1)
491: B12D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B12D(IMIN+1)
492: TEMP = RWORK(IU2CS+IMIN-1)*B21E(IMIN) +
493: $ RWORK(IU2SN+IMIN-1)*B21D(IMIN+1)
494: B21D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21D(IMIN+1) -
495: $ RWORK(IU2SN+IMIN-1)*B21E(IMIN)
496: B21E(IMIN) = TEMP
497: IF( IMAX .GT. IMIN+1 ) THEN
498: B21BULGE = RWORK(IU2SN+IMIN-1)*B21E(IMIN+1)
499: B21E(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21E(IMIN+1)
500: END IF
501: TEMP = RWORK(IU2CS+IMIN-1)*B22D(IMIN) +
502: $ RWORK(IU2SN+IMIN-1)*B22E(IMIN)
503: B22E(IMIN) = RWORK(IU2CS+IMIN-1)*B22E(IMIN) -
504: $ RWORK(IU2SN+IMIN-1)*B22D(IMIN)
505: B22D(IMIN) = TEMP
506: B22BULGE = RWORK(IU2SN+IMIN-1)*B22D(IMIN+1)
507: B22D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B22D(IMIN+1)
508: *
509: * Inner loop: chase bulges from B11(IMIN,IMIN+2),
510: * B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
511: * bottom-right
512: *
513: DO I = IMIN+1, IMAX-1
514: *
515: * Compute PHI(I-1)
516: *
517: X1 = SIN(THETA(I-1))*B11E(I-1) + COS(THETA(I-1))*B21E(I-1)
518: X2 = SIN(THETA(I-1))*B11BULGE + COS(THETA(I-1))*B21BULGE
519: Y1 = SIN(THETA(I-1))*B12D(I-1) + COS(THETA(I-1))*B22D(I-1)
520: Y2 = SIN(THETA(I-1))*B12BULGE + COS(THETA(I-1))*B22BULGE
521: *
522: PHI(I-1) = ATAN2( SQRT(X1**2+X2**2), SQRT(Y1**2+Y2**2) )
523: *
524: * Determine if there are bulges to chase or if a new direct
525: * summand has been reached
526: *
527: RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2
528: RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2
529: RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2
530: RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2
531: *
532: * If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
533: * B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
534: * chasing by applying the original shift again.
535: *
536: IF( .NOT. RESTART11 .AND. .NOT. RESTART21 ) THEN
537: CALL DLARTGP( X2, X1, RWORK(IV1TSN+I-1),
538: $ RWORK(IV1TCS+I-1), R )
539: ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN
540: CALL DLARTGP( B11BULGE, B11E(I-1), RWORK(IV1TSN+I-1),
541: $ RWORK(IV1TCS+I-1), R )
542: ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN
543: CALL DLARTGP( B21BULGE, B21E(I-1), RWORK(IV1TSN+I-1),
544: $ RWORK(IV1TCS+I-1), R )
545: ELSE IF( MU .LE. NU ) THEN
546: CALL DLARTGS( B11D(I), B11E(I), MU, RWORK(IV1TCS+I-1),
547: $ RWORK(IV1TSN+I-1) )
548: ELSE
549: CALL DLARTGS( B21D(I), B21E(I), NU, RWORK(IV1TCS+I-1),
550: $ RWORK(IV1TSN+I-1) )
551: END IF
552: RWORK(IV1TCS+I-1) = -RWORK(IV1TCS+I-1)
553: RWORK(IV1TSN+I-1) = -RWORK(IV1TSN+I-1)
554: IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
555: CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+I-1-1),
556: $ RWORK(IV2TCS+I-1-1), R )
557: ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
558: CALL DLARTGP( B12BULGE, B12D(I-1), RWORK(IV2TSN+I-1-1),
559: $ RWORK(IV2TCS+I-1-1), R )
560: ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
561: CALL DLARTGP( B22BULGE, B22D(I-1), RWORK(IV2TSN+I-1-1),
562: $ RWORK(IV2TCS+I-1-1), R )
563: ELSE IF( NU .LT. MU ) THEN
564: CALL DLARTGS( B12E(I-1), B12D(I), NU,
565: $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) )
566: ELSE
567: CALL DLARTGS( B22E(I-1), B22D(I), MU,
568: $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) )
569: END IF
570: *
571: TEMP = RWORK(IV1TCS+I-1)*B11D(I) + RWORK(IV1TSN+I-1)*B11E(I)
572: B11E(I) = RWORK(IV1TCS+I-1)*B11E(I) -
573: $ RWORK(IV1TSN+I-1)*B11D(I)
574: B11D(I) = TEMP
575: B11BULGE = RWORK(IV1TSN+I-1)*B11D(I+1)
576: B11D(I+1) = RWORK(IV1TCS+I-1)*B11D(I+1)
577: TEMP = RWORK(IV1TCS+I-1)*B21D(I) + RWORK(IV1TSN+I-1)*B21E(I)
578: B21E(I) = RWORK(IV1TCS+I-1)*B21E(I) -
579: $ RWORK(IV1TSN+I-1)*B21D(I)
580: B21D(I) = TEMP
581: B21BULGE = RWORK(IV1TSN+I-1)*B21D(I+1)
582: B21D(I+1) = RWORK(IV1TCS+I-1)*B21D(I+1)
583: TEMP = RWORK(IV2TCS+I-1-1)*B12E(I-1) +
584: $ RWORK(IV2TSN+I-1-1)*B12D(I)
585: B12D(I) = RWORK(IV2TCS+I-1-1)*B12D(I) -
586: $ RWORK(IV2TSN+I-1-1)*B12E(I-1)
587: B12E(I-1) = TEMP
588: B12BULGE = RWORK(IV2TSN+I-1-1)*B12E(I)
589: B12E(I) = RWORK(IV2TCS+I-1-1)*B12E(I)
590: TEMP = RWORK(IV2TCS+I-1-1)*B22E(I-1) +
591: $ RWORK(IV2TSN+I-1-1)*B22D(I)
592: B22D(I) = RWORK(IV2TCS+I-1-1)*B22D(I) -
593: $ RWORK(IV2TSN+I-1-1)*B22E(I-1)
594: B22E(I-1) = TEMP
595: B22BULGE = RWORK(IV2TSN+I-1-1)*B22E(I)
596: B22E(I) = RWORK(IV2TCS+I-1-1)*B22E(I)
597: *
598: * Compute THETA(I)
599: *
600: X1 = COS(PHI(I-1))*B11D(I) + SIN(PHI(I-1))*B12E(I-1)
601: X2 = COS(PHI(I-1))*B11BULGE + SIN(PHI(I-1))*B12BULGE
602: Y1 = COS(PHI(I-1))*B21D(I) + SIN(PHI(I-1))*B22E(I-1)
603: Y2 = COS(PHI(I-1))*B21BULGE + SIN(PHI(I-1))*B22BULGE
604: *
605: THETA(I) = ATAN2( SQRT(Y1**2+Y2**2), SQRT(X1**2+X2**2) )
606: *
607: * Determine if there are bulges to chase or if a new direct
608: * summand has been reached
609: *
610: RESTART11 = B11D(I)**2 + B11BULGE**2 .LE. THRESH**2
611: RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2
612: RESTART21 = B21D(I)**2 + B21BULGE**2 .LE. THRESH**2
613: RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2
614: *
615: * If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
616: * B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
617: * chasing by applying the original shift again.
618: *
619: IF( .NOT. RESTART11 .AND. .NOT. RESTART12 ) THEN
620: CALL DLARTGP( X2, X1, RWORK(IU1SN+I-1), RWORK(IU1CS+I-1),
621: $ R )
622: ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN
623: CALL DLARTGP( B11BULGE, B11D(I), RWORK(IU1SN+I-1),
624: $ RWORK(IU1CS+I-1), R )
625: ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN
626: CALL DLARTGP( B12BULGE, B12E(I-1), RWORK(IU1SN+I-1),
627: $ RWORK(IU1CS+I-1), R )
628: ELSE IF( MU .LE. NU ) THEN
629: CALL DLARTGS( B11E(I), B11D(I+1), MU, RWORK(IU1CS+I-1),
630: $ RWORK(IU1SN+I-1) )
631: ELSE
632: CALL DLARTGS( B12D(I), B12E(I), NU, RWORK(IU1CS+I-1),
633: $ RWORK(IU1SN+I-1) )
634: END IF
635: IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN
636: CALL DLARTGP( Y2, Y1, RWORK(IU2SN+I-1), RWORK(IU2CS+I-1),
637: $ R )
638: ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN
639: CALL DLARTGP( B21BULGE, B21D(I), RWORK(IU2SN+I-1),
640: $ RWORK(IU2CS+I-1), R )
641: ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN
642: CALL DLARTGP( B22BULGE, B22E(I-1), RWORK(IU2SN+I-1),
643: $ RWORK(IU2CS+I-1), R )
644: ELSE IF( NU .LT. MU ) THEN
645: CALL DLARTGS( B21E(I), B21E(I+1), NU, RWORK(IU2CS+I-1),
646: $ RWORK(IU2SN+I-1) )
647: ELSE
648: CALL DLARTGS( B22D(I), B22E(I), MU, RWORK(IU2CS+I-1),
649: $ RWORK(IU2SN+I-1) )
650: END IF
651: RWORK(IU2CS+I-1) = -RWORK(IU2CS+I-1)
652: RWORK(IU2SN+I-1) = -RWORK(IU2SN+I-1)
653: *
654: TEMP = RWORK(IU1CS+I-1)*B11E(I) + RWORK(IU1SN+I-1)*B11D(I+1)
655: B11D(I+1) = RWORK(IU1CS+I-1)*B11D(I+1) -
656: $ RWORK(IU1SN+I-1)*B11E(I)
657: B11E(I) = TEMP
658: IF( I .LT. IMAX - 1 ) THEN
659: B11BULGE = RWORK(IU1SN+I-1)*B11E(I+1)
660: B11E(I+1) = RWORK(IU1CS+I-1)*B11E(I+1)
661: END IF
662: TEMP = RWORK(IU2CS+I-1)*B21E(I) + RWORK(IU2SN+I-1)*B21D(I+1)
663: B21D(I+1) = RWORK(IU2CS+I-1)*B21D(I+1) -
664: $ RWORK(IU2SN+I-1)*B21E(I)
665: B21E(I) = TEMP
666: IF( I .LT. IMAX - 1 ) THEN
667: B21BULGE = RWORK(IU2SN+I-1)*B21E(I+1)
668: B21E(I+1) = RWORK(IU2CS+I-1)*B21E(I+1)
669: END IF
670: TEMP = RWORK(IU1CS+I-1)*B12D(I) + RWORK(IU1SN+I-1)*B12E(I)
671: B12E(I) = RWORK(IU1CS+I-1)*B12E(I) -
672: $ RWORK(IU1SN+I-1)*B12D(I)
673: B12D(I) = TEMP
674: B12BULGE = RWORK(IU1SN+I-1)*B12D(I+1)
675: B12D(I+1) = RWORK(IU1CS+I-1)*B12D(I+1)
676: TEMP = RWORK(IU2CS+I-1)*B22D(I) + RWORK(IU2SN+I-1)*B22E(I)
677: B22E(I) = RWORK(IU2CS+I-1)*B22E(I) -
678: $ RWORK(IU2SN+I-1)*B22D(I)
679: B22D(I) = TEMP
680: B22BULGE = RWORK(IU2SN+I-1)*B22D(I+1)
681: B22D(I+1) = RWORK(IU2CS+I-1)*B22D(I+1)
682: *
683: END DO
684: *
685: * Compute PHI(IMAX-1)
686: *
687: X1 = SIN(THETA(IMAX-1))*B11E(IMAX-1) +
688: $ COS(THETA(IMAX-1))*B21E(IMAX-1)
689: Y1 = SIN(THETA(IMAX-1))*B12D(IMAX-1) +
690: $ COS(THETA(IMAX-1))*B22D(IMAX-1)
691: Y2 = SIN(THETA(IMAX-1))*B12BULGE + COS(THETA(IMAX-1))*B22BULGE
692: *
693: PHI(IMAX-1) = ATAN2( ABS(X1), SQRT(Y1**2+Y2**2) )
694: *
695: * Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
696: *
697: RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2
698: RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2
699: *
700: IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
701: CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+IMAX-1-1),
702: $ RWORK(IV2TCS+IMAX-1-1), R )
703: ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
704: CALL DLARTGP( B12BULGE, B12D(IMAX-1),
705: $ RWORK(IV2TSN+IMAX-1-1),
706: $ RWORK(IV2TCS+IMAX-1-1), R )
707: ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
708: CALL DLARTGP( B22BULGE, B22D(IMAX-1),
709: $ RWORK(IV2TSN+IMAX-1-1),
710: $ RWORK(IV2TCS+IMAX-1-1), R )
711: ELSE IF( NU .LT. MU ) THEN
712: CALL DLARTGS( B12E(IMAX-1), B12D(IMAX), NU,
713: $ RWORK(IV2TCS+IMAX-1-1),
714: $ RWORK(IV2TSN+IMAX-1-1) )
715: ELSE
716: CALL DLARTGS( B22E(IMAX-1), B22D(IMAX), MU,
717: $ RWORK(IV2TCS+IMAX-1-1),
718: $ RWORK(IV2TSN+IMAX-1-1) )
719: END IF
720: *
721: TEMP = RWORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) +
722: $ RWORK(IV2TSN+IMAX-1-1)*B12D(IMAX)
723: B12D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B12D(IMAX) -
724: $ RWORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1)
725: B12E(IMAX-1) = TEMP
726: TEMP = RWORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) +
727: $ RWORK(IV2TSN+IMAX-1-1)*B22D(IMAX)
728: B22D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B22D(IMAX) -
729: $ RWORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1)
730: B22E(IMAX-1) = TEMP
731: *
732: * Update singular vectors
733: *
734: IF( WANTU1 ) THEN
735: IF( COLMAJOR ) THEN
736: CALL ZLASR( 'R', 'V', 'F', P, IMAX-IMIN+1,
737: $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1),
738: $ U1(1,IMIN), LDU1 )
739: ELSE
740: CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, P,
741: $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1),
742: $ U1(IMIN,1), LDU1 )
743: END IF
744: END IF
745: IF( WANTU2 ) THEN
746: IF( COLMAJOR ) THEN
747: CALL ZLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1,
748: $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1),
749: $ U2(1,IMIN), LDU2 )
750: ELSE
751: CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P,
752: $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1),
753: $ U2(IMIN,1), LDU2 )
754: END IF
755: END IF
756: IF( WANTV1T ) THEN
757: IF( COLMAJOR ) THEN
758: CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q,
759: $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1),
760: $ V1T(IMIN,1), LDV1T )
761: ELSE
762: CALL ZLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1,
763: $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1),
764: $ V1T(1,IMIN), LDV1T )
765: END IF
766: END IF
767: IF( WANTV2T ) THEN
768: IF( COLMAJOR ) THEN
769: CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q,
770: $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1),
771: $ V2T(IMIN,1), LDV2T )
772: ELSE
773: CALL ZLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1,
774: $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1),
775: $ V2T(1,IMIN), LDV2T )
776: END IF
777: END IF
778: *
779: * Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
780: *
781: IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN
782: B11D(IMAX) = -B11D(IMAX)
783: B21D(IMAX) = -B21D(IMAX)
784: IF( WANTV1T ) THEN
785: IF( COLMAJOR ) THEN
786: CALL ZSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T )
787: ELSE
788: CALL ZSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 )
789: END IF
790: END IF
791: END IF
792: *
793: * Compute THETA(IMAX)
794: *
795: X1 = COS(PHI(IMAX-1))*B11D(IMAX) +
796: $ SIN(PHI(IMAX-1))*B12E(IMAX-1)
797: Y1 = COS(PHI(IMAX-1))*B21D(IMAX) +
798: $ SIN(PHI(IMAX-1))*B22E(IMAX-1)
799: *
800: THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) )
801: *
802: * Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
803: * and B22(IMAX,IMAX-1)
804: *
805: IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN
806: B12D(IMAX) = -B12D(IMAX)
807: IF( WANTU1 ) THEN
808: IF( COLMAJOR ) THEN
809: CALL ZSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 )
810: ELSE
811: CALL ZSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 )
812: END IF
813: END IF
814: END IF
815: IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN
816: B22D(IMAX) = -B22D(IMAX)
817: IF( WANTU2 ) THEN
818: IF( COLMAJOR ) THEN
819: CALL ZSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 )
820: ELSE
821: CALL ZSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 )
822: END IF
823: END IF
824: END IF
825: *
826: * Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
827: *
828: IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN
829: IF( WANTV2T ) THEN
830: IF( COLMAJOR ) THEN
831: CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T )
832: ELSE
833: CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 )
834: END IF
835: END IF
836: END IF
837: *
838: * Test for negligible sines or cosines
839: *
840: DO I = IMIN, IMAX
841: IF( THETA(I) .LT. THRESH ) THEN
842: THETA(I) = ZERO
843: ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
844: THETA(I) = PIOVER2
845: END IF
846: END DO
847: DO I = IMIN, IMAX-1
848: IF( PHI(I) .LT. THRESH ) THEN
849: PHI(I) = ZERO
850: ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
851: PHI(I) = PIOVER2
852: END IF
853: END DO
854: *
855: * Deflate
856: *
857: IF (IMAX .GT. 1) THEN
858: DO WHILE( PHI(IMAX-1) .EQ. ZERO )
859: IMAX = IMAX - 1
860: IF (IMAX .LE. 1) EXIT
861: END DO
862: END IF
863: IF( IMIN .GT. IMAX - 1 )
864: $ IMIN = IMAX - 1
865: IF (IMIN .GT. 1) THEN
866: DO WHILE (PHI(IMIN-1) .NE. ZERO)
867: IMIN = IMIN - 1
868: IF (IMIN .LE. 1) EXIT
869: END DO
870: END IF
871: *
872: * Repeat main iteration loop
873: *
874: END DO
875: *
876: * Postprocessing: order THETA from least to greatest
877: *
878: DO I = 1, Q
879: *
880: MINI = I
881: THETAMIN = THETA(I)
882: DO J = I+1, Q
883: IF( THETA(J) .LT. THETAMIN ) THEN
884: MINI = J
885: THETAMIN = THETA(J)
886: END IF
887: END DO
888: *
889: IF( MINI .NE. I ) THEN
890: THETA(MINI) = THETA(I)
891: THETA(I) = THETAMIN
892: IF( COLMAJOR ) THEN
893: IF( WANTU1 )
894: $ CALL ZSWAP( P, U1(1,I), 1, U1(1,MINI), 1 )
895: IF( WANTU2 )
896: $ CALL ZSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 )
897: IF( WANTV1T )
898: $ CALL ZSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T )
899: IF( WANTV2T )
900: $ CALL ZSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1),
901: $ LDV2T )
902: ELSE
903: IF( WANTU1 )
904: $ CALL ZSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 )
905: IF( WANTU2 )
906: $ CALL ZSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 )
907: IF( WANTV1T )
908: $ CALL ZSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 )
909: IF( WANTV2T )
910: $ CALL ZSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 )
911: END IF
912: END IF
913: *
914: END DO
915: *
916: RETURN
917: *
918: * End of ZBBCSD
919: *
920: END
921:
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zbbcsd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack