LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dbbcsd.f
Go to the documentation of this file.
1 *> \brief \b DBBCSD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DBBCSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbbcsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbbcsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbbcsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DBBCSD( 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, WORK, LWORK, INFO )
25 *
26 * .. Scalar Arguments ..
27 * CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
28 * INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
29 * ..
30 * .. Array Arguments ..
31 * DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
32 * $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
33 * $ PHI( * ), THETA( * ), WORK( * )
34 * DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
35 * $ V2T( LDV2T, * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> DBBCSD computes the CS decomposition of an orthogonal 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 | ]**T
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 DORCSD 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 orthogonal 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 orthogonal 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDV1T,Q)
180 *> On entry, a Q-by-Q matrix. On exit, V1T is premultiplied
181 *> by the 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 DOUBLE PRECISION array, dimension (LDV2T,M-Q)
194 *> On entry, an (M-Q)-by-(M-Q) matrix. On exit, V2T is
195 *> premultiplied by the 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 DBBCSD converges, B11D contains the cosines of THETA(1),
210 *> ..., THETA(Q). If DBBCSD 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 DBBCSD converges, B11E contains zeros. If DBBCSD 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 DBBCSD converges, B12D contains the negative sines of
227 *> THETA(1), ..., THETA(Q). If DBBCSD 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 DBBCSD converges, B12E contains zeros. If DBBCSD 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 DBBCSD converges, B21D contains the negative sines of
244 *> THETA(1), ..., THETA(Q). If DBBCSD 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 DBBCSD converges, B21E contains zeros. If DBBCSD 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 DBBCSD converges, B22D contains the negative sines of
261 *> THETA(1), ..., THETA(Q). If DBBCSD 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 DBBCSD converges, B22E contains zeros. If DBBCSD fails
270 *> to converge, then B22E contains the subdiagonal of the
271 *> partially reduced bottom-right block.
272 *> \endverbatim
273 *>
274 *> \param[out] WORK
275 *> \verbatim
276 *> WORK 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] LWORK
281 *> \verbatim
282 *> LWORK is INTEGER
283 *> The dimension of the array WORK. LWORK >= MAX(1,8*Q).
284 *>
285 *> If LWORK = -1, then a workspace query is assumed; the
286 *> routine only calculates the optimal size of the WORK array,
287 *> returns this value as the first entry of the work array, and
288 *> no error message related to LWORK 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 DBBCSD 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 *> \ingroup doubleOTHERcomputational
326 *
327 * =====================================================================
328  SUBROUTINE dbbcsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
329  $ THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T,
330  $ V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E,
331  $ B22D, B22E, WORK, LWORK, INFO )
332 *
333 * -- LAPACK computational routine --
334 * -- LAPACK is a software package provided by Univ. of Tennessee, --
335 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
336 *
337 * .. Scalar Arguments ..
338  CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
339  INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
340 * ..
341 * .. Array Arguments ..
342  DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
343  $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
344  $ PHI( * ), THETA( * ), WORK( * )
345  DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
346  $ V2T( LDV2T, * )
347 * ..
348 *
349 * ===================================================================
350 *
351 * .. Parameters ..
352  INTEGER MAXITR
353  PARAMETER ( MAXITR = 6 )
354  double precision hundred, meighth, one, ten, zero
355  parameter( hundred = 100.0d0, meighth = -0.125d0,
356  $ one = 1.0d0, ten = 10.0d0, zero = 0.0d0 )
357  DOUBLE PRECISION NEGONE
358  parameter( negone = -1.0d0 )
359  DOUBLE PRECISION PIOVER2
360  parameter( piover2 = 1.57079632679489661923132169163975144210d0 )
361 * ..
362 * .. Local Scalars ..
363  LOGICAL COLMAJOR, LQUERY, RESTART11, RESTART12,
364  $ RESTART21, RESTART22, WANTU1, WANTU2, WANTV1T,
365  $ WANTV2T
366  INTEGER I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS,
367  $ IU2SN, IV1TCS, IV1TSN, IV2TCS, IV2TSN, J,
368  $ LWORKMIN, LWORKOPT, MAXIT, MINI
369  DOUBLE PRECISION B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY,
370  $ EPS, MU, NU, R, SIGMA11, SIGMA21,
371  $ TEMP, THETAMAX, THETAMIN, THRESH, TOL, TOLMUL,
372  $ unfl, x1, x2, y1, y2
373 *
374 * .. External Subroutines ..
375  EXTERNAL dlasr, dscal, dswap, dlartgp, dlartgs, dlas2,
376  $ xerbla
377 * ..
378 * .. External Functions ..
379  DOUBLE PRECISION DLAMCH
380  LOGICAL LSAME
381  EXTERNAL LSAME, DLAMCH
382 * ..
383 * .. Intrinsic Functions ..
384  INTRINSIC abs, atan2, cos, max, min, sin, sqrt
385 * ..
386 * .. Executable Statements ..
387 *
388 * Test input arguments
389 *
390  info = 0
391  lquery = lwork .EQ. -1
392  wantu1 = lsame( jobu1, 'Y' )
393  wantu2 = lsame( jobu2, 'Y' )
394  wantv1t = lsame( jobv1t, 'Y' )
395  wantv2t = lsame( jobv2t, 'Y' )
396  colmajor = .NOT. lsame( trans, 'T' )
397 *
398  IF( m .LT. 0 ) THEN
399  info = -6
400  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
401  info = -7
402  ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
403  info = -8
404  ELSE IF( q .GT. p .OR. q .GT. m-p .OR. q .GT. m-q ) THEN
405  info = -8
406  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
407  info = -12
408  ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
409  info = -14
410  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
411  info = -16
412  ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
413  info = -18
414  END IF
415 *
416 * Quick return if Q = 0
417 *
418  IF( info .EQ. 0 .AND. q .EQ. 0 ) THEN
419  lworkmin = 1
420  work(1) = lworkmin
421  RETURN
422  END IF
423 *
424 * Compute workspace
425 *
426  IF( info .EQ. 0 ) THEN
427  iu1cs = 1
428  iu1sn = iu1cs + q
429  iu2cs = iu1sn + q
430  iu2sn = iu2cs + q
431  iv1tcs = iu2sn + q
432  iv1tsn = iv1tcs + q
433  iv2tcs = iv1tsn + q
434  iv2tsn = iv2tcs + q
435  lworkopt = iv2tsn + q - 1
436  lworkmin = lworkopt
437  work(1) = lworkopt
438  IF( lwork .LT. lworkmin .AND. .NOT. lquery ) THEN
439  info = -28
440  END IF
441  END IF
442 *
443  IF( info .NE. 0 ) THEN
444  CALL xerbla( 'DBBCSD', -info )
445  RETURN
446  ELSE IF( lquery ) THEN
447  RETURN
448  END IF
449 *
450 * Get machine constants
451 *
452  eps = dlamch( 'Epsilon' )
453  unfl = dlamch( 'Safe minimum' )
454  tolmul = max( ten, min( hundred, eps**meighth ) )
455  tol = tolmul*eps
456  thresh = max( tol, maxitr*q*q*unfl )
457 *
458 * Test for negligible sines or cosines
459 *
460  DO i = 1, q
461  IF( theta(i) .LT. thresh ) THEN
462  theta(i) = zero
463  ELSE IF( theta(i) .GT. piover2-thresh ) THEN
464  theta(i) = piover2
465  END IF
466  END DO
467  DO i = 1, q-1
468  IF( phi(i) .LT. thresh ) THEN
469  phi(i) = zero
470  ELSE IF( phi(i) .GT. piover2-thresh ) THEN
471  phi(i) = piover2
472  END IF
473  END DO
474 *
475 * Initial deflation
476 *
477  imax = q
478  DO WHILE( imax .GT. 1 )
479  IF( phi(imax-1) .NE. zero ) THEN
480  EXIT
481  END IF
482  imax = imax - 1
483  END DO
484  imin = imax - 1
485  IF ( imin .GT. 1 ) THEN
486  DO WHILE( phi(imin-1) .NE. zero )
487  imin = imin - 1
488  IF ( imin .LE. 1 ) EXIT
489  END DO
490  END IF
491 *
492 * Initialize iteration counter
493 *
494  maxit = maxitr*q*q
495  iter = 0
496 *
497 * Begin main iteration loop
498 *
499  DO WHILE( imax .GT. 1 )
500 *
501 * Compute the matrix entries
502 *
503  b11d(imin) = cos( theta(imin) )
504  b21d(imin) = -sin( theta(imin) )
505  DO i = imin, imax - 1
506  b11e(i) = -sin( theta(i) ) * sin( phi(i) )
507  b11d(i+1) = cos( theta(i+1) ) * cos( phi(i) )
508  b12d(i) = sin( theta(i) ) * cos( phi(i) )
509  b12e(i) = cos( theta(i+1) ) * sin( phi(i) )
510  b21e(i) = -cos( theta(i) ) * sin( phi(i) )
511  b21d(i+1) = -sin( theta(i+1) ) * cos( phi(i) )
512  b22d(i) = cos( theta(i) ) * cos( phi(i) )
513  b22e(i) = -sin( theta(i+1) ) * sin( phi(i) )
514  END DO
515  b12d(imax) = sin( theta(imax) )
516  b22d(imax) = cos( theta(imax) )
517 *
518 * Abort if not converging; otherwise, increment ITER
519 *
520  IF( iter .GT. maxit ) THEN
521  info = 0
522  DO i = 1, q
523  IF( phi(i) .NE. zero )
524  $ info = info + 1
525  END DO
526  RETURN
527  END IF
528 *
529  iter = iter + imax - imin
530 *
531 * Compute shifts
532 *
533  thetamax = theta(imin)
534  thetamin = theta(imin)
535  DO i = imin+1, imax
536  IF( theta(i) > thetamax )
537  $ thetamax = theta(i)
538  IF( theta(i) < thetamin )
539  $ thetamin = theta(i)
540  END DO
541 *
542  IF( thetamax .GT. piover2 - thresh ) THEN
543 *
544 * Zero on diagonals of B11 and B22; induce deflation with a
545 * zero shift
546 *
547  mu = zero
548  nu = one
549 *
550  ELSE IF( thetamin .LT. thresh ) THEN
551 *
552 * Zero on diagonals of B12 and B22; induce deflation with a
553 * zero shift
554 *
555  mu = one
556  nu = zero
557 *
558  ELSE
559 *
560 * Compute shifts for B11 and B21 and use the lesser
561 *
562  CALL dlas2( b11d(imax-1), b11e(imax-1), b11d(imax), sigma11,
563  $ dummy )
564  CALL dlas2( b21d(imax-1), b21e(imax-1), b21d(imax), sigma21,
565  $ dummy )
566 *
567  IF( sigma11 .LE. sigma21 ) THEN
568  mu = sigma11
569  nu = sqrt( one - mu**2 )
570  IF( mu .LT. thresh ) THEN
571  mu = zero
572  nu = one
573  END IF
574  ELSE
575  nu = sigma21
576  mu = sqrt( 1.0 - nu**2 )
577  IF( nu .LT. thresh ) THEN
578  mu = one
579  nu = zero
580  END IF
581  END IF
582  END IF
583 *
584 * Rotate to produce bulges in B11 and B21
585 *
586  IF( mu .LE. nu ) THEN
587  CALL dlartgs( b11d(imin), b11e(imin), mu,
588  $ work(iv1tcs+imin-1), work(iv1tsn+imin-1) )
589  ELSE
590  CALL dlartgs( b21d(imin), b21e(imin), nu,
591  $ work(iv1tcs+imin-1), work(iv1tsn+imin-1) )
592  END IF
593 *
594  temp = work(iv1tcs+imin-1)*b11d(imin) +
595  $ work(iv1tsn+imin-1)*b11e(imin)
596  b11e(imin) = work(iv1tcs+imin-1)*b11e(imin) -
597  $ work(iv1tsn+imin-1)*b11d(imin)
598  b11d(imin) = temp
599  b11bulge = work(iv1tsn+imin-1)*b11d(imin+1)
600  b11d(imin+1) = work(iv1tcs+imin-1)*b11d(imin+1)
601  temp = work(iv1tcs+imin-1)*b21d(imin) +
602  $ work(iv1tsn+imin-1)*b21e(imin)
603  b21e(imin) = work(iv1tcs+imin-1)*b21e(imin) -
604  $ work(iv1tsn+imin-1)*b21d(imin)
605  b21d(imin) = temp
606  b21bulge = work(iv1tsn+imin-1)*b21d(imin+1)
607  b21d(imin+1) = work(iv1tcs+imin-1)*b21d(imin+1)
608 *
609 * Compute THETA(IMIN)
610 *
611  theta( imin ) = atan2( sqrt( b21d(imin)**2+b21bulge**2 ),
612  $ sqrt( b11d(imin)**2+b11bulge**2 ) )
613 *
614 * Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
615 *
616  IF( b11d(imin)**2+b11bulge**2 .GT. thresh**2 ) THEN
617  CALL dlartgp( b11bulge, b11d(imin), work(iu1sn+imin-1),
618  $ work(iu1cs+imin-1), r )
619  ELSE IF( mu .LE. nu ) THEN
620  CALL dlartgs( b11e( imin ), b11d( imin + 1 ), mu,
621  $ work(iu1cs+imin-1), work(iu1sn+imin-1) )
622  ELSE
623  CALL dlartgs( b12d( imin ), b12e( imin ), nu,
624  $ work(iu1cs+imin-1), work(iu1sn+imin-1) )
625  END IF
626  IF( b21d(imin)**2+b21bulge**2 .GT. thresh**2 ) THEN
627  CALL dlartgp( b21bulge, b21d(imin), work(iu2sn+imin-1),
628  $ work(iu2cs+imin-1), r )
629  ELSE IF( nu .LT. mu ) THEN
630  CALL dlartgs( b21e( imin ), b21d( imin + 1 ), nu,
631  $ work(iu2cs+imin-1), work(iu2sn+imin-1) )
632  ELSE
633  CALL dlartgs( b22d(imin), b22e(imin), mu,
634  $ work(iu2cs+imin-1), work(iu2sn+imin-1) )
635  END IF
636  work(iu2cs+imin-1) = -work(iu2cs+imin-1)
637  work(iu2sn+imin-1) = -work(iu2sn+imin-1)
638 *
639  temp = work(iu1cs+imin-1)*b11e(imin) +
640  $ work(iu1sn+imin-1)*b11d(imin+1)
641  b11d(imin+1) = work(iu1cs+imin-1)*b11d(imin+1) -
642  $ work(iu1sn+imin-1)*b11e(imin)
643  b11e(imin) = temp
644  IF( imax .GT. imin+1 ) THEN
645  b11bulge = work(iu1sn+imin-1)*b11e(imin+1)
646  b11e(imin+1) = work(iu1cs+imin-1)*b11e(imin+1)
647  END IF
648  temp = work(iu1cs+imin-1)*b12d(imin) +
649  $ work(iu1sn+imin-1)*b12e(imin)
650  b12e(imin) = work(iu1cs+imin-1)*b12e(imin) -
651  $ work(iu1sn+imin-1)*b12d(imin)
652  b12d(imin) = temp
653  b12bulge = work(iu1sn+imin-1)*b12d(imin+1)
654  b12d(imin+1) = work(iu1cs+imin-1)*b12d(imin+1)
655  temp = work(iu2cs+imin-1)*b21e(imin) +
656  $ work(iu2sn+imin-1)*b21d(imin+1)
657  b21d(imin+1) = work(iu2cs+imin-1)*b21d(imin+1) -
658  $ work(iu2sn+imin-1)*b21e(imin)
659  b21e(imin) = temp
660  IF( imax .GT. imin+1 ) THEN
661  b21bulge = work(iu2sn+imin-1)*b21e(imin+1)
662  b21e(imin+1) = work(iu2cs+imin-1)*b21e(imin+1)
663  END IF
664  temp = work(iu2cs+imin-1)*b22d(imin) +
665  $ work(iu2sn+imin-1)*b22e(imin)
666  b22e(imin) = work(iu2cs+imin-1)*b22e(imin) -
667  $ work(iu2sn+imin-1)*b22d(imin)
668  b22d(imin) = temp
669  b22bulge = work(iu2sn+imin-1)*b22d(imin+1)
670  b22d(imin+1) = work(iu2cs+imin-1)*b22d(imin+1)
671 *
672 * Inner loop: chase bulges from B11(IMIN,IMIN+2),
673 * B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
674 * bottom-right
675 *
676  DO i = imin+1, imax-1
677 *
678 * Compute PHI(I-1)
679 *
680  x1 = sin(theta(i-1))*b11e(i-1) + cos(theta(i-1))*b21e(i-1)
681  x2 = sin(theta(i-1))*b11bulge + cos(theta(i-1))*b21bulge
682  y1 = sin(theta(i-1))*b12d(i-1) + cos(theta(i-1))*b22d(i-1)
683  y2 = sin(theta(i-1))*b12bulge + cos(theta(i-1))*b22bulge
684 *
685  phi(i-1) = atan2( sqrt(x1**2+x2**2), sqrt(y1**2+y2**2) )
686 *
687 * Determine if there are bulges to chase or if a new direct
688 * summand has been reached
689 *
690  restart11 = b11e(i-1)**2 + b11bulge**2 .LE. thresh**2
691  restart21 = b21e(i-1)**2 + b21bulge**2 .LE. thresh**2
692  restart12 = b12d(i-1)**2 + b12bulge**2 .LE. thresh**2
693  restart22 = b22d(i-1)**2 + b22bulge**2 .LE. thresh**2
694 *
695 * If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
696 * B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
697 * chasing by applying the original shift again.
698 *
699  IF( .NOT. restart11 .AND. .NOT. restart21 ) THEN
700  CALL dlartgp( x2, x1, work(iv1tsn+i-1), work(iv1tcs+i-1),
701  $ r )
702  ELSE IF( .NOT. restart11 .AND. restart21 ) THEN
703  CALL dlartgp( b11bulge, b11e(i-1), work(iv1tsn+i-1),
704  $ work(iv1tcs+i-1), r )
705  ELSE IF( restart11 .AND. .NOT. restart21 ) THEN
706  CALL dlartgp( b21bulge, b21e(i-1), work(iv1tsn+i-1),
707  $ work(iv1tcs+i-1), r )
708  ELSE IF( mu .LE. nu ) THEN
709  CALL dlartgs( b11d(i), b11e(i), mu, work(iv1tcs+i-1),
710  $ work(iv1tsn+i-1) )
711  ELSE
712  CALL dlartgs( b21d(i), b21e(i), nu, work(iv1tcs+i-1),
713  $ work(iv1tsn+i-1) )
714  END IF
715  work(iv1tcs+i-1) = -work(iv1tcs+i-1)
716  work(iv1tsn+i-1) = -work(iv1tsn+i-1)
717  IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
718  CALL dlartgp( y2, y1, work(iv2tsn+i-1-1),
719  $ work(iv2tcs+i-1-1), r )
720  ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
721  CALL dlartgp( b12bulge, b12d(i-1), work(iv2tsn+i-1-1),
722  $ work(iv2tcs+i-1-1), r )
723  ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
724  CALL dlartgp( b22bulge, b22d(i-1), work(iv2tsn+i-1-1),
725  $ work(iv2tcs+i-1-1), r )
726  ELSE IF( nu .LT. mu ) THEN
727  CALL dlartgs( b12e(i-1), b12d(i), nu, work(iv2tcs+i-1-1),
728  $ work(iv2tsn+i-1-1) )
729  ELSE
730  CALL dlartgs( b22e(i-1), b22d(i), mu, work(iv2tcs+i-1-1),
731  $ work(iv2tsn+i-1-1) )
732  END IF
733 *
734  temp = work(iv1tcs+i-1)*b11d(i) + work(iv1tsn+i-1)*b11e(i)
735  b11e(i) = work(iv1tcs+i-1)*b11e(i) -
736  $ work(iv1tsn+i-1)*b11d(i)
737  b11d(i) = temp
738  b11bulge = work(iv1tsn+i-1)*b11d(i+1)
739  b11d(i+1) = work(iv1tcs+i-1)*b11d(i+1)
740  temp = work(iv1tcs+i-1)*b21d(i) + work(iv1tsn+i-1)*b21e(i)
741  b21e(i) = work(iv1tcs+i-1)*b21e(i) -
742  $ work(iv1tsn+i-1)*b21d(i)
743  b21d(i) = temp
744  b21bulge = work(iv1tsn+i-1)*b21d(i+1)
745  b21d(i+1) = work(iv1tcs+i-1)*b21d(i+1)
746  temp = work(iv2tcs+i-1-1)*b12e(i-1) +
747  $ work(iv2tsn+i-1-1)*b12d(i)
748  b12d(i) = work(iv2tcs+i-1-1)*b12d(i) -
749  $ work(iv2tsn+i-1-1)*b12e(i-1)
750  b12e(i-1) = temp
751  b12bulge = work(iv2tsn+i-1-1)*b12e(i)
752  b12e(i) = work(iv2tcs+i-1-1)*b12e(i)
753  temp = work(iv2tcs+i-1-1)*b22e(i-1) +
754  $ work(iv2tsn+i-1-1)*b22d(i)
755  b22d(i) = work(iv2tcs+i-1-1)*b22d(i) -
756  $ work(iv2tsn+i-1-1)*b22e(i-1)
757  b22e(i-1) = temp
758  b22bulge = work(iv2tsn+i-1-1)*b22e(i)
759  b22e(i) = work(iv2tcs+i-1-1)*b22e(i)
760 *
761 * Compute THETA(I)
762 *
763  x1 = cos(phi(i-1))*b11d(i) + sin(phi(i-1))*b12e(i-1)
764  x2 = cos(phi(i-1))*b11bulge + sin(phi(i-1))*b12bulge
765  y1 = cos(phi(i-1))*b21d(i) + sin(phi(i-1))*b22e(i-1)
766  y2 = cos(phi(i-1))*b21bulge + sin(phi(i-1))*b22bulge
767 *
768  theta(i) = atan2( sqrt(y1**2+y2**2), sqrt(x1**2+x2**2) )
769 *
770 * Determine if there are bulges to chase or if a new direct
771 * summand has been reached
772 *
773  restart11 = b11d(i)**2 + b11bulge**2 .LE. thresh**2
774  restart12 = b12e(i-1)**2 + b12bulge**2 .LE. thresh**2
775  restart21 = b21d(i)**2 + b21bulge**2 .LE. thresh**2
776  restart22 = b22e(i-1)**2 + b22bulge**2 .LE. thresh**2
777 *
778 * If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
779 * B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
780 * chasing by applying the original shift again.
781 *
782  IF( .NOT. restart11 .AND. .NOT. restart12 ) THEN
783  CALL dlartgp( x2, x1, work(iu1sn+i-1), work(iu1cs+i-1),
784  $ r )
785  ELSE IF( .NOT. restart11 .AND. restart12 ) THEN
786  CALL dlartgp( b11bulge, b11d(i), work(iu1sn+i-1),
787  $ work(iu1cs+i-1), r )
788  ELSE IF( restart11 .AND. .NOT. restart12 ) THEN
789  CALL dlartgp( b12bulge, b12e(i-1), work(iu1sn+i-1),
790  $ work(iu1cs+i-1), r )
791  ELSE IF( mu .LE. nu ) THEN
792  CALL dlartgs( b11e(i), b11d(i+1), mu, work(iu1cs+i-1),
793  $ work(iu1sn+i-1) )
794  ELSE
795  CALL dlartgs( b12d(i), b12e(i), nu, work(iu1cs+i-1),
796  $ work(iu1sn+i-1) )
797  END IF
798  IF( .NOT. restart21 .AND. .NOT. restart22 ) THEN
799  CALL dlartgp( y2, y1, work(iu2sn+i-1), work(iu2cs+i-1),
800  $ r )
801  ELSE IF( .NOT. restart21 .AND. restart22 ) THEN
802  CALL dlartgp( b21bulge, b21d(i), work(iu2sn+i-1),
803  $ work(iu2cs+i-1), r )
804  ELSE IF( restart21 .AND. .NOT. restart22 ) THEN
805  CALL dlartgp( b22bulge, b22e(i-1), work(iu2sn+i-1),
806  $ work(iu2cs+i-1), r )
807  ELSE IF( nu .LT. mu ) THEN
808  CALL dlartgs( b21e(i), b21e(i+1), nu, work(iu2cs+i-1),
809  $ work(iu2sn+i-1) )
810  ELSE
811  CALL dlartgs( b22d(i), b22e(i), mu, work(iu2cs+i-1),
812  $ work(iu2sn+i-1) )
813  END IF
814  work(iu2cs+i-1) = -work(iu2cs+i-1)
815  work(iu2sn+i-1) = -work(iu2sn+i-1)
816 *
817  temp = work(iu1cs+i-1)*b11e(i) + work(iu1sn+i-1)*b11d(i+1)
818  b11d(i+1) = work(iu1cs+i-1)*b11d(i+1) -
819  $ work(iu1sn+i-1)*b11e(i)
820  b11e(i) = temp
821  IF( i .LT. imax - 1 ) THEN
822  b11bulge = work(iu1sn+i-1)*b11e(i+1)
823  b11e(i+1) = work(iu1cs+i-1)*b11e(i+1)
824  END IF
825  temp = work(iu2cs+i-1)*b21e(i) + work(iu2sn+i-1)*b21d(i+1)
826  b21d(i+1) = work(iu2cs+i-1)*b21d(i+1) -
827  $ work(iu2sn+i-1)*b21e(i)
828  b21e(i) = temp
829  IF( i .LT. imax - 1 ) THEN
830  b21bulge = work(iu2sn+i-1)*b21e(i+1)
831  b21e(i+1) = work(iu2cs+i-1)*b21e(i+1)
832  END IF
833  temp = work(iu1cs+i-1)*b12d(i) + work(iu1sn+i-1)*b12e(i)
834  b12e(i) = work(iu1cs+i-1)*b12e(i) - work(iu1sn+i-1)*b12d(i)
835  b12d(i) = temp
836  b12bulge = work(iu1sn+i-1)*b12d(i+1)
837  b12d(i+1) = work(iu1cs+i-1)*b12d(i+1)
838  temp = work(iu2cs+i-1)*b22d(i) + work(iu2sn+i-1)*b22e(i)
839  b22e(i) = work(iu2cs+i-1)*b22e(i) - work(iu2sn+i-1)*b22d(i)
840  b22d(i) = temp
841  b22bulge = work(iu2sn+i-1)*b22d(i+1)
842  b22d(i+1) = work(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, work(iv2tsn+imax-1-1),
863  $ work(iv2tcs+imax-1-1), r )
864  ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
865  CALL dlartgp( b12bulge, b12d(imax-1), work(iv2tsn+imax-1-1),
866  $ work(iv2tcs+imax-1-1), r )
867  ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
868  CALL dlartgp( b22bulge, b22d(imax-1), work(iv2tsn+imax-1-1),
869  $ work(iv2tcs+imax-1-1), r )
870  ELSE IF( nu .LT. mu ) THEN
871  CALL dlartgs( b12e(imax-1), b12d(imax), nu,
872  $ work(iv2tcs+imax-1-1), work(iv2tsn+imax-1-1) )
873  ELSE
874  CALL dlartgs( b22e(imax-1), b22d(imax), mu,
875  $ work(iv2tcs+imax-1-1), work(iv2tsn+imax-1-1) )
876  END IF
877 *
878  temp = work(iv2tcs+imax-1-1)*b12e(imax-1) +
879  $ work(iv2tsn+imax-1-1)*b12d(imax)
880  b12d(imax) = work(iv2tcs+imax-1-1)*b12d(imax) -
881  $ work(iv2tsn+imax-1-1)*b12e(imax-1)
882  b12e(imax-1) = temp
883  temp = work(iv2tcs+imax-1-1)*b22e(imax-1) +
884  $ work(iv2tsn+imax-1-1)*b22d(imax)
885  b22d(imax) = work(iv2tcs+imax-1-1)*b22d(imax) -
886  $ work(iv2tsn+imax-1-1)*b22e(imax-1)
887  b22e(imax-1) = temp
888 *
889 * Update singular vectors
890 *
891  IF( wantu1 ) THEN
892  IF( colmajor ) THEN
893  CALL dlasr( 'R', 'V', 'F', p, imax-imin+1,
894  $ work(iu1cs+imin-1), work(iu1sn+imin-1),
895  $ u1(1,imin), ldu1 )
896  ELSE
897  CALL dlasr( 'L', 'V', 'F', imax-imin+1, p,
898  $ work(iu1cs+imin-1), work(iu1sn+imin-1),
899  $ u1(imin,1), ldu1 )
900  END IF
901  END IF
902  IF( wantu2 ) THEN
903  IF( colmajor ) THEN
904  CALL dlasr( 'R', 'V', 'F', m-p, imax-imin+1,
905  $ work(iu2cs+imin-1), work(iu2sn+imin-1),
906  $ u2(1,imin), ldu2 )
907  ELSE
908  CALL dlasr( 'L', 'V', 'F', imax-imin+1, m-p,
909  $ work(iu2cs+imin-1), work(iu2sn+imin-1),
910  $ u2(imin,1), ldu2 )
911  END IF
912  END IF
913  IF( wantv1t ) THEN
914  IF( colmajor ) THEN
915  CALL dlasr( 'L', 'V', 'F', imax-imin+1, q,
916  $ work(iv1tcs+imin-1), work(iv1tsn+imin-1),
917  $ v1t(imin,1), ldv1t )
918  ELSE
919  CALL dlasr( 'R', 'V', 'F', q, imax-imin+1,
920  $ work(iv1tcs+imin-1), work(iv1tsn+imin-1),
921  $ v1t(1,imin), ldv1t )
922  END IF
923  END IF
924  IF( wantv2t ) THEN
925  IF( colmajor ) THEN
926  CALL dlasr( 'L', 'V', 'F', imax-imin+1, m-q,
927  $ work(iv2tcs+imin-1), work(iv2tsn+imin-1),
928  $ v2t(imin,1), ldv2t )
929  ELSE
930  CALL dlasr( 'R', 'V', 'F', m-q, imax-imin+1,
931  $ work(iv2tcs+imin-1), work(iv2tsn+imin-1),
932  $ v2t(1,imin), ldv2t )
933  END IF
934  END IF
935 *
936 * Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
937 *
938  IF( b11e(imax-1)+b21e(imax-1) .GT. 0 ) THEN
939  b11d(imax) = -b11d(imax)
940  b21d(imax) = -b21d(imax)
941  IF( wantv1t ) THEN
942  IF( colmajor ) THEN
943  CALL dscal( q, negone, v1t(imax,1), ldv1t )
944  ELSE
945  CALL dscal( q, negone, v1t(1,imax), 1 )
946  END IF
947  END IF
948  END IF
949 *
950 * Compute THETA(IMAX)
951 *
952  x1 = cos(phi(imax-1))*b11d(imax) +
953  $ sin(phi(imax-1))*b12e(imax-1)
954  y1 = cos(phi(imax-1))*b21d(imax) +
955  $ sin(phi(imax-1))*b22e(imax-1)
956 *
957  theta(imax) = atan2( abs(y1), abs(x1) )
958 *
959 * Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
960 * and B22(IMAX,IMAX-1)
961 *
962  IF( b11d(imax)+b12e(imax-1) .LT. 0 ) THEN
963  b12d(imax) = -b12d(imax)
964  IF( wantu1 ) THEN
965  IF( colmajor ) THEN
966  CALL dscal( p, negone, u1(1,imax), 1 )
967  ELSE
968  CALL dscal( p, negone, u1(imax,1), ldu1 )
969  END IF
970  END IF
971  END IF
972  IF( b21d(imax)+b22e(imax-1) .GT. 0 ) THEN
973  b22d(imax) = -b22d(imax)
974  IF( wantu2 ) THEN
975  IF( colmajor ) THEN
976  CALL dscal( m-p, negone, u2(1,imax), 1 )
977  ELSE
978  CALL dscal( m-p, negone, u2(imax,1), ldu2 )
979  END IF
980  END IF
981  END IF
982 *
983 * Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
984 *
985  IF( b12d(imax)+b22d(imax) .LT. 0 ) THEN
986  IF( wantv2t ) THEN
987  IF( colmajor ) THEN
988  CALL dscal( m-q, negone, v2t(imax,1), ldv2t )
989  ELSE
990  CALL dscal( m-q, negone, v2t(1,imax), 1 )
991  END IF
992  END IF
993  END IF
994 *
995 * Test for negligible sines or cosines
996 *
997  DO i = imin, imax
998  IF( theta(i) .LT. thresh ) THEN
999  theta(i) = zero
1000  ELSE IF( theta(i) .GT. piover2-thresh ) THEN
1001  theta(i) = piover2
1002  END IF
1003  END DO
1004  DO i = imin, imax-1
1005  IF( phi(i) .LT. thresh ) THEN
1006  phi(i) = zero
1007  ELSE IF( phi(i) .GT. piover2-thresh ) THEN
1008  phi(i) = piover2
1009  END IF
1010  END DO
1011 *
1012 * Deflate
1013 *
1014  IF (imax .GT. 1) THEN
1015  DO WHILE( phi(imax-1) .EQ. zero )
1016  imax = imax - 1
1017  IF (imax .LE. 1) EXIT
1018  END DO
1019  END IF
1020  IF( imin .GT. imax - 1 )
1021  $ imin = imax - 1
1022  IF (imin .GT. 1) THEN
1023  DO WHILE (phi(imin-1) .NE. zero)
1024  imin = imin - 1
1025  IF (imin .LE. 1) EXIT
1026  END DO
1027  END IF
1028 *
1029 * Repeat main iteration loop
1030 *
1031  END DO
1032 *
1033 * Postprocessing: order THETA from least to greatest
1034 *
1035  DO i = 1, q
1036 *
1037  mini = i
1038  thetamin = theta(i)
1039  DO j = i+1, q
1040  IF( theta(j) .LT. thetamin ) THEN
1041  mini = j
1042  thetamin = theta(j)
1043  END IF
1044  END DO
1045 *
1046  IF( mini .NE. i ) THEN
1047  theta(mini) = theta(i)
1048  theta(i) = thetamin
1049  IF( colmajor ) THEN
1050  IF( wantu1 )
1051  $ CALL dswap( p, u1(1,i), 1, u1(1,mini), 1 )
1052  IF( wantu2 )
1053  $ CALL dswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1054  IF( wantv1t )
1055  $ CALL dswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1056  IF( wantv2t )
1057  $ CALL dswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1058  $ ldv2t )
1059  ELSE
1060  IF( wantu1 )
1061  $ CALL dswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1062  IF( wantu2 )
1063  $ CALL dswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1064  IF( wantv1t )
1065  $ CALL dswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1066  IF( wantv2t )
1067  $ CALL dswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
1068  END IF
1069  END IF
1070 *
1071  END DO
1072 *
1073  RETURN
1074 *
1075 * End of DBBCSD
1076 *
1077  END
1078 
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: dlas2.f:107
subroutine dlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: dlasr.f:199
subroutine dlartgp(F, G, CS, SN, R)
DLARTGP generates a plane rotation so that the diagonal is nonnegative.
Definition: dlartgp.f:95
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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:90
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dbbcsd(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, WORK, LWORK, INFO)
DBBCSD
Definition: dbbcsd.f:332