LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
cbbcsd.f
Go to the documentation of this file.
1 *> \brief \b CBBCSD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CBBCSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cbbcsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cbbcsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cbbcsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CBBCSD( 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 * REAL B11D( * ), B11E( * ), B12D( * ), B12E( * ),
32 * $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
33 * $ PHI( * ), THETA( * ), RWORK( * )
34 * COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
35 * $ V2T( LDV2T, * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> CBBCSD 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 CUNCSD 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 REAL 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 REAL 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 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 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 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 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 REAL array, dimension (Q)
209 *> When CBBCSD converges, B11D contains the cosines of THETA(1),
210 *> ..., THETA(Q). If CBBCSD 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 REAL array, dimension (Q-1)
218 *> When CBBCSD converges, B11E contains zeros. If CBBCSD 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 REAL array, dimension (Q)
226 *> When CBBCSD converges, B12D contains the negative sines of
227 *> THETA(1), ..., THETA(Q). If CBBCSD 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 REAL array, dimension (Q-1)
235 *> When CBBCSD converges, B12E contains zeros. If CBBCSD 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 CBBCSD 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 REAL, 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 complexOTHERcomputational
326 *
327 * =====================================================================
328  SUBROUTINE cbbcsd( 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, RWORK, LRWORK, 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, LRWORK, M, P, Q
340 * ..
341 * .. Array Arguments ..
342  REAL B11D( * ), B11E( * ), B12D( * ), B12E( * ),
343  $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
344  $ PHI( * ), THETA( * ), RWORK( * )
345  COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
346  $ V2T( LDV2T, * )
347 * ..
348 *
349 * ===================================================================
350 *
351 * .. Parameters ..
352  INTEGER MAXITR
353  PARAMETER ( MAXITR = 6 )
354  real hundred, meighth, one, ten, zero
355  parameter( hundred = 100.0e0, meighth = -0.125e0,
356  $ one = 1.0e0, ten = 10.0e0, zero = 0.0e0 )
357  COMPLEX NEGONECOMPLEX
358  parameter( negonecomplex = (-1.0e0,0.0e0) )
359  REAL PIOVER2
360  parameter( piover2 = 1.57079632679489661923132169163975144210e0 )
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  $ LRWORKMIN, LRWORKOPT, MAXIT, MINI
369  REAL 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 clasr, cscal, cswap, slartgp, slartgs, slas2,
376  $ xerbla
377 * ..
378 * .. External Functions ..
379  REAL SLAMCH
380  LOGICAL LSAME
381  EXTERNAL LSAME, SLAMCH
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 = lrwork .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  lrworkmin = 1
420  rwork(1) = lrworkmin
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  lrworkopt = iv2tsn + q - 1
436  lrworkmin = lrworkopt
437  rwork(1) = lrworkopt
438  IF( lrwork .LT. lrworkmin .AND. .NOT. lquery ) THEN
439  info = -28
440  END IF
441  END IF
442 *
443  IF( info .NE. 0 ) THEN
444  CALL xerbla( 'CBBCSD', -info )
445  RETURN
446  ELSE IF( lquery ) THEN
447  RETURN
448  END IF
449 *
450 * Get machine constants
451 *
452  eps = slamch( 'Epsilon' )
453  unfl = slamch( '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 slas2( b11d(imax-1), b11e(imax-1), b11d(imax), sigma11,
563  $ dummy )
564  CALL slas2( 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 slartgs( b11d(imin), b11e(imin), mu,
588  $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
589  ELSE
590  CALL slartgs( b21d(imin), b21e(imin), nu,
591  $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
592  END IF
593 *
594  temp = rwork(iv1tcs+imin-1)*b11d(imin) +
595  $ rwork(iv1tsn+imin-1)*b11e(imin)
596  b11e(imin) = rwork(iv1tcs+imin-1)*b11e(imin) -
597  $ rwork(iv1tsn+imin-1)*b11d(imin)
598  b11d(imin) = temp
599  b11bulge = rwork(iv1tsn+imin-1)*b11d(imin+1)
600  b11d(imin+1) = rwork(iv1tcs+imin-1)*b11d(imin+1)
601  temp = rwork(iv1tcs+imin-1)*b21d(imin) +
602  $ rwork(iv1tsn+imin-1)*b21e(imin)
603  b21e(imin) = rwork(iv1tcs+imin-1)*b21e(imin) -
604  $ rwork(iv1tsn+imin-1)*b21d(imin)
605  b21d(imin) = temp
606  b21bulge = rwork(iv1tsn+imin-1)*b21d(imin+1)
607  b21d(imin+1) = rwork(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 slartgp( b11bulge, b11d(imin), rwork(iu1sn+imin-1),
618  $ rwork(iu1cs+imin-1), r )
619  ELSE IF( mu .LE. nu ) THEN
620  CALL slartgs( b11e( imin ), b11d( imin + 1 ), mu,
621  $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
622  ELSE
623  CALL slartgs( b12d( imin ), b12e( imin ), nu,
624  $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
625  END IF
626  IF( b21d(imin)**2+b21bulge**2 .GT. thresh**2 ) THEN
627  CALL slartgp( b21bulge, b21d(imin), rwork(iu2sn+imin-1),
628  $ rwork(iu2cs+imin-1), r )
629  ELSE IF( nu .LT. mu ) THEN
630  CALL slartgs( b21e( imin ), b21d( imin + 1 ), nu,
631  $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
632  ELSE
633  CALL slartgs( b22d(imin), b22e(imin), mu,
634  $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
635  END IF
636  rwork(iu2cs+imin-1) = -rwork(iu2cs+imin-1)
637  rwork(iu2sn+imin-1) = -rwork(iu2sn+imin-1)
638 *
639  temp = rwork(iu1cs+imin-1)*b11e(imin) +
640  $ rwork(iu1sn+imin-1)*b11d(imin+1)
641  b11d(imin+1) = rwork(iu1cs+imin-1)*b11d(imin+1) -
642  $ rwork(iu1sn+imin-1)*b11e(imin)
643  b11e(imin) = temp
644  IF( imax .GT. imin+1 ) THEN
645  b11bulge = rwork(iu1sn+imin-1)*b11e(imin+1)
646  b11e(imin+1) = rwork(iu1cs+imin-1)*b11e(imin+1)
647  END IF
648  temp = rwork(iu1cs+imin-1)*b12d(imin) +
649  $ rwork(iu1sn+imin-1)*b12e(imin)
650  b12e(imin) = rwork(iu1cs+imin-1)*b12e(imin) -
651  $ rwork(iu1sn+imin-1)*b12d(imin)
652  b12d(imin) = temp
653  b12bulge = rwork(iu1sn+imin-1)*b12d(imin+1)
654  b12d(imin+1) = rwork(iu1cs+imin-1)*b12d(imin+1)
655  temp = rwork(iu2cs+imin-1)*b21e(imin) +
656  $ rwork(iu2sn+imin-1)*b21d(imin+1)
657  b21d(imin+1) = rwork(iu2cs+imin-1)*b21d(imin+1) -
658  $ rwork(iu2sn+imin-1)*b21e(imin)
659  b21e(imin) = temp
660  IF( imax .GT. imin+1 ) THEN
661  b21bulge = rwork(iu2sn+imin-1)*b21e(imin+1)
662  b21e(imin+1) = rwork(iu2cs+imin-1)*b21e(imin+1)
663  END IF
664  temp = rwork(iu2cs+imin-1)*b22d(imin) +
665  $ rwork(iu2sn+imin-1)*b22e(imin)
666  b22e(imin) = rwork(iu2cs+imin-1)*b22e(imin) -
667  $ rwork(iu2sn+imin-1)*b22d(imin)
668  b22d(imin) = temp
669  b22bulge = rwork(iu2sn+imin-1)*b22d(imin+1)
670  b22d(imin+1) = rwork(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 slartgp( x2, x1, rwork(iv1tsn+i-1),
701  $ rwork(iv1tcs+i-1), r )
702  ELSE IF( .NOT. restart11 .AND. restart21 ) THEN
703  CALL slartgp( b11bulge, b11e(i-1), rwork(iv1tsn+i-1),
704  $ rwork(iv1tcs+i-1), r )
705  ELSE IF( restart11 .AND. .NOT. restart21 ) THEN
706  CALL slartgp( b21bulge, b21e(i-1), rwork(iv1tsn+i-1),
707  $ rwork(iv1tcs+i-1), r )
708  ELSE IF( mu .LE. nu ) THEN
709  CALL slartgs( b11d(i), b11e(i), mu, rwork(iv1tcs+i-1),
710  $ rwork(iv1tsn+i-1) )
711  ELSE
712  CALL slartgs( b21d(i), b21e(i), nu, rwork(iv1tcs+i-1),
713  $ rwork(iv1tsn+i-1) )
714  END IF
715  rwork(iv1tcs+i-1) = -rwork(iv1tcs+i-1)
716  rwork(iv1tsn+i-1) = -rwork(iv1tsn+i-1)
717  IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
718  CALL slartgp( y2, y1, rwork(iv2tsn+i-1-1),
719  $ rwork(iv2tcs+i-1-1), r )
720  ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
721  CALL slartgp( b12bulge, b12d(i-1), rwork(iv2tsn+i-1-1),
722  $ rwork(iv2tcs+i-1-1), r )
723  ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
724  CALL slartgp( b22bulge, b22d(i-1), rwork(iv2tsn+i-1-1),
725  $ rwork(iv2tcs+i-1-1), r )
726  ELSE IF( nu .LT. mu ) THEN
727  CALL slartgs( b12e(i-1), b12d(i), nu,
728  $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
729  ELSE
730  CALL slartgs( b22e(i-1), b22d(i), mu,
731  $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
732  END IF
733 *
734  temp = rwork(iv1tcs+i-1)*b11d(i) + rwork(iv1tsn+i-1)*b11e(i)
735  b11e(i) = rwork(iv1tcs+i-1)*b11e(i) -
736  $ rwork(iv1tsn+i-1)*b11d(i)
737  b11d(i) = temp
738  b11bulge = rwork(iv1tsn+i-1)*b11d(i+1)
739  b11d(i+1) = rwork(iv1tcs+i-1)*b11d(i+1)
740  temp = rwork(iv1tcs+i-1)*b21d(i) + rwork(iv1tsn+i-1)*b21e(i)
741  b21e(i) = rwork(iv1tcs+i-1)*b21e(i) -
742  $ rwork(iv1tsn+i-1)*b21d(i)
743  b21d(i) = temp
744  b21bulge = rwork(iv1tsn+i-1)*b21d(i+1)
745  b21d(i+1) = rwork(iv1tcs+i-1)*b21d(i+1)
746  temp = rwork(iv2tcs+i-1-1)*b12e(i-1) +
747  $ rwork(iv2tsn+i-1-1)*b12d(i)
748  b12d(i) = rwork(iv2tcs+i-1-1)*b12d(i) -
749  $ rwork(iv2tsn+i-1-1)*b12e(i-1)
750  b12e(i-1) = temp
751  b12bulge = rwork(iv2tsn+i-1-1)*b12e(i)
752  b12e(i) = rwork(iv2tcs+i-1-1)*b12e(i)
753  temp = rwork(iv2tcs+i-1-1)*b22e(i-1) +
754  $ rwork(iv2tsn+i-1-1)*b22d(i)
755  b22d(i) = rwork(iv2tcs+i-1-1)*b22d(i) -
756  $ rwork(iv2tsn+i-1-1)*b22e(i-1)
757  b22e(i-1) = temp
758  b22bulge = rwork(iv2tsn+i-1-1)*b22e(i)
759  b22e(i) = rwork(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 slartgp( x2, x1, rwork(iu1sn+i-1), rwork(iu1cs+i-1),
784  $ r )
785  ELSE IF( .NOT. restart11 .AND. restart12 ) THEN
786  CALL slartgp( b11bulge, b11d(i), rwork(iu1sn+i-1),
787  $ rwork(iu1cs+i-1), r )
788  ELSE IF( restart11 .AND. .NOT. restart12 ) THEN
789  CALL slartgp( b12bulge, b12e(i-1), rwork(iu1sn+i-1),
790  $ rwork(iu1cs+i-1), r )
791  ELSE IF( mu .LE. nu ) THEN
792  CALL slartgs( b11e(i), b11d(i+1), mu, rwork(iu1cs+i-1),
793  $ rwork(iu1sn+i-1) )
794  ELSE
795  CALL slartgs( b12d(i), b12e(i), nu, rwork(iu1cs+i-1),
796  $ rwork(iu1sn+i-1) )
797  END IF
798  IF( .NOT. restart21 .AND. .NOT. restart22 ) THEN
799  CALL slartgp( y2, y1, rwork(iu2sn+i-1), rwork(iu2cs+i-1),
800  $ r )
801  ELSE IF( .NOT. restart21 .AND. restart22 ) THEN
802  CALL slartgp( b21bulge, b21d(i), rwork(iu2sn+i-1),
803  $ rwork(iu2cs+i-1), r )
804  ELSE IF( restart21 .AND. .NOT. restart22 ) THEN
805  CALL slartgp( b22bulge, b22e(i-1), rwork(iu2sn+i-1),
806  $ rwork(iu2cs+i-1), r )
807  ELSE IF( nu .LT. mu ) THEN
808  CALL slartgs( b21e(i), b21e(i+1), nu, rwork(iu2cs+i-1),
809  $ rwork(iu2sn+i-1) )
810  ELSE
811  CALL slartgs( b22d(i), b22e(i), mu, rwork(iu2cs+i-1),
812  $ rwork(iu2sn+i-1) )
813  END IF
814  rwork(iu2cs+i-1) = -rwork(iu2cs+i-1)
815  rwork(iu2sn+i-1) = -rwork(iu2sn+i-1)
816 *
817  temp = rwork(iu1cs+i-1)*b11e(i) + rwork(iu1sn+i-1)*b11d(i+1)
818  b11d(i+1) = rwork(iu1cs+i-1)*b11d(i+1) -
819  $ rwork(iu1sn+i-1)*b11e(i)
820  b11e(i) = temp
821  IF( i .LT. imax - 1 ) THEN
822  b11bulge = rwork(iu1sn+i-1)*b11e(i+1)
823  b11e(i+1) = rwork(iu1cs+i-1)*b11e(i+1)
824  END IF
825  temp = rwork(iu2cs+i-1)*b21e(i) + rwork(iu2sn+i-1)*b21d(i+1)
826  b21d(i+1) = rwork(iu2cs+i-1)*b21d(i+1) -
827  $ rwork(iu2sn+i-1)*b21e(i)
828  b21e(i) = temp
829  IF( i .LT. imax - 1 ) THEN
830  b21bulge = rwork(iu2sn+i-1)*b21e(i+1)
831  b21e(i+1) = rwork(iu2cs+i-1)*b21e(i+1)
832  END IF
833  temp = rwork(iu1cs+i-1)*b12d(i) + rwork(iu1sn+i-1)*b12e(i)
834  b12e(i) = rwork(iu1cs+i-1)*b12e(i) -
835  $ rwork(iu1sn+i-1)*b12d(i)
836  b12d(i) = temp
837  b12bulge = rwork(iu1sn+i-1)*b12d(i+1)
838  b12d(i+1) = rwork(iu1cs+i-1)*b12d(i+1)
839  temp = rwork(iu2cs+i-1)*b22d(i) + rwork(iu2sn+i-1)*b22e(i)
840  b22e(i) = rwork(iu2cs+i-1)*b22e(i) -
841  $ rwork(iu2sn+i-1)*b22d(i)
842  b22d(i) = temp
843  b22bulge = rwork(iu2sn+i-1)*b22d(i+1)
844  b22d(i+1) = rwork(iu2cs+i-1)*b22d(i+1)
845 *
846  END DO
847 *
848 * Compute PHI(IMAX-1)
849 *
850  x1 = sin(theta(imax-1))*b11e(imax-1) +
851  $ cos(theta(imax-1))*b21e(imax-1)
852  y1 = sin(theta(imax-1))*b12d(imax-1) +
853  $ cos(theta(imax-1))*b22d(imax-1)
854  y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
855 *
856  phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
857 *
858 * Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
859 *
860  restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
861  restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
862 *
863  IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
864  CALL slartgp( y2, y1, rwork(iv2tsn+imax-1-1),
865  $ rwork(iv2tcs+imax-1-1), r )
866  ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
867  CALL slartgp( b12bulge, b12d(imax-1),
868  $ rwork(iv2tsn+imax-1-1),
869  $ rwork(iv2tcs+imax-1-1), r )
870  ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
871  CALL slartgp( b22bulge, b22d(imax-1),
872  $ rwork(iv2tsn+imax-1-1),
873  $ rwork(iv2tcs+imax-1-1), r )
874  ELSE IF( nu .LT. mu ) THEN
875  CALL slartgs( b12e(imax-1), b12d(imax), nu,
876  $ rwork(iv2tcs+imax-1-1),
877  $ rwork(iv2tsn+imax-1-1) )
878  ELSE
879  CALL slartgs( b22e(imax-1), b22d(imax), mu,
880  $ rwork(iv2tcs+imax-1-1),
881  $ rwork(iv2tsn+imax-1-1) )
882  END IF
883 *
884  temp = rwork(iv2tcs+imax-1-1)*b12e(imax-1) +
885  $ rwork(iv2tsn+imax-1-1)*b12d(imax)
886  b12d(imax) = rwork(iv2tcs+imax-1-1)*b12d(imax) -
887  $ rwork(iv2tsn+imax-1-1)*b12e(imax-1)
888  b12e(imax-1) = temp
889  temp = rwork(iv2tcs+imax-1-1)*b22e(imax-1) +
890  $ rwork(iv2tsn+imax-1-1)*b22d(imax)
891  b22d(imax) = rwork(iv2tcs+imax-1-1)*b22d(imax) -
892  $ rwork(iv2tsn+imax-1-1)*b22e(imax-1)
893  b22e(imax-1) = temp
894 *
895 * Update singular vectors
896 *
897  IF( wantu1 ) THEN
898  IF( colmajor ) THEN
899  CALL clasr( 'R', 'V', 'F', p, imax-imin+1,
900  $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
901  $ u1(1,imin), ldu1 )
902  ELSE
903  CALL clasr( 'L', 'V', 'F', imax-imin+1, p,
904  $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
905  $ u1(imin,1), ldu1 )
906  END IF
907  END IF
908  IF( wantu2 ) THEN
909  IF( colmajor ) THEN
910  CALL clasr( 'R', 'V', 'F', m-p, imax-imin+1,
911  $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
912  $ u2(1,imin), ldu2 )
913  ELSE
914  CALL clasr( 'L', 'V', 'F', imax-imin+1, m-p,
915  $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
916  $ u2(imin,1), ldu2 )
917  END IF
918  END IF
919  IF( wantv1t ) THEN
920  IF( colmajor ) THEN
921  CALL clasr( 'L', 'V', 'F', imax-imin+1, q,
922  $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
923  $ v1t(imin,1), ldv1t )
924  ELSE
925  CALL clasr( 'R', 'V', 'F', q, imax-imin+1,
926  $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
927  $ v1t(1,imin), ldv1t )
928  END IF
929  END IF
930  IF( wantv2t ) THEN
931  IF( colmajor ) THEN
932  CALL clasr( 'L', 'V', 'F', imax-imin+1, m-q,
933  $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
934  $ v2t(imin,1), ldv2t )
935  ELSE
936  CALL clasr( 'R', 'V', 'F', m-q, imax-imin+1,
937  $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
938  $ v2t(1,imin), ldv2t )
939  END IF
940  END IF
941 *
942 * Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
943 *
944  IF( b11e(imax-1)+b21e(imax-1) .GT. 0 ) THEN
945  b11d(imax) = -b11d(imax)
946  b21d(imax) = -b21d(imax)
947  IF( wantv1t ) THEN
948  IF( colmajor ) THEN
949  CALL cscal( q, negonecomplex, v1t(imax,1), ldv1t )
950  ELSE
951  CALL cscal( q, negonecomplex, v1t(1,imax), 1 )
952  END IF
953  END IF
954  END IF
955 *
956 * Compute THETA(IMAX)
957 *
958  x1 = cos(phi(imax-1))*b11d(imax) +
959  $ sin(phi(imax-1))*b12e(imax-1)
960  y1 = cos(phi(imax-1))*b21d(imax) +
961  $ sin(phi(imax-1))*b22e(imax-1)
962 *
963  theta(imax) = atan2( abs(y1), abs(x1) )
964 *
965 * Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
966 * and B22(IMAX,IMAX-1)
967 *
968  IF( b11d(imax)+b12e(imax-1) .LT. 0 ) THEN
969  b12d(imax) = -b12d(imax)
970  IF( wantu1 ) THEN
971  IF( colmajor ) THEN
972  CALL cscal( p, negonecomplex, u1(1,imax), 1 )
973  ELSE
974  CALL cscal( p, negonecomplex, u1(imax,1), ldu1 )
975  END IF
976  END IF
977  END IF
978  IF( b21d(imax)+b22e(imax-1) .GT. 0 ) THEN
979  b22d(imax) = -b22d(imax)
980  IF( wantu2 ) THEN
981  IF( colmajor ) THEN
982  CALL cscal( m-p, negonecomplex, u2(1,imax), 1 )
983  ELSE
984  CALL cscal( m-p, negonecomplex, u2(imax,1), ldu2 )
985  END IF
986  END IF
987  END IF
988 *
989 * Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
990 *
991  IF( b12d(imax)+b22d(imax) .LT. 0 ) THEN
992  IF( wantv2t ) THEN
993  IF( colmajor ) THEN
994  CALL cscal( m-q, negonecomplex, v2t(imax,1), ldv2t )
995  ELSE
996  CALL cscal( m-q, negonecomplex, v2t(1,imax), 1 )
997  END IF
998  END IF
999  END IF
1000 *
1001 * Test for negligible sines or cosines
1002 *
1003  DO i = imin, imax
1004  IF( theta(i) .LT. thresh ) THEN
1005  theta(i) = zero
1006  ELSE IF( theta(i) .GT. piover2-thresh ) THEN
1007  theta(i) = piover2
1008  END IF
1009  END DO
1010  DO i = imin, imax-1
1011  IF( phi(i) .LT. thresh ) THEN
1012  phi(i) = zero
1013  ELSE IF( phi(i) .GT. piover2-thresh ) THEN
1014  phi(i) = piover2
1015  END IF
1016  END DO
1017 *
1018 * Deflate
1019 *
1020  IF (imax .GT. 1) THEN
1021  DO WHILE( phi(imax-1) .EQ. zero )
1022  imax = imax - 1
1023  IF (imax .LE. 1) EXIT
1024  END DO
1025  END IF
1026  IF( imin .GT. imax - 1 )
1027  $ imin = imax - 1
1028  IF (imin .GT. 1) THEN
1029  DO WHILE (phi(imin-1) .NE. zero)
1030  imin = imin - 1
1031  IF (imin .LE. 1) EXIT
1032  END DO
1033  END IF
1034 *
1035 * Repeat main iteration loop
1036 *
1037  END DO
1038 *
1039 * Postprocessing: order THETA from least to greatest
1040 *
1041  DO i = 1, q
1042 *
1043  mini = i
1044  thetamin = theta(i)
1045  DO j = i+1, q
1046  IF( theta(j) .LT. thetamin ) THEN
1047  mini = j
1048  thetamin = theta(j)
1049  END IF
1050  END DO
1051 *
1052  IF( mini .NE. i ) THEN
1053  theta(mini) = theta(i)
1054  theta(i) = thetamin
1055  IF( colmajor ) THEN
1056  IF( wantu1 )
1057  $ CALL cswap( p, u1(1,i), 1, u1(1,mini), 1 )
1058  IF( wantu2 )
1059  $ CALL cswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1060  IF( wantv1t )
1061  $ CALL cswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1062  IF( wantv2t )
1063  $ CALL cswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1064  $ ldv2t )
1065  ELSE
1066  IF( wantu1 )
1067  $ CALL cswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1068  IF( wantu2 )
1069  $ CALL cswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1070  IF( wantv1t )
1071  $ CALL cswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1072  IF( wantv2t )
1073  $ CALL cswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
1074  END IF
1075  END IF
1076 *
1077  END DO
1078 *
1079  RETURN
1080 *
1081 * End of CBBCSD
1082 *
1083  END
1084 
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: slas2.f:107
subroutine slartgp(F, G, CS, SN, R)
SLARTGP generates a plane rotation so that the diagonal is nonnegative.
Definition: slartgp.f:95
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slartgs(X, Y, SIGMA, CS, SN)
SLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bid...
Definition: slartgs.f:90
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine clasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
CLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: clasr.f:200
subroutine cbbcsd(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)
CBBCSD
Definition: cbbcsd.f:332