LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 2011
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.4.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 2011
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