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