LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
zbbcsd.f
Go to the documentation of this file.
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 2013
326 *
327 *> \ingroup complex16OTHERcomputational
328 *
329 * =====================================================================
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 *
335 * -- LAPACK computational routine (version 3.5.0) --
336 * -- LAPACK is a software package provided by Univ. of Tennessee, --
337 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
338 * November 2013
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 )
480  IF( phi(imax-1) .NE. zero ) THEN
481  EXIT
482  END IF
483  imax = imax - 1
484  END DO
485  imin = imax - 1
486  IF ( imin .GT. 1 ) THEN
487  DO WHILE( phi(imin-1) .NE. zero )
488  imin = imin - 1
489  IF ( imin .LE. 1 ) EXIT
490  END DO
491  END IF
492 *
493 * Initialize iteration counter
494 *
495  maxit = maxitr*q*q
496  iter = 0
497 *
498 * Begin main iteration loop
499 *
500  DO WHILE( imax .GT. 1 )
501 *
502 * Compute the matrix entries
503 *
504  b11d(imin) = cos( theta(imin) )
505  b21d(imin) = -sin( theta(imin) )
506  DO i = imin, imax - 1
507  b11e(i) = -sin( theta(i) ) * sin( phi(i) )
508  b11d(i+1) = cos( theta(i+1) ) * cos( phi(i) )
509  b12d(i) = sin( theta(i) ) * cos( phi(i) )
510  b12e(i) = cos( theta(i+1) ) * sin( phi(i) )
511  b21e(i) = -cos( theta(i) ) * sin( phi(i) )
512  b21d(i+1) = -sin( theta(i+1) ) * cos( phi(i) )
513  b22d(i) = cos( theta(i) ) * cos( phi(i) )
514  b22e(i) = -sin( theta(i+1) ) * sin( phi(i) )
515  END DO
516  b12d(imax) = sin( theta(imax) )
517  b22d(imax) = cos( theta(imax) )
518 *
519 * Abort if not converging; otherwise, increment ITER
520 *
521  IF( iter .GT. maxit ) THEN
522  info = 0
523  DO i = 1, q
524  IF( phi(i) .NE. zero )
525  $ info = info + 1
526  END DO
527  RETURN
528  END IF
529 *
530  iter = iter + imax - imin
531 *
532 * Compute shifts
533 *
534  thetamax = theta(imin)
535  thetamin = theta(imin)
536  DO i = imin+1, imax
537  IF( theta(i) > thetamax )
538  $ thetamax = theta(i)
539  IF( theta(i) < thetamin )
540  $ thetamin = theta(i)
541  END DO
542 *
543  IF( thetamax .GT. piover2 - thresh ) THEN
544 *
545 * Zero on diagonals of B11 and B22; induce deflation with a
546 * zero shift
547 *
548  mu = zero
549  nu = one
550 *
551  ELSE IF( thetamin .LT. thresh ) THEN
552 *
553 * Zero on diagonals of B12 and B22; induce deflation with a
554 * zero shift
555 *
556  mu = one
557  nu = zero
558 *
559  ELSE
560 *
561 * Compute shifts for B11 and B21 and use the lesser
562 *
563  CALL dlas2( b11d(imax-1), b11e(imax-1), b11d(imax), sigma11,
564  $ dummy )
565  CALL dlas2( b21d(imax-1), b21e(imax-1), b21d(imax), sigma21,
566  $ dummy )
567 *
568  IF( sigma11 .LE. sigma21 ) THEN
569  mu = sigma11
570  nu = sqrt( one - mu**2 )
571  IF( mu .LT. thresh ) THEN
572  mu = zero
573  nu = one
574  END IF
575  ELSE
576  nu = sigma21
577  mu = sqrt( 1.0 - nu**2 )
578  IF( nu .LT. thresh ) THEN
579  mu = one
580  nu = zero
581  END IF
582  END IF
583  END IF
584 *
585 * Rotate to produce bulges in B11 and B21
586 *
587  IF( mu .LE. nu ) THEN
588  CALL dlartgs( b11d(imin), b11e(imin), mu,
589  $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
590  ELSE
591  CALL dlartgs( b21d(imin), b21e(imin), nu,
592  $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
593  END IF
594 *
595  temp = rwork(iv1tcs+imin-1)*b11d(imin) +
596  $ rwork(iv1tsn+imin-1)*b11e(imin)
597  b11e(imin) = rwork(iv1tcs+imin-1)*b11e(imin) -
598  $ rwork(iv1tsn+imin-1)*b11d(imin)
599  b11d(imin) = temp
600  b11bulge = rwork(iv1tsn+imin-1)*b11d(imin+1)
601  b11d(imin+1) = rwork(iv1tcs+imin-1)*b11d(imin+1)
602  temp = rwork(iv1tcs+imin-1)*b21d(imin) +
603  $ rwork(iv1tsn+imin-1)*b21e(imin)
604  b21e(imin) = rwork(iv1tcs+imin-1)*b21e(imin) -
605  $ rwork(iv1tsn+imin-1)*b21d(imin)
606  b21d(imin) = temp
607  b21bulge = rwork(iv1tsn+imin-1)*b21d(imin+1)
608  b21d(imin+1) = rwork(iv1tcs+imin-1)*b21d(imin+1)
609 *
610 * Compute THETA(IMIN)
611 *
612  theta( imin ) = atan2( sqrt( b21d(imin)**2+b21bulge**2 ),
613  $ sqrt( b11d(imin)**2+b11bulge**2 ) )
614 *
615 * Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
616 *
617  IF( b11d(imin)**2+b11bulge**2 .GT. thresh**2 ) THEN
618  CALL dlartgp( b11bulge, b11d(imin), rwork(iu1sn+imin-1),
619  $ rwork(iu1cs+imin-1), r )
620  ELSE IF( mu .LE. nu ) THEN
621  CALL dlartgs( b11e( imin ), b11d( imin + 1 ), mu,
622  $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
623  ELSE
624  CALL dlartgs( b12d( imin ), b12e( imin ), nu,
625  $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
626  END IF
627  IF( b21d(imin)**2+b21bulge**2 .GT. thresh**2 ) THEN
628  CALL dlartgp( b21bulge, b21d(imin), rwork(iu2sn+imin-1),
629  $ rwork(iu2cs+imin-1), r )
630  ELSE IF( nu .LT. mu ) THEN
631  CALL dlartgs( b21e( imin ), b21d( imin + 1 ), nu,
632  $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
633  ELSE
634  CALL dlartgs( b22d(imin), b22e(imin), mu,
635  $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
636  END IF
637  rwork(iu2cs+imin-1) = -rwork(iu2cs+imin-1)
638  rwork(iu2sn+imin-1) = -rwork(iu2sn+imin-1)
639 *
640  temp = rwork(iu1cs+imin-1)*b11e(imin) +
641  $ rwork(iu1sn+imin-1)*b11d(imin+1)
642  b11d(imin+1) = rwork(iu1cs+imin-1)*b11d(imin+1) -
643  $ rwork(iu1sn+imin-1)*b11e(imin)
644  b11e(imin) = temp
645  IF( imax .GT. imin+1 ) THEN
646  b11bulge = rwork(iu1sn+imin-1)*b11e(imin+1)
647  b11e(imin+1) = rwork(iu1cs+imin-1)*b11e(imin+1)
648  END IF
649  temp = rwork(iu1cs+imin-1)*b12d(imin) +
650  $ rwork(iu1sn+imin-1)*b12e(imin)
651  b12e(imin) = rwork(iu1cs+imin-1)*b12e(imin) -
652  $ rwork(iu1sn+imin-1)*b12d(imin)
653  b12d(imin) = temp
654  b12bulge = rwork(iu1sn+imin-1)*b12d(imin+1)
655  b12d(imin+1) = rwork(iu1cs+imin-1)*b12d(imin+1)
656  temp = rwork(iu2cs+imin-1)*b21e(imin) +
657  $ rwork(iu2sn+imin-1)*b21d(imin+1)
658  b21d(imin+1) = rwork(iu2cs+imin-1)*b21d(imin+1) -
659  $ rwork(iu2sn+imin-1)*b21e(imin)
660  b21e(imin) = temp
661  IF( imax .GT. imin+1 ) THEN
662  b21bulge = rwork(iu2sn+imin-1)*b21e(imin+1)
663  b21e(imin+1) = rwork(iu2cs+imin-1)*b21e(imin+1)
664  END IF
665  temp = rwork(iu2cs+imin-1)*b22d(imin) +
666  $ rwork(iu2sn+imin-1)*b22e(imin)
667  b22e(imin) = rwork(iu2cs+imin-1)*b22e(imin) -
668  $ rwork(iu2sn+imin-1)*b22d(imin)
669  b22d(imin) = temp
670  b22bulge = rwork(iu2sn+imin-1)*b22d(imin+1)
671  b22d(imin+1) = rwork(iu2cs+imin-1)*b22d(imin+1)
672 *
673 * Inner loop: chase bulges from B11(IMIN,IMIN+2),
674 * B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
675 * bottom-right
676 *
677  DO i = imin+1, imax-1
678 *
679 * Compute PHI(I-1)
680 *
681  x1 = sin(theta(i-1))*b11e(i-1) + cos(theta(i-1))*b21e(i-1)
682  x2 = sin(theta(i-1))*b11bulge + cos(theta(i-1))*b21bulge
683  y1 = sin(theta(i-1))*b12d(i-1) + cos(theta(i-1))*b22d(i-1)
684  y2 = sin(theta(i-1))*b12bulge + cos(theta(i-1))*b22bulge
685 *
686  phi(i-1) = atan2( sqrt(x1**2+x2**2), sqrt(y1**2+y2**2) )
687 *
688 * Determine if there are bulges to chase or if a new direct
689 * summand has been reached
690 *
691  restart11 = b11e(i-1)**2 + b11bulge**2 .LE. thresh**2
692  restart21 = b21e(i-1)**2 + b21bulge**2 .LE. thresh**2
693  restart12 = b12d(i-1)**2 + b12bulge**2 .LE. thresh**2
694  restart22 = b22d(i-1)**2 + b22bulge**2 .LE. thresh**2
695 *
696 * If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
697 * B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
698 * chasing by applying the original shift again.
699 *
700  IF( .NOT. restart11 .AND. .NOT. restart21 ) THEN
701  CALL dlartgp( x2, x1, rwork(iv1tsn+i-1),
702  $ rwork(iv1tcs+i-1), r )
703  ELSE IF( .NOT. restart11 .AND. restart21 ) THEN
704  CALL dlartgp( b11bulge, b11e(i-1), rwork(iv1tsn+i-1),
705  $ rwork(iv1tcs+i-1), r )
706  ELSE IF( restart11 .AND. .NOT. restart21 ) THEN
707  CALL dlartgp( b21bulge, b21e(i-1), rwork(iv1tsn+i-1),
708  $ rwork(iv1tcs+i-1), r )
709  ELSE IF( mu .LE. nu ) THEN
710  CALL dlartgs( b11d(i), b11e(i), mu, rwork(iv1tcs+i-1),
711  $ rwork(iv1tsn+i-1) )
712  ELSE
713  CALL dlartgs( b21d(i), b21e(i), nu, rwork(iv1tcs+i-1),
714  $ rwork(iv1tsn+i-1) )
715  END IF
716  rwork(iv1tcs+i-1) = -rwork(iv1tcs+i-1)
717  rwork(iv1tsn+i-1) = -rwork(iv1tsn+i-1)
718  IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
719  CALL dlartgp( y2, y1, rwork(iv2tsn+i-1-1),
720  $ rwork(iv2tcs+i-1-1), r )
721  ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
722  CALL dlartgp( b12bulge, b12d(i-1), rwork(iv2tsn+i-1-1),
723  $ rwork(iv2tcs+i-1-1), r )
724  ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
725  CALL dlartgp( b22bulge, b22d(i-1), rwork(iv2tsn+i-1-1),
726  $ rwork(iv2tcs+i-1-1), r )
727  ELSE IF( nu .LT. mu ) THEN
728  CALL dlartgs( b12e(i-1), b12d(i), nu,
729  $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
730  ELSE
731  CALL dlartgs( b22e(i-1), b22d(i), mu,
732  $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
733  END IF
734 *
735  temp = rwork(iv1tcs+i-1)*b11d(i) + rwork(iv1tsn+i-1)*b11e(i)
736  b11e(i) = rwork(iv1tcs+i-1)*b11e(i) -
737  $ rwork(iv1tsn+i-1)*b11d(i)
738  b11d(i) = temp
739  b11bulge = rwork(iv1tsn+i-1)*b11d(i+1)
740  b11d(i+1) = rwork(iv1tcs+i-1)*b11d(i+1)
741  temp = rwork(iv1tcs+i-1)*b21d(i) + rwork(iv1tsn+i-1)*b21e(i)
742  b21e(i) = rwork(iv1tcs+i-1)*b21e(i) -
743  $ rwork(iv1tsn+i-1)*b21d(i)
744  b21d(i) = temp
745  b21bulge = rwork(iv1tsn+i-1)*b21d(i+1)
746  b21d(i+1) = rwork(iv1tcs+i-1)*b21d(i+1)
747  temp = rwork(iv2tcs+i-1-1)*b12e(i-1) +
748  $ rwork(iv2tsn+i-1-1)*b12d(i)
749  b12d(i) = rwork(iv2tcs+i-1-1)*b12d(i) -
750  $ rwork(iv2tsn+i-1-1)*b12e(i-1)
751  b12e(i-1) = temp
752  b12bulge = rwork(iv2tsn+i-1-1)*b12e(i)
753  b12e(i) = rwork(iv2tcs+i-1-1)*b12e(i)
754  temp = rwork(iv2tcs+i-1-1)*b22e(i-1) +
755  $ rwork(iv2tsn+i-1-1)*b22d(i)
756  b22d(i) = rwork(iv2tcs+i-1-1)*b22d(i) -
757  $ rwork(iv2tsn+i-1-1)*b22e(i-1)
758  b22e(i-1) = temp
759  b22bulge = rwork(iv2tsn+i-1-1)*b22e(i)
760  b22e(i) = rwork(iv2tcs+i-1-1)*b22e(i)
761 *
762 * Compute THETA(I)
763 *
764  x1 = cos(phi(i-1))*b11d(i) + sin(phi(i-1))*b12e(i-1)
765  x2 = cos(phi(i-1))*b11bulge + sin(phi(i-1))*b12bulge
766  y1 = cos(phi(i-1))*b21d(i) + sin(phi(i-1))*b22e(i-1)
767  y2 = cos(phi(i-1))*b21bulge + sin(phi(i-1))*b22bulge
768 *
769  theta(i) = atan2( sqrt(y1**2+y2**2), sqrt(x1**2+x2**2) )
770 *
771 * Determine if there are bulges to chase or if a new direct
772 * summand has been reached
773 *
774  restart11 = b11d(i)**2 + b11bulge**2 .LE. thresh**2
775  restart12 = b12e(i-1)**2 + b12bulge**2 .LE. thresh**2
776  restart21 = b21d(i)**2 + b21bulge**2 .LE. thresh**2
777  restart22 = b22e(i-1)**2 + b22bulge**2 .LE. thresh**2
778 *
779 * If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
780 * B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
781 * chasing by applying the original shift again.
782 *
783  IF( .NOT. restart11 .AND. .NOT. restart12 ) THEN
784  CALL dlartgp( x2, x1, rwork(iu1sn+i-1), rwork(iu1cs+i-1),
785  $ r )
786  ELSE IF( .NOT. restart11 .AND. restart12 ) THEN
787  CALL dlartgp( b11bulge, b11d(i), rwork(iu1sn+i-1),
788  $ rwork(iu1cs+i-1), r )
789  ELSE IF( restart11 .AND. .NOT. restart12 ) THEN
790  CALL dlartgp( b12bulge, b12e(i-1), rwork(iu1sn+i-1),
791  $ rwork(iu1cs+i-1), r )
792  ELSE IF( mu .LE. nu ) THEN
793  CALL dlartgs( b11e(i), b11d(i+1), mu, rwork(iu1cs+i-1),
794  $ rwork(iu1sn+i-1) )
795  ELSE
796  CALL dlartgs( b12d(i), b12e(i), nu, rwork(iu1cs+i-1),
797  $ rwork(iu1sn+i-1) )
798  END IF
799  IF( .NOT. restart21 .AND. .NOT. restart22 ) THEN
800  CALL dlartgp( y2, y1, rwork(iu2sn+i-1), rwork(iu2cs+i-1),
801  $ r )
802  ELSE IF( .NOT. restart21 .AND. restart22 ) THEN
803  CALL dlartgp( b21bulge, b21d(i), rwork(iu2sn+i-1),
804  $ rwork(iu2cs+i-1), r )
805  ELSE IF( restart21 .AND. .NOT. restart22 ) THEN
806  CALL dlartgp( b22bulge, b22e(i-1), rwork(iu2sn+i-1),
807  $ rwork(iu2cs+i-1), r )
808  ELSE IF( nu .LT. mu ) THEN
809  CALL dlartgs( b21e(i), b21e(i+1), nu, rwork(iu2cs+i-1),
810  $ rwork(iu2sn+i-1) )
811  ELSE
812  CALL dlartgs( b22d(i), b22e(i), mu, rwork(iu2cs+i-1),
813  $ rwork(iu2sn+i-1) )
814  END IF
815  rwork(iu2cs+i-1) = -rwork(iu2cs+i-1)
816  rwork(iu2sn+i-1) = -rwork(iu2sn+i-1)
817 *
818  temp = rwork(iu1cs+i-1)*b11e(i) + rwork(iu1sn+i-1)*b11d(i+1)
819  b11d(i+1) = rwork(iu1cs+i-1)*b11d(i+1) -
820  $ rwork(iu1sn+i-1)*b11e(i)
821  b11e(i) = temp
822  IF( i .LT. imax - 1 ) THEN
823  b11bulge = rwork(iu1sn+i-1)*b11e(i+1)
824  b11e(i+1) = rwork(iu1cs+i-1)*b11e(i+1)
825  END IF
826  temp = rwork(iu2cs+i-1)*b21e(i) + rwork(iu2sn+i-1)*b21d(i+1)
827  b21d(i+1) = rwork(iu2cs+i-1)*b21d(i+1) -
828  $ rwork(iu2sn+i-1)*b21e(i)
829  b21e(i) = temp
830  IF( i .LT. imax - 1 ) THEN
831  b21bulge = rwork(iu2sn+i-1)*b21e(i+1)
832  b21e(i+1) = rwork(iu2cs+i-1)*b21e(i+1)
833  END IF
834  temp = rwork(iu1cs+i-1)*b12d(i) + rwork(iu1sn+i-1)*b12e(i)
835  b12e(i) = rwork(iu1cs+i-1)*b12e(i) -
836  $ rwork(iu1sn+i-1)*b12d(i)
837  b12d(i) = temp
838  b12bulge = rwork(iu1sn+i-1)*b12d(i+1)
839  b12d(i+1) = rwork(iu1cs+i-1)*b12d(i+1)
840  temp = rwork(iu2cs+i-1)*b22d(i) + rwork(iu2sn+i-1)*b22e(i)
841  b22e(i) = rwork(iu2cs+i-1)*b22e(i) -
842  $ rwork(iu2sn+i-1)*b22d(i)
843  b22d(i) = temp
844  b22bulge = rwork(iu2sn+i-1)*b22d(i+1)
845  b22d(i+1) = rwork(iu2cs+i-1)*b22d(i+1)
846 *
847  END DO
848 *
849 * Compute PHI(IMAX-1)
850 *
851  x1 = sin(theta(imax-1))*b11e(imax-1) +
852  $ cos(theta(imax-1))*b21e(imax-1)
853  y1 = sin(theta(imax-1))*b12d(imax-1) +
854  $ cos(theta(imax-1))*b22d(imax-1)
855  y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
856 *
857  phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
858 *
859 * Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
860 *
861  restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
862  restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
863 *
864  IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
865  CALL dlartgp( y2, y1, rwork(iv2tsn+imax-1-1),
866  $ rwork(iv2tcs+imax-1-1), r )
867  ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
868  CALL dlartgp( b12bulge, b12d(imax-1),
869  $ rwork(iv2tsn+imax-1-1),
870  $ rwork(iv2tcs+imax-1-1), r )
871  ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
872  CALL dlartgp( b22bulge, b22d(imax-1),
873  $ rwork(iv2tsn+imax-1-1),
874  $ rwork(iv2tcs+imax-1-1), r )
875  ELSE IF( nu .LT. mu ) THEN
876  CALL dlartgs( b12e(imax-1), b12d(imax), nu,
877  $ rwork(iv2tcs+imax-1-1),
878  $ rwork(iv2tsn+imax-1-1) )
879  ELSE
880  CALL dlartgs( b22e(imax-1), b22d(imax), mu,
881  $ rwork(iv2tcs+imax-1-1),
882  $ rwork(iv2tsn+imax-1-1) )
883  END IF
884 *
885  temp = rwork(iv2tcs+imax-1-1)*b12e(imax-1) +
886  $ rwork(iv2tsn+imax-1-1)*b12d(imax)
887  b12d(imax) = rwork(iv2tcs+imax-1-1)*b12d(imax) -
888  $ rwork(iv2tsn+imax-1-1)*b12e(imax-1)
889  b12e(imax-1) = temp
890  temp = rwork(iv2tcs+imax-1-1)*b22e(imax-1) +
891  $ rwork(iv2tsn+imax-1-1)*b22d(imax)
892  b22d(imax) = rwork(iv2tcs+imax-1-1)*b22d(imax) -
893  $ rwork(iv2tsn+imax-1-1)*b22e(imax-1)
894  b22e(imax-1) = temp
895 *
896 * Update singular vectors
897 *
898  IF( wantu1 ) THEN
899  IF( colmajor ) THEN
900  CALL zlasr( 'R', 'V', 'F', p, imax-imin+1,
901  $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
902  $ u1(1,imin), ldu1 )
903  ELSE
904  CALL zlasr( 'L', 'V', 'F', imax-imin+1, p,
905  $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
906  $ u1(imin,1), ldu1 )
907  END IF
908  END IF
909  IF( wantu2 ) THEN
910  IF( colmajor ) THEN
911  CALL zlasr( 'R', 'V', 'F', m-p, imax-imin+1,
912  $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
913  $ u2(1,imin), ldu2 )
914  ELSE
915  CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-p,
916  $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
917  $ u2(imin,1), ldu2 )
918  END IF
919  END IF
920  IF( wantv1t ) THEN
921  IF( colmajor ) THEN
922  CALL zlasr( 'L', 'V', 'F', imax-imin+1, q,
923  $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
924  $ v1t(imin,1), ldv1t )
925  ELSE
926  CALL zlasr( 'R', 'V', 'F', q, imax-imin+1,
927  $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
928  $ v1t(1,imin), ldv1t )
929  END IF
930  END IF
931  IF( wantv2t ) THEN
932  IF( colmajor ) THEN
933  CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-q,
934  $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
935  $ v2t(imin,1), ldv2t )
936  ELSE
937  CALL zlasr( 'R', 'V', 'F', m-q, imax-imin+1,
938  $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
939  $ v2t(1,imin), ldv2t )
940  END IF
941  END IF
942 *
943 * Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
944 *
945  IF( b11e(imax-1)+b21e(imax-1) .GT. 0 ) THEN
946  b11d(imax) = -b11d(imax)
947  b21d(imax) = -b21d(imax)
948  IF( wantv1t ) THEN
949  IF( colmajor ) THEN
950  CALL zscal( q, negonecomplex, v1t(imax,1), ldv1t )
951  ELSE
952  CALL zscal( q, negonecomplex, v1t(1,imax), 1 )
953  END IF
954  END IF
955  END IF
956 *
957 * Compute THETA(IMAX)
958 *
959  x1 = cos(phi(imax-1))*b11d(imax) +
960  $ sin(phi(imax-1))*b12e(imax-1)
961  y1 = cos(phi(imax-1))*b21d(imax) +
962  $ sin(phi(imax-1))*b22e(imax-1)
963 *
964  theta(imax) = atan2( abs(y1), abs(x1) )
965 *
966 * Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
967 * and B22(IMAX,IMAX-1)
968 *
969  IF( b11d(imax)+b12e(imax-1) .LT. 0 ) THEN
970  b12d(imax) = -b12d(imax)
971  IF( wantu1 ) THEN
972  IF( colmajor ) THEN
973  CALL zscal( p, negonecomplex, u1(1,imax), 1 )
974  ELSE
975  CALL zscal( p, negonecomplex, u1(imax,1), ldu1 )
976  END IF
977  END IF
978  END IF
979  IF( b21d(imax)+b22e(imax-1) .GT. 0 ) THEN
980  b22d(imax) = -b22d(imax)
981  IF( wantu2 ) THEN
982  IF( colmajor ) THEN
983  CALL zscal( m-p, negonecomplex, u2(1,imax), 1 )
984  ELSE
985  CALL zscal( m-p, negonecomplex, u2(imax,1), ldu2 )
986  END IF
987  END IF
988  END IF
989 *
990 * Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
991 *
992  IF( b12d(imax)+b22d(imax) .LT. 0 ) THEN
993  IF( wantv2t ) THEN
994  IF( colmajor ) THEN
995  CALL zscal( m-q, negonecomplex, v2t(imax,1), ldv2t )
996  ELSE
997  CALL zscal( m-q, negonecomplex, v2t(1,imax), 1 )
998  END IF
999  END IF
1000  END IF
1001 *
1002 * Test for negligible sines or cosines
1003 *
1004  DO i = imin, imax
1005  IF( theta(i) .LT. thresh ) THEN
1006  theta(i) = zero
1007  ELSE IF( theta(i) .GT. piover2-thresh ) THEN
1008  theta(i) = piover2
1009  END IF
1010  END DO
1011  DO i = imin, imax-1
1012  IF( phi(i) .LT. thresh ) THEN
1013  phi(i) = zero
1014  ELSE IF( phi(i) .GT. piover2-thresh ) THEN
1015  phi(i) = piover2
1016  END IF
1017  END DO
1018 *
1019 * Deflate
1020 *
1021  IF (imax .GT. 1) THEN
1022  DO WHILE( phi(imax-1) .EQ. zero )
1023  imax = imax - 1
1024  IF (imax .LE. 1) EXIT
1025  END DO
1026  END IF
1027  IF( imin .GT. imax - 1 )
1028  $ imin = imax - 1
1029  IF (imin .GT. 1) THEN
1030  DO WHILE (phi(imin-1) .NE. zero)
1031  imin = imin - 1
1032  IF (imin .LE. 1) EXIT
1033  END DO
1034  END IF
1035 *
1036 * Repeat main iteration loop
1037 *
1038  END DO
1039 *
1040 * Postprocessing: order THETA from least to greatest
1041 *
1042  DO i = 1, q
1043 *
1044  mini = i
1045  thetamin = theta(i)
1046  DO j = i+1, q
1047  IF( theta(j) .LT. thetamin ) THEN
1048  mini = j
1049  thetamin = theta(j)
1050  END IF
1051  END DO
1052 *
1053  IF( mini .NE. i ) THEN
1054  theta(mini) = theta(i)
1055  theta(i) = thetamin
1056  IF( colmajor ) THEN
1057  IF( wantu1 )
1058  $ CALL zswap( p, u1(1,i), 1, u1(1,mini), 1 )
1059  IF( wantu2 )
1060  $ CALL zswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1061  IF( wantv1t )
1062  $ CALL zswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1063  IF( wantv2t )
1064  $ CALL zswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1065  $ ldv2t )
1066  ELSE
1067  IF( wantu1 )
1068  $ CALL zswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1069  IF( wantu2 )
1070  $ CALL zswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1071  IF( wantv1t )
1072  $ CALL zswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1073  IF( wantv2t )
1074  $ CALL zswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
1075  END IF
1076  END IF
1077 *
1078  END DO
1079 *
1080  RETURN
1081 *
1082 * End of ZBBCSD
1083 *
1084  END
1085