LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> Download ZBBCSD + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zbbcsd.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zbbcsd.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zbbcsd.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
20* THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T,
21* V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E,
22* B22D, B22E, RWORK, LRWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
26* INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
30* $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
31* $ PHI( * ), THETA( * ), RWORK( * )
32* COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
33* $ V2T( LDV2T, * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> ZBBCSD computes the CS decomposition of a unitary matrix in
43*> bidiagonal-block form,
44*>
45*>
46*> [ B11 | B12 0 0 ]
47*> [ 0 | 0 -I 0 ]
48*> X = [----------------]
49*> [ B21 | B22 0 0 ]
50*> [ 0 | 0 0 I ]
51*>
52*> [ C | -S 0 0 ]
53*> [ U1 | ] [ 0 | 0 -I 0 ] [ V1 | ]**H
54*> = [---------] [---------------] [---------] .
55*> [ | U2 ] [ S | C 0 0 ] [ | V2 ]
56*> [ 0 | 0 0 I ]
57*>
58*> X is M-by-M, its top-left block is P-by-Q, and Q must be no larger
59*> than P, M-P, or M-Q. (If Q is not the smallest index, then X must be
60*> transposed and/or permuted. This can be done in constant time using
61*> the TRANS and SIGNS options. See ZUNCSD for details.)
62*>
63*> The bidiagonal matrices B11, B12, B21, and B22 are represented
64*> implicitly by angles THETA(1:Q) and PHI(1:Q-1).
65*>
66*> The unitary matrices U1, U2, V1T, and V2T are input/output.
67*> The input matrices are pre- or post-multiplied by the appropriate
68*> singular vector matrices.
69*> \endverbatim
70*
71* Arguments:
72* ==========
73*
74*> \param[in] JOBU1
75*> \verbatim
76*> JOBU1 is CHARACTER
77*> = 'Y': U1 is updated;
78*> otherwise: U1 is not updated.
79*> \endverbatim
80*>
81*> \param[in] JOBU2
82*> \verbatim
83*> JOBU2 is CHARACTER
84*> = 'Y': U2 is updated;
85*> otherwise: U2 is not updated.
86*> \endverbatim
87*>
88*> \param[in] JOBV1T
89*> \verbatim
90*> JOBV1T is CHARACTER
91*> = 'Y': V1T is updated;
92*> otherwise: V1T is not updated.
93*> \endverbatim
94*>
95*> \param[in] JOBV2T
96*> \verbatim
97*> JOBV2T is CHARACTER
98*> = 'Y': V2T is updated;
99*> otherwise: V2T is not updated.
100*> \endverbatim
101*>
102*> \param[in] TRANS
103*> \verbatim
104*> TRANS is CHARACTER
105*> = 'T': X, U1, U2, V1T, and V2T are stored in row-major
106*> order;
107*> otherwise: X, U1, U2, V1T, and V2T are stored in column-
108*> major order.
109*> \endverbatim
110*>
111*> \param[in] M
112*> \verbatim
113*> M is INTEGER
114*> The number of rows and columns in X, the unitary matrix in
115*> bidiagonal-block form.
116*> \endverbatim
117*>
118*> \param[in] P
119*> \verbatim
120*> P is INTEGER
121*> The number of rows in the top-left block of X. 0 <= P <= M.
122*> \endverbatim
123*>
124*> \param[in] Q
125*> \verbatim
126*> Q is INTEGER
127*> The number of columns in the top-left block of X.
128*> 0 <= Q <= MIN(P,M-P,M-Q).
129*> \endverbatim
130*>
131*> \param[in,out] THETA
132*> \verbatim
133*> THETA is DOUBLE PRECISION array, dimension (Q)
134*> On entry, the angles THETA(1),...,THETA(Q) that, along with
135*> PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block
136*> form. On exit, the angles whose cosines and sines define the
137*> diagonal blocks in the CS decomposition.
138*> \endverbatim
139*>
140*> \param[in,out] PHI
141*> \verbatim
142*> PHI is DOUBLE PRECISION array, dimension (Q-1)
143*> The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),...,
144*> THETA(Q), define the matrix in bidiagonal-block form.
145*> \endverbatim
146*>
147*> \param[in,out] U1
148*> \verbatim
149*> U1 is COMPLEX*16 array, dimension (LDU1,P)
150*> On entry, a P-by-P matrix. On exit, U1 is postmultiplied
151*> by the left singular vector matrix common to [ B11 ; 0 ] and
152*> [ B12 0 0 ; 0 -I 0 0 ].
153*> \endverbatim
154*>
155*> \param[in] LDU1
156*> \verbatim
157*> LDU1 is INTEGER
158*> The leading dimension of the array U1, LDU1 >= MAX(1,P).
159*> \endverbatim
160*>
161*> \param[in,out] U2
162*> \verbatim
163*> U2 is COMPLEX*16 array, dimension (LDU2,M-P)
164*> On entry, an (M-P)-by-(M-P) matrix. On exit, U2 is
165*> postmultiplied by the left singular vector matrix common to
166*> [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ].
167*> \endverbatim
168*>
169*> \param[in] LDU2
170*> \verbatim
171*> LDU2 is INTEGER
172*> The leading dimension of the array U2, LDU2 >= MAX(1,M-P).
173*> \endverbatim
174*>
175*> \param[in,out] V1T
176*> \verbatim
177*> V1T is COMPLEX*16 array, dimension (LDV1T,Q)
178*> On entry, a Q-by-Q matrix. On exit, V1T is premultiplied
179*> by the conjugate transpose of the right singular vector
180*> matrix common to [ B11 ; 0 ] and [ B21 ; 0 ].
181*> \endverbatim
182*>
183*> \param[in] LDV1T
184*> \verbatim
185*> LDV1T is INTEGER
186*> The leading dimension of the array V1T, LDV1T >= MAX(1,Q).
187*> \endverbatim
188*>
189*> \param[in,out] V2T
190*> \verbatim
191*> V2T is COMPLEX*16 array, dimension (LDV2T,M-Q)
192*> On entry, an (M-Q)-by-(M-Q) matrix. On exit, V2T is
193*> premultiplied by the conjugate transpose of the right
194*> singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and
195*> [ B22 0 0 ; 0 0 I ].
196*> \endverbatim
197*>
198*> \param[in] LDV2T
199*> \verbatim
200*> LDV2T is INTEGER
201*> The leading dimension of the array V2T, LDV2T >= MAX(1,M-Q).
202*> \endverbatim
203*>
204*> \param[out] B11D
205*> \verbatim
206*> B11D is DOUBLE PRECISION array, dimension (Q)
207*> When ZBBCSD converges, B11D contains the cosines of THETA(1),
208*> ..., THETA(Q). If ZBBCSD fails to converge, then B11D
209*> contains the diagonal of the partially reduced top-left
210*> block.
211*> \endverbatim
212*>
213*> \param[out] B11E
214*> \verbatim
215*> B11E is DOUBLE PRECISION array, dimension (Q-1)
216*> When ZBBCSD converges, B11E contains zeros. If ZBBCSD fails
217*> to converge, then B11E contains the superdiagonal of the
218*> partially reduced top-left block.
219*> \endverbatim
220*>
221*> \param[out] B12D
222*> \verbatim
223*> B12D is DOUBLE PRECISION array, dimension (Q)
224*> When ZBBCSD converges, B12D contains the negative sines of
225*> THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
226*> B12D contains the diagonal of the partially reduced top-right
227*> block.
228*> \endverbatim
229*>
230*> \param[out] B12E
231*> \verbatim
232*> B12E is DOUBLE PRECISION array, dimension (Q-1)
233*> When ZBBCSD converges, B12E contains zeros. If ZBBCSD fails
234*> to converge, then B12E contains the subdiagonal of the
235*> partially reduced top-right block.
236*> \endverbatim
237*>
238*> \param[out] B21D
239*> \verbatim
240*> B21D is DOUBLE PRECISION array, dimension (Q)
241*> When ZBBCSD converges, B21D contains the negative sines of
242*> THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
243*> B21D contains the diagonal of the partially reduced bottom-left
244*> block.
245*> \endverbatim
246*>
247*> \param[out] B21E
248*> \verbatim
249*> B21E is DOUBLE PRECISION array, dimension (Q-1)
250*> When ZBBCSD converges, B21E contains zeros. If ZBBCSD fails
251*> to converge, then B21E contains the subdiagonal of the
252*> partially reduced bottom-left block.
253*> \endverbatim
254*>
255*> \param[out] B22D
256*> \verbatim
257*> B22D is DOUBLE PRECISION array, dimension (Q)
258*> When ZBBCSD converges, B22D contains the negative sines of
259*> THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
260*> B22D contains the diagonal of the partially reduced bottom-right
261*> block.
262*> \endverbatim
263*>
264*> \param[out] B22E
265*> \verbatim
266*> B22E is DOUBLE PRECISION array, dimension (Q-1)
267*> When ZBBCSD converges, B22E contains zeros. If ZBBCSD fails
268*> to converge, then B22E contains the subdiagonal of the
269*> partially reduced bottom-right block.
270*> \endverbatim
271*>
272*> \param[out] RWORK
273*> \verbatim
274*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
275*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
276*> \endverbatim
277*>
278*> \param[in] LRWORK
279*> \verbatim
280*> LRWORK is INTEGER
281*> The dimension of the array RWORK. LRWORK >= MAX(1,8*Q).
282*>
283*> If LRWORK = -1, then a workspace query is assumed; the
284*> routine only calculates the optimal size of the RWORK array,
285*> returns this value as the first entry of the work array, and
286*> no error message related to LRWORK is issued by XERBLA.
287*> \endverbatim
288*>
289*> \param[out] INFO
290*> \verbatim
291*> INFO is INTEGER
292*> = 0: successful exit.
293*> < 0: if INFO = -i, the i-th argument had an illegal value.
294*> > 0: if ZBBCSD did not converge, INFO specifies the number
295*> of nonzero entries in PHI, and B11D, B11E, etc.,
296*> contain the partially reduced matrix.
297*> \endverbatim
298*
299*> \par Internal Parameters:
300* =========================
301*>
302*> \verbatim
303*> TOLMUL DOUBLE PRECISION, default = MAX(10,MIN(100,EPS**(-1/8)))
304*> TOLMUL controls the convergence criterion of the QR loop.
305*> Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they
306*> are within TOLMUL*EPS of either bound.
307*> \endverbatim
308*
309*> \par References:
310* ================
311*>
312*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
313*> Algorithms, 50(1):33-65, 2009.
314*
315* Authors:
316* ========
317*
318*> \author Univ. of Tennessee
319*> \author Univ. of California Berkeley
320*> \author Univ. of Colorado Denver
321*> \author NAG Ltd.
322*
323*> \ingroup bbcsd
324*
325* =====================================================================
326 SUBROUTINE zbbcsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P,
327 $ Q,
328 $ THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T,
329 $ V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E,
330 $ B22D, B22E, RWORK, LRWORK, INFO )
331*
332* -- LAPACK computational routine --
333* -- LAPACK is a software package provided by Univ. of Tennessee, --
334* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
335*
336* .. Scalar Arguments ..
337 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
338 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q
339* ..
340* .. Array Arguments ..
341 DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
342 $ b21d( * ), b21e( * ), b22d( * ), b22e( * ),
343 $ phi( * ), theta( * ), rwork( * )
344 COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
345 $ v2t( ldv2t, * )
346* ..
347*
348* ===================================================================
349*
350* .. Parameters ..
351 INTEGER MAXITR
352 PARAMETER ( MAXITR = 6 )
353 double precision hundred, meighth, one, ten, zero
354 parameter( hundred = 100.0d0, meighth = -0.125d0,
355 $ one = 1.0d0, ten = 10.0d0, zero = 0.0d0 )
356 COMPLEX*16 NEGONECOMPLEX
357 parameter( negonecomplex = (-1.0d0,0.0d0) )
358 DOUBLE PRECISION PIOVER2
359 PARAMETER ( PIOVER2 = 1.57079632679489661923132169163975144210d0 )
360* ..
361* .. Local Scalars ..
362 LOGICAL COLMAJOR, LQUERY, RESTART11, RESTART12,
363 $ restart21, restart22, wantu1, wantu2, wantv1t,
364 $ wantv2t
365 INTEGER I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS,
366 $ iu2sn, iv1tcs, iv1tsn, iv2tcs, iv2tsn, j,
367 $ lrworkmin, lrworkopt, maxit, mini
368 DOUBLE PRECISION B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY,
369 $ eps, mu, nu, r, sigma11, sigma21,
370 $ temp, thetamax, thetamin, thresh, tol, tolmul,
371 $ unfl, x1, x2, y1, y2
372*
373 EXTERNAL dlartgp, dlartgs, dlas2, xerbla, zlasr,
374 $ zscal,
375 $ zswap
376* ..
377* .. External Functions ..
378 DOUBLE PRECISION DLAMCH
379 LOGICAL LSAME
380 EXTERNAL LSAME, DLAMCH
381* ..
382* .. Intrinsic Functions ..
383 INTRINSIC abs, atan2, cos, max, min, sin, sqrt
384* ..
385* .. Executable Statements ..
386*
387* Test input arguments
388*
389 info = 0
390 lquery = lrwork .EQ. -1
391 wantu1 = lsame( jobu1, 'Y' )
392 wantu2 = lsame( jobu2, 'Y' )
393 wantv1t = lsame( jobv1t, 'Y' )
394 wantv2t = lsame( jobv2t, 'Y' )
395 colmajor = .NOT. lsame( trans, 'T' )
396*
397 IF( m .LT. 0 ) THEN
398 info = -6
399 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
400 info = -7
401 ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
402 info = -8
403 ELSE IF( q .GT. p .OR. q .GT. m-p .OR. q .GT. m-q ) THEN
404 info = -8
405 ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
406 info = -12
407 ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
408 info = -14
409 ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
410 info = -16
411 ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
412 info = -18
413 END IF
414*
415* Quick return if Q = 0
416*
417 IF( info .EQ. 0 .AND. q .EQ. 0 ) THEN
418 lrworkmin = 1
419 rwork(1) = lrworkmin
420 RETURN
421 END IF
422*
423* Compute workspace
424*
425 IF( info .EQ. 0 ) THEN
426 iu1cs = 1
427 iu1sn = iu1cs + q
428 iu2cs = iu1sn + q
429 iu2sn = iu2cs + q
430 iv1tcs = iu2sn + q
431 iv1tsn = iv1tcs + q
432 iv2tcs = iv1tsn + q
433 iv2tsn = iv2tcs + q
434 lrworkopt = iv2tsn + q - 1
435 lrworkmin = lrworkopt
436 rwork(1) = lrworkopt
437 IF( lrwork .LT. lrworkmin .AND. .NOT. lquery ) THEN
438 info = -28
439 END IF
440 END IF
441*
442 IF( info .NE. 0 ) THEN
443 CALL xerbla( 'ZBBCSD', -info )
444 RETURN
445 ELSE IF( lquery ) THEN
446 RETURN
447 END IF
448*
449* Get machine constants
450*
451 eps = dlamch( 'Epsilon' )
452 unfl = dlamch( 'Safe minimum' )
453 tolmul = max( ten, min( hundred, eps**meighth ) )
454 tol = tolmul*eps
455 thresh = max( tol, maxitr*q*q*unfl )
456*
457* Test for negligible sines or cosines
458*
459 DO i = 1, q
460 IF( theta(i) .LT. thresh ) THEN
461 theta(i) = zero
462 ELSE IF( theta(i) .GT. piover2-thresh ) THEN
463 theta(i) = piover2
464 END IF
465 END DO
466 DO i = 1, q-1
467 IF( phi(i) .LT. thresh ) THEN
468 phi(i) = zero
469 ELSE IF( phi(i) .GT. piover2-thresh ) THEN
470 phi(i) = piover2
471 END IF
472 END DO
473*
474* Initial deflation
475*
476 imax = q
477 DO WHILE( imax .GT. 1 )
478 IF( phi(imax-1) .NE. zero ) THEN
479 EXIT
480 END IF
481 imax = imax - 1
482 END DO
483 imin = imax - 1
484 IF ( imin .GT. 1 ) THEN
485 DO WHILE( phi(imin-1) .NE. zero )
486 imin = imin - 1
487 IF ( imin .LE. 1 ) EXIT
488 END DO
489 END IF
490*
491* Initialize iteration counter
492*
493 maxit = maxitr*q*q
494 iter = 0
495*
496* Begin main iteration loop
497*
498 DO WHILE( imax .GT. 1 )
499*
500* Compute the matrix entries
501*
502 b11d(imin) = cos( theta(imin) )
503 b21d(imin) = -sin( theta(imin) )
504 DO i = imin, imax - 1
505 b11e(i) = -sin( theta(i) ) * sin( phi(i) )
506 b11d(i+1) = cos( theta(i+1) ) * cos( phi(i) )
507 b12d(i) = sin( theta(i) ) * cos( phi(i) )
508 b12e(i) = cos( theta(i+1) ) * sin( phi(i) )
509 b21e(i) = -cos( theta(i) ) * sin( phi(i) )
510 b21d(i+1) = -sin( theta(i+1) ) * cos( phi(i) )
511 b22d(i) = cos( theta(i) ) * cos( phi(i) )
512 b22e(i) = -sin( theta(i+1) ) * sin( phi(i) )
513 END DO
514 b12d(imax) = sin( theta(imax) )
515 b22d(imax) = cos( theta(imax) )
516*
517* Abort if not converging; otherwise, increment ITER
518*
519 IF( iter .GT. maxit ) THEN
520 info = 0
521 DO i = 1, q
522 IF( phi(i) .NE. zero )
523 $ info = info + 1
524 END DO
525 RETURN
526 END IF
527*
528 iter = iter + imax - imin
529*
530* Compute shifts
531*
532 thetamax = theta(imin)
533 thetamin = theta(imin)
534 DO i = imin+1, imax
535 IF( theta(i) > thetamax )
536 $ thetamax = theta(i)
537 IF( theta(i) < thetamin )
538 $ thetamin = theta(i)
539 END DO
540*
541 IF( thetamax .GT. piover2 - thresh ) THEN
542*
543* Zero on diagonals of B11 and B22; induce deflation with a
544* zero shift
545*
546 mu = zero
547 nu = one
548*
549 ELSE IF( thetamin .LT. thresh ) THEN
550*
551* Zero on diagonals of B12 and B22; induce deflation with a
552* zero shift
553*
554 mu = one
555 nu = zero
556*
557 ELSE
558*
559* Compute shifts for B11 and B21 and use the lesser
560*
561 CALL dlas2( b11d(imax-1), b11e(imax-1), b11d(imax),
562 $ sigma11,
563 $ dummy )
564 CALL dlas2( b21d(imax-1), b21e(imax-1), b21d(imax),
565 $ 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),
723 $ rwork(iv2tsn+i-1-1),
724 $ rwork(iv2tcs+i-1-1), r )
725 ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
726 CALL dlartgp( b22bulge, b22d(i-1),
727 $ rwork(iv2tsn+i-1-1),
728 $ rwork(iv2tcs+i-1-1), r )
729 ELSE IF( nu .LT. mu ) THEN
730 CALL dlartgs( b12e(i-1), b12d(i), nu,
731 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
732 ELSE
733 CALL dlartgs( b22e(i-1), b22d(i), mu,
734 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
735 END IF
736*
737 temp = rwork(iv1tcs+i-1)*b11d(i) + rwork(iv1tsn+i-1)*b11e(i)
738 b11e(i) = rwork(iv1tcs+i-1)*b11e(i) -
739 $ rwork(iv1tsn+i-1)*b11d(i)
740 b11d(i) = temp
741 b11bulge = rwork(iv1tsn+i-1)*b11d(i+1)
742 b11d(i+1) = rwork(iv1tcs+i-1)*b11d(i+1)
743 temp = rwork(iv1tcs+i-1)*b21d(i) + rwork(iv1tsn+i-1)*b21e(i)
744 b21e(i) = rwork(iv1tcs+i-1)*b21e(i) -
745 $ rwork(iv1tsn+i-1)*b21d(i)
746 b21d(i) = temp
747 b21bulge = rwork(iv1tsn+i-1)*b21d(i+1)
748 b21d(i+1) = rwork(iv1tcs+i-1)*b21d(i+1)
749 temp = rwork(iv2tcs+i-1-1)*b12e(i-1) +
750 $ rwork(iv2tsn+i-1-1)*b12d(i)
751 b12d(i) = rwork(iv2tcs+i-1-1)*b12d(i) -
752 $ rwork(iv2tsn+i-1-1)*b12e(i-1)
753 b12e(i-1) = temp
754 b12bulge = rwork(iv2tsn+i-1-1)*b12e(i)
755 b12e(i) = rwork(iv2tcs+i-1-1)*b12e(i)
756 temp = rwork(iv2tcs+i-1-1)*b22e(i-1) +
757 $ rwork(iv2tsn+i-1-1)*b22d(i)
758 b22d(i) = rwork(iv2tcs+i-1-1)*b22d(i) -
759 $ rwork(iv2tsn+i-1-1)*b22e(i-1)
760 b22e(i-1) = temp
761 b22bulge = rwork(iv2tsn+i-1-1)*b22e(i)
762 b22e(i) = rwork(iv2tcs+i-1-1)*b22e(i)
763*
764* Compute THETA(I)
765*
766 x1 = cos(phi(i-1))*b11d(i) + sin(phi(i-1))*b12e(i-1)
767 x2 = cos(phi(i-1))*b11bulge + sin(phi(i-1))*b12bulge
768 y1 = cos(phi(i-1))*b21d(i) + sin(phi(i-1))*b22e(i-1)
769 y2 = cos(phi(i-1))*b21bulge + sin(phi(i-1))*b22bulge
770*
771 theta(i) = atan2( sqrt(y1**2+y2**2), sqrt(x1**2+x2**2) )
772*
773* Determine if there are bulges to chase or if a new direct
774* summand has been reached
775*
776 restart11 = b11d(i)**2 + b11bulge**2 .LE. thresh**2
777 restart12 = b12e(i-1)**2 + b12bulge**2 .LE. thresh**2
778 restart21 = b21d(i)**2 + b21bulge**2 .LE. thresh**2
779 restart22 = b22e(i-1)**2 + b22bulge**2 .LE. thresh**2
780*
781* If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
782* B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
783* chasing by applying the original shift again.
784*
785 IF( .NOT. restart11 .AND. .NOT. restart12 ) THEN
786 CALL dlartgp( x2, x1, rwork(iu1sn+i-1),
787 $ rwork(iu1cs+i-1),
788 $ r )
789 ELSE IF( .NOT. restart11 .AND. restart12 ) THEN
790 CALL dlartgp( b11bulge, b11d(i), rwork(iu1sn+i-1),
791 $ rwork(iu1cs+i-1), r )
792 ELSE IF( restart11 .AND. .NOT. restart12 ) THEN
793 CALL dlartgp( b12bulge, b12e(i-1), rwork(iu1sn+i-1),
794 $ rwork(iu1cs+i-1), r )
795 ELSE IF( mu .LE. nu ) THEN
796 CALL dlartgs( b11e(i), b11d(i+1), mu,
797 $ rwork(iu1cs+i-1),
798 $ rwork(iu1sn+i-1) )
799 ELSE
800 CALL dlartgs( b12d(i), b12e(i), nu, rwork(iu1cs+i-1),
801 $ rwork(iu1sn+i-1) )
802 END IF
803 IF( .NOT. restart21 .AND. .NOT. restart22 ) THEN
804 CALL dlartgp( y2, y1, rwork(iu2sn+i-1),
805 $ rwork(iu2cs+i-1),
806 $ r )
807 ELSE IF( .NOT. restart21 .AND. restart22 ) THEN
808 CALL dlartgp( b21bulge, b21d(i), rwork(iu2sn+i-1),
809 $ rwork(iu2cs+i-1), r )
810 ELSE IF( restart21 .AND. .NOT. restart22 ) THEN
811 CALL dlartgp( b22bulge, b22e(i-1), rwork(iu2sn+i-1),
812 $ rwork(iu2cs+i-1), r )
813 ELSE IF( nu .LT. mu ) THEN
814 CALL dlartgs( b21e(i), b21e(i+1), nu,
815 $ rwork(iu2cs+i-1),
816 $ rwork(iu2sn+i-1) )
817 ELSE
818 CALL dlartgs( b22d(i), b22e(i), mu, rwork(iu2cs+i-1),
819 $ rwork(iu2sn+i-1) )
820 END IF
821 rwork(iu2cs+i-1) = -rwork(iu2cs+i-1)
822 rwork(iu2sn+i-1) = -rwork(iu2sn+i-1)
823*
824 temp = rwork(iu1cs+i-1)*b11e(i) + rwork(iu1sn+i-1)*b11d(i+1)
825 b11d(i+1) = rwork(iu1cs+i-1)*b11d(i+1) -
826 $ rwork(iu1sn+i-1)*b11e(i)
827 b11e(i) = temp
828 IF( i .LT. imax - 1 ) THEN
829 b11bulge = rwork(iu1sn+i-1)*b11e(i+1)
830 b11e(i+1) = rwork(iu1cs+i-1)*b11e(i+1)
831 END IF
832 temp = rwork(iu2cs+i-1)*b21e(i) + rwork(iu2sn+i-1)*b21d(i+1)
833 b21d(i+1) = rwork(iu2cs+i-1)*b21d(i+1) -
834 $ rwork(iu2sn+i-1)*b21e(i)
835 b21e(i) = temp
836 IF( i .LT. imax - 1 ) THEN
837 b21bulge = rwork(iu2sn+i-1)*b21e(i+1)
838 b21e(i+1) = rwork(iu2cs+i-1)*b21e(i+1)
839 END IF
840 temp = rwork(iu1cs+i-1)*b12d(i) + rwork(iu1sn+i-1)*b12e(i)
841 b12e(i) = rwork(iu1cs+i-1)*b12e(i) -
842 $ rwork(iu1sn+i-1)*b12d(i)
843 b12d(i) = temp
844 b12bulge = rwork(iu1sn+i-1)*b12d(i+1)
845 b12d(i+1) = rwork(iu1cs+i-1)*b12d(i+1)
846 temp = rwork(iu2cs+i-1)*b22d(i) + rwork(iu2sn+i-1)*b22e(i)
847 b22e(i) = rwork(iu2cs+i-1)*b22e(i) -
848 $ rwork(iu2sn+i-1)*b22d(i)
849 b22d(i) = temp
850 b22bulge = rwork(iu2sn+i-1)*b22d(i+1)
851 b22d(i+1) = rwork(iu2cs+i-1)*b22d(i+1)
852*
853 END DO
854*
855* Compute PHI(IMAX-1)
856*
857 x1 = sin(theta(imax-1))*b11e(imax-1) +
858 $ cos(theta(imax-1))*b21e(imax-1)
859 y1 = sin(theta(imax-1))*b12d(imax-1) +
860 $ cos(theta(imax-1))*b22d(imax-1)
861 y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
862*
863 phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
864*
865* Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
866*
867 restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
868 restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
869*
870 IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
871 CALL dlartgp( y2, y1, rwork(iv2tsn+imax-1-1),
872 $ rwork(iv2tcs+imax-1-1), r )
873 ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
874 CALL dlartgp( b12bulge, b12d(imax-1),
875 $ rwork(iv2tsn+imax-1-1),
876 $ rwork(iv2tcs+imax-1-1), r )
877 ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
878 CALL dlartgp( b22bulge, b22d(imax-1),
879 $ rwork(iv2tsn+imax-1-1),
880 $ rwork(iv2tcs+imax-1-1), r )
881 ELSE IF( nu .LT. mu ) THEN
882 CALL dlartgs( b12e(imax-1), b12d(imax), nu,
883 $ rwork(iv2tcs+imax-1-1),
884 $ rwork(iv2tsn+imax-1-1) )
885 ELSE
886 CALL dlartgs( b22e(imax-1), b22d(imax), mu,
887 $ rwork(iv2tcs+imax-1-1),
888 $ rwork(iv2tsn+imax-1-1) )
889 END IF
890*
891 temp = rwork(iv2tcs+imax-1-1)*b12e(imax-1) +
892 $ rwork(iv2tsn+imax-1-1)*b12d(imax)
893 b12d(imax) = rwork(iv2tcs+imax-1-1)*b12d(imax) -
894 $ rwork(iv2tsn+imax-1-1)*b12e(imax-1)
895 b12e(imax-1) = temp
896 temp = rwork(iv2tcs+imax-1-1)*b22e(imax-1) +
897 $ rwork(iv2tsn+imax-1-1)*b22d(imax)
898 b22d(imax) = rwork(iv2tcs+imax-1-1)*b22d(imax) -
899 $ rwork(iv2tsn+imax-1-1)*b22e(imax-1)
900 b22e(imax-1) = temp
901*
902* Update singular vectors
903*
904 IF( wantu1 ) THEN
905 IF( colmajor ) THEN
906 CALL zlasr( 'R', 'V', 'F', p, imax-imin+1,
907 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
908 $ u1(1,imin), ldu1 )
909 ELSE
910 CALL zlasr( 'L', 'V', 'F', imax-imin+1, p,
911 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
912 $ u1(imin,1), ldu1 )
913 END IF
914 END IF
915 IF( wantu2 ) THEN
916 IF( colmajor ) THEN
917 CALL zlasr( 'R', 'V', 'F', m-p, imax-imin+1,
918 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
919 $ u2(1,imin), ldu2 )
920 ELSE
921 CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-p,
922 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
923 $ u2(imin,1), ldu2 )
924 END IF
925 END IF
926 IF( wantv1t ) THEN
927 IF( colmajor ) THEN
928 CALL zlasr( 'L', 'V', 'F', imax-imin+1, q,
929 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
930 $ v1t(imin,1), ldv1t )
931 ELSE
932 CALL zlasr( 'R', 'V', 'F', q, imax-imin+1,
933 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
934 $ v1t(1,imin), ldv1t )
935 END IF
936 END IF
937 IF( wantv2t ) THEN
938 IF( colmajor ) THEN
939 CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-q,
940 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
941 $ v2t(imin,1), ldv2t )
942 ELSE
943 CALL zlasr( 'R', 'V', 'F', m-q, imax-imin+1,
944 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
945 $ v2t(1,imin), ldv2t )
946 END IF
947 END IF
948*
949* Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
950*
951 IF( b11e(imax-1)+b21e(imax-1) .GT. 0 ) THEN
952 b11d(imax) = -b11d(imax)
953 b21d(imax) = -b21d(imax)
954 IF( wantv1t ) THEN
955 IF( colmajor ) THEN
956 CALL zscal( q, negonecomplex, v1t(imax,1), ldv1t )
957 ELSE
958 CALL zscal( q, negonecomplex, v1t(1,imax), 1 )
959 END IF
960 END IF
961 END IF
962*
963* Compute THETA(IMAX)
964*
965 x1 = cos(phi(imax-1))*b11d(imax) +
966 $ sin(phi(imax-1))*b12e(imax-1)
967 y1 = cos(phi(imax-1))*b21d(imax) +
968 $ sin(phi(imax-1))*b22e(imax-1)
969*
970 theta(imax) = atan2( abs(y1), abs(x1) )
971*
972* Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
973* and B22(IMAX,IMAX-1)
974*
975 IF( b11d(imax)+b12e(imax-1) .LT. 0 ) THEN
976 b12d(imax) = -b12d(imax)
977 IF( wantu1 ) THEN
978 IF( colmajor ) THEN
979 CALL zscal( p, negonecomplex, u1(1,imax), 1 )
980 ELSE
981 CALL zscal( p, negonecomplex, u1(imax,1), ldu1 )
982 END IF
983 END IF
984 END IF
985 IF( b21d(imax)+b22e(imax-1) .GT. 0 ) THEN
986 b22d(imax) = -b22d(imax)
987 IF( wantu2 ) THEN
988 IF( colmajor ) THEN
989 CALL zscal( m-p, negonecomplex, u2(1,imax), 1 )
990 ELSE
991 CALL zscal( m-p, negonecomplex, u2(imax,1), ldu2 )
992 END IF
993 END IF
994 END IF
995*
996* Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
997*
998 IF( b12d(imax)+b22d(imax) .LT. 0 ) THEN
999 IF( wantv2t ) THEN
1000 IF( colmajor ) THEN
1001 CALL zscal( m-q, negonecomplex, v2t(imax,1),
1002 $ ldv2t )
1003 ELSE
1004 CALL zscal( m-q, negonecomplex, v2t(1,imax), 1 )
1005 END IF
1006 END IF
1007 END IF
1008*
1009* Test for negligible sines or cosines
1010*
1011 DO i = imin, imax
1012 IF( theta(i) .LT. thresh ) THEN
1013 theta(i) = zero
1014 ELSE IF( theta(i) .GT. piover2-thresh ) THEN
1015 theta(i) = piover2
1016 END IF
1017 END DO
1018 DO i = imin, imax-1
1019 IF( phi(i) .LT. thresh ) THEN
1020 phi(i) = zero
1021 ELSE IF( phi(i) .GT. piover2-thresh ) THEN
1022 phi(i) = piover2
1023 END IF
1024 END DO
1025*
1026* Deflate
1027*
1028 IF (imax .GT. 1) THEN
1029 DO WHILE( phi(imax-1) .EQ. zero )
1030 imax = imax - 1
1031 IF (imax .LE. 1) EXIT
1032 END DO
1033 END IF
1034 IF( imin .GT. imax - 1 )
1035 $ imin = imax - 1
1036 IF (imin .GT. 1) THEN
1037 DO WHILE (phi(imin-1) .NE. zero)
1038 imin = imin - 1
1039 IF (imin .LE. 1) EXIT
1040 END DO
1041 END IF
1042*
1043* Repeat main iteration loop
1044*
1045 END DO
1046*
1047* Postprocessing: order THETA from least to greatest
1048*
1049 DO i = 1, q
1050*
1051 mini = i
1052 thetamin = theta(i)
1053 DO j = i+1, q
1054 IF( theta(j) .LT. thetamin ) THEN
1055 mini = j
1056 thetamin = theta(j)
1057 END IF
1058 END DO
1059*
1060 IF( mini .NE. i ) THEN
1061 theta(mini) = theta(i)
1062 theta(i) = thetamin
1063 IF( colmajor ) THEN
1064 IF( wantu1 )
1065 $ CALL zswap( p, u1(1,i), 1, u1(1,mini), 1 )
1066 IF( wantu2 )
1067 $ CALL zswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1068 IF( wantv1t )
1069 $ CALL zswap( q, v1t(i,1), ldv1t, v1t(mini,1),
1070 $ ldv1t )
1071 IF( wantv2t )
1072 $ CALL zswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1073 $ ldv2t )
1074 ELSE
1075 IF( wantu1 )
1076 $ CALL zswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1077 IF( wantu2 )
1078 $ CALL zswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1079 IF( wantv1t )
1080 $ CALL zswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1081 IF( wantv2t )
1082 $ CALL zswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
1083 END IF
1084 END IF
1085*
1086 END DO
1087*
1088 RETURN
1089*
1090* End of ZBBCSD
1091*
1092 END
1093
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zbbcsd
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:331
dlartgp
subroutine dlartgp(f, g, cs, sn, r)
DLARTGP generates a plane rotation so that the diagonal is nonnegative.
Definition dlartgp.f:93
dlartgs
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:88
dlas2
subroutine dlas2(f, g, h, ssmin, ssmax)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition dlas2.f:103
zlasr
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:198
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81