LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
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, a P-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, LDU1 >= MAX(1,P).
161 *> \endverbatim
162 *>
163 *> \param[in,out] U2
164 *> \verbatim
165 *> U2 is COMPLEX*16 array, dimension (LDU2,M-P)
166 *> On entry, an (M-P)-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, LDU2 >= MAX(1,M-P).
175 *> \endverbatim
176 *>
177 *> \param[in,out] V1T
178 *> \verbatim
179 *> V1T is COMPLEX*16 array, dimension (LDV1T,Q)
180 *> On entry, a Q-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, LDV1T >= MAX(1,Q).
189 *> \endverbatim
190 *>
191 *> \param[in,out] V2T
192 *> \verbatim
193 *> V2T is COMPLEX*16 array, dimension (LDV2T,M-Q)
194 *> On entry, an (M-Q)-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, LDV2T >= MAX(1,M-Q).
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 ZBBCSD converges, B21D contains the negative sines of
244 *> THETA(1), ..., THETA(Q). If ZBBCSD 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 ZBBCSD converges, B21E contains zeros. If ZBBCSD 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 ZBBCSD converges, B22D contains the negative sines of
261 *> THETA(1), ..., THETA(Q). If ZBBCSD 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 ZBBCSD converges, B22E contains zeros. If ZBBCSD 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,LRWORK))
277 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
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 June 2016
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.7.1) --
336 * -- LAPACK is a software package provided by Univ. of Tennessee, --
337 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
338 * June 2016
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 
subroutine zlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
ZLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: zlasr.f:202
subroutine dlartgs(X, Y, SIGMA, CS, SN)
DLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bid...
Definition: dlartgs.f:92
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:83
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlartgp(F, G, CS, SN, R)
DLARTGP generates a plane rotation so that the diagonal is nonnegative.
Definition: dlartgp.f:97
subroutine zbbcsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E, B22D, B22E, RWORK, LRWORK, INFO)
ZBBCSD
Definition: zbbcsd.f:334
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: dlas2.f:109