OpenRadioss 2025.1.11
OpenRadioss project
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*> \htmlonly
9*> Download ZBBCSD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zbbcsd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zbbcsd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zbbcsd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
22* THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T,
23* V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E,
24* B22D, B22E, RWORK, LRWORK, INFO )
25*
26* .. Scalar Arguments ..
27* CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
28* INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q
29* ..
30* .. Array Arguments ..
31* DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
32* $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
33* $ PHI( * ), THETA( * ), RWORK( * )
34* COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
35* $ V2T( LDV2T, * )
36* ..
37*
38*
39*> \par Purpose:
40* =============
41*>
42*> \verbatim
43*>
44*> ZBBCSD computes the CS decomposition of a unitary matrix in
45*> bidiagonal-block form,
46*>
47*>
48*> [ B11 | B12 0 0 ]
49*> [ 0 | 0 -I 0 ]
50*> X = [----------------]
51*> [ B21 | B22 0 0 ]
52*> [ 0 | 0 0 I ]
53*>
54*> [ C | -S 0 0 ]
55*> [ U1 | ] [ 0 | 0 -I 0 ] [ V1 | ]**H
56*> = [---------] [---------------] [---------] .
57*> [ | U2 ] [ S | C 0 0 ] [ | V2 ]
58*> [ 0 | 0 0 I ]
59*>
60*> X is M-by-M, its top-left block is P-by-Q, and Q must be no larger
61*> than P, M-P, or M-Q. (If Q is not the smallest index, then X must be
62*> transposed and/or permuted. This can be done in constant time using
63*> the TRANS and SIGNS options. See ZUNCSD for details.)
64*>
65*> The bidiagonal matrices B11, B12, B21, and B22 are represented
66*> implicitly by angles THETA(1:Q) and PHI(1:Q-1).
67*>
68*> The unitary matrices U1, U2, V1T, and V2T are input/output.
69*> The input matrices are pre- or post-multiplied by the appropriate
70*> singular vector matrices.
71*> \endverbatim
72*
73* Arguments:
74* ==========
75*
76*> \param[in] JOBU1
77*> \verbatim
78*> JOBU1 is CHARACTER
79*> = 'Y': U1 is updated;
80*> otherwise: U1 is not updated.
81*> \endverbatim
82*>
83*> \param[in] JOBU2
84*> \verbatim
85*> JOBU2 is CHARACTER
86*> = 'Y': U2 is updated;
87*> otherwise: U2 is not updated.
88*> \endverbatim
89*>
90*> \param[in] JOBV1T
91*> \verbatim
92*> JOBV1T is CHARACTER
93*> = 'Y': V1T is updated;
94*> otherwise: V1T is not updated.
95*> \endverbatim
96*>
97*> \param[in] JOBV2T
98*> \verbatim
99*> JOBV2T is CHARACTER
100*> = 'Y': V2T is updated;
101*> otherwise: V2T is not updated.
102*> \endverbatim
103*>
104*> \param[in] TRANS
105*> \verbatim
106*> TRANS is CHARACTER
107*> = 'T': X, U1, U2, V1T, and V2T are stored in row-major
108*> order;
109*> otherwise: X, U1, U2, V1T, and V2T are stored in column-
110*> major order.
111*> \endverbatim
112*>
113*> \param[in] M
114*> \verbatim
115*> M is INTEGER
116*> The number of rows and columns in X, the unitary matrix in
117*> bidiagonal-block form.
118*> \endverbatim
119*>
120*> \param[in] P
121*> \verbatim
122*> P is INTEGER
123*> The number of rows in the top-left block of X. 0 <= P <= M.
124*> \endverbatim
125*>
126*> \param[in] Q
127*> \verbatim
128*> Q is INTEGER
129*> The number of columns in the top-left block of X.
130*> 0 <= Q <= MIN(P,M-P,M-Q).
131*> \endverbatim
132*>
133*> \param[in,out] THETA
134*> \verbatim
135*> THETA is DOUBLE PRECISION array, dimension (Q)
136*> On entry, the angles THETA(1),...,THETA(Q) that, along with
137*> PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block
138*> form. On exit, the angles whose cosines and sines define the
139*> diagonal blocks in the CS decomposition.
140*> \endverbatim
141*>
142*> \param[in,out] PHI
143*> \verbatim
144*> PHI is DOUBLE PRECISION array, dimension (Q-1)
145*> The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),...,
146*> THETA(Q), define the matrix in bidiagonal-block form.
147*> \endverbatim
148*>
149*> \param[in,out] U1
150*> \verbatim
151*> U1 is COMPLEX*16 array, dimension (LDU1,P)
152*> On entry, a P-by-P matrix. On exit, U1 is postmultiplied
153*> by the left singular vector matrix common to [ B11 ; 0 ] and
154*> [ B12 0 0 ; 0 -I 0 0 ].
155*> \endverbatim
156*>
157*> \param[in] LDU1
158*> \verbatim
159*> LDU1 is INTEGER
160*> The leading dimension of the array U1, LDU1 >= MAX(1,P).
161*> \endverbatim
162*>
163*> \param[in,out] U2
164*> \verbatim
165*> U2 is COMPLEX*16 array, dimension (LDU2,M-P)
166*> On entry, an (M-P)-by-(M-P) matrix. On exit, U2 is
167*> postmultiplied by the left singular vector matrix common to
168*> [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ].
169*> \endverbatim
170*>
171*> \param[in] LDU2
172*> \verbatim
173*> LDU2 is INTEGER
174*> The leading dimension of the array U2, LDU2 >= MAX(1,M-P).
175*> \endverbatim
176*>
177*> \param[in,out] V1T
178*> \verbatim
179*> V1T is COMPLEX*16 array, dimension (LDV1T,Q)
180*> On entry, a Q-by-Q matrix. On exit, V1T is premultiplied
181*> by the conjugate transpose of the right singular vector
182*> matrix common to [ B11 ; 0 ] and [ B21 ; 0 ].
183*> \endverbatim
184*>
185*> \param[in] LDV1T
186*> \verbatim
187*> LDV1T is INTEGER
188*> The leading dimension of the array V1T, LDV1T >= MAX(1,Q).
189*> \endverbatim
190*>
191*> \param[in,out] V2T
192*> \verbatim
193*> V2T is COMPLEX*16 array, dimension (LDV2T,M-Q)
194*> On entry, an (M-Q)-by-(M-Q) matrix. On exit, V2T is
195*> premultiplied by the conjugate transpose of the right
196*> singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and
197*> [ B22 0 0 ; 0 0 I ].
198*> \endverbatim
199*>
200*> \param[in] LDV2T
201*> \verbatim
202*> LDV2T is INTEGER
203*> The leading dimension of the array V2T, LDV2T >= MAX(1,M-Q).
204*> \endverbatim
205*>
206*> \param[out] B11D
207*> \verbatim
208*> B11D is DOUBLE PRECISION array, dimension (Q)
209*> When ZBBCSD converges, B11D contains the cosines of THETA(1),
210*> ..., THETA(Q). If ZBBCSD fails to converge, then B11D
211*> contains the diagonal of the partially reduced top-left
212*> block.
213*> \endverbatim
214*>
215*> \param[out] B11E
216*> \verbatim
217*> B11E is DOUBLE PRECISION array, dimension (Q-1)
218*> When ZBBCSD converges, B11E contains zeros. If ZBBCSD fails
219*> to converge, then B11E contains the superdiagonal of the
220*> partially reduced top-left block.
221*> \endverbatim
222*>
223*> \param[out] B12D
224*> \verbatim
225*> B12D is DOUBLE PRECISION array, dimension (Q)
226*> When ZBBCSD converges, B12D contains the negative sines of
227*> THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
228*> B12D contains the diagonal of the partially reduced top-right
229*> block.
230*> \endverbatim
231*>
232*> \param[out] B12E
233*> \verbatim
234*> B12E is DOUBLE PRECISION array, dimension (Q-1)
235*> When ZBBCSD converges, B12E contains zeros. If ZBBCSD fails
236*> to converge, then B12E contains the subdiagonal of the
237*> partially reduced top-right block.
238*> \endverbatim
239*>
240*> \param[out] B21D
241*> \verbatim
242*> B21D is DOUBLE PRECISION array, dimension (Q)
243*> When ZBBCSD converges, B21D contains the negative sines of
244*> THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
245*> B21D contains the diagonal of the partially reduced bottom-left
246*> block.
247*> \endverbatim
248*>
249*> \param[out] B21E
250*> \verbatim
251*> B21E is DOUBLE PRECISION array, dimension (Q-1)
252*> When ZBBCSD converges, B21E contains zeros. If ZBBCSD fails
253*> to converge, then B21E contains the subdiagonal of the
254*> partially reduced bottom-left block.
255*> \endverbatim
256*>
257*> \param[out] B22D
258*> \verbatim
259*> B22D is DOUBLE PRECISION array, dimension (Q)
260*> When ZBBCSD converges, B22D contains the negative sines of
261*> THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
262*> B22D contains the diagonal of the partially reduced bottom-right
263*> block.
264*> \endverbatim
265*>
266*> \param[out] B22E
267*> \verbatim
268*> B22E is DOUBLE PRECISION array, dimension (Q-1)
269*> When ZBBCSD converges, B22E contains zeros. If ZBBCSD fails
270*> to converge, then B22E contains the subdiagonal of the
271*> partially reduced bottom-right block.
272*> \endverbatim
273*>
274*> \param[out] RWORK
275*> \verbatim
276*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
277*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
278*> \endverbatim
279*>
280*> \param[in] LRWORK
281*> \verbatim
282*> LRWORK is INTEGER
283*> The dimension of the array RWORK. LRWORK >= MAX(1,8*Q).
284*>
285*> If LRWORK = -1, then a workspace query is assumed; the
286*> routine only calculates the optimal size of the RWORK array,
287*> returns this value as the first entry of the work array, and
288*> no error message related to LRWORK is issued by XERBLA.
289*> \endverbatim
290*>
291*> \param[out] INFO
292*> \verbatim
293*> INFO is INTEGER
294*> = 0: successful exit.
295*> < 0: if INFO = -i, the i-th argument had an illegal value.
296*> > 0: if ZBBCSD did not converge, INFO specifies the number
297*> of nonzero entries in PHI, and B11D, B11E, etc.,
298*> contain the partially reduced matrix.
299*> \endverbatim
300*
301*> \par Internal Parameters:
302* =========================
303*>
304*> \verbatim
305*> TOLMUL DOUBLE PRECISION, default = MAX(10,MIN(100,EPS**(-1/8)))
306*> TOLMUL controls the convergence criterion of the QR loop.
307*> Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they
308*> are within TOLMUL*EPS of either bound.
309*> \endverbatim
310*
311*> \par References:
312* ================
313*>
314*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
315*> Algorithms, 50(1):33-65, 2009.
316*
317* Authors:
318* ========
319*
320*> \author Univ. of Tennessee
321*> \author Univ. of California Berkeley
322*> \author Univ. of Colorado Denver
323*> \author NAG Ltd.
324*
325*> \ingroup complex16OTHERcomputational
326*
327* =====================================================================
328 SUBROUTINE zbbcsd( 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 DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
343 $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
344 $ PHI( * ), THETA( * ), RWORK( * )
345 COMPLEX*16 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 COMPLEX*16 NEGONECOMPLEX
358 parameter( negonecomplex = (-1.0d0,0.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 $ LRWORKMIN, LRWORKOPT, 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 dlartgp, dlartgs, dlas2, xerbla, zlasr, 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), sigma11,
562 $ dummy )
563 CALL dlas2( b21d(imax-1), b21e(imax-1), b21d(imax), sigma21,
564 $ dummy )
565*
566 IF( sigma11 .LE. sigma21 ) THEN
567 mu = sigma11
568 nu = sqrt( one - mu**2 )
569 IF( mu .LT. thresh ) THEN
570 mu = zero
571 nu = one
572 END IF
573 ELSE
574 nu = sigma21
575 mu = sqrt( 1.0 - nu**2 )
576 IF( nu .LT. thresh ) THEN
577 mu = one
578 nu = zero
579 END IF
580 END IF
581 END IF
582*
583* Rotate to produce bulges in B11 and B21
584*
585 IF( mu .LE. nu ) THEN
586 CALL dlartgs( b11d(imin), b11e(imin), mu,
587 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
588 ELSE
589 CALL dlartgs( b21d(imin), b21e(imin), nu,
590 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
591 END IF
592*
593 temp = rwork(iv1tcs+imin-1)*b11d(imin) +
594 $ rwork(iv1tsn+imin-1)*b11e(imin)
595 b11e(imin) = rwork(iv1tcs+imin-1)*b11e(imin) -
596 $ rwork(iv1tsn+imin-1)*b11d(imin)
597 b11d(imin) = temp
598 b11bulge = rwork(iv1tsn+imin-1)*b11d(imin+1)
599 b11d(imin+1) = rwork(iv1tcs+imin-1)*b11d(imin+1)
600 temp = rwork(iv1tcs+imin-1)*b21d(imin) +
601 $ rwork(iv1tsn+imin-1)*b21e(imin)
602 b21e(imin) = rwork(iv1tcs+imin-1)*b21e(imin) -
603 $ rwork(iv1tsn+imin-1)*b21d(imin)
604 b21d(imin) = temp
605 b21bulge = rwork(iv1tsn+imin-1)*b21d(imin+1)
606 b21d(imin+1) = rwork(iv1tcs+imin-1)*b21d(imin+1)
607*
608* Compute THETA(IMIN)
609*
610 theta( imin ) = atan2( sqrt( b21d(imin)**2+b21bulge**2 ),
611 $ sqrt( b11d(imin)**2+b11bulge**2 ) )
612*
613* Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
614*
615 IF( b11d(imin)**2+b11bulge**2 .GT. thresh**2 ) THEN
616 CALL dlartgp( b11bulge, b11d(imin), rwork(iu1sn+imin-1),
617 $ rwork(iu1cs+imin-1), r )
618 ELSE IF( mu .LE. nu ) THEN
619 CALL dlartgs( b11e( imin ), b11d( imin + 1 ), mu,
620 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
621 ELSE
622 CALL dlartgs( b12d( imin ), b12e( imin ), nu,
623 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
624 END IF
625 IF( b21d(imin)**2+b21bulge**2 .GT. thresh**2 ) THEN
626 CALL dlartgp( b21bulge, b21d(imin), rwork(iu2sn+imin-1),
627 $ rwork(iu2cs+imin-1), r )
628 ELSE IF( nu .LT. mu ) THEN
629 CALL dlartgs( b21e( imin ), b21d( imin + 1 ), nu,
630 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
631 ELSE
632 CALL dlartgs( b22d(imin), b22e(imin), mu,
633 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
634 END IF
635 rwork(iu2cs+imin-1) = -rwork(iu2cs+imin-1)
636 rwork(iu2sn+imin-1) = -rwork(iu2sn+imin-1)
637*
638 temp = rwork(iu1cs+imin-1)*b11e(imin) +
639 $ rwork(iu1sn+imin-1)*b11d(imin+1)
640 b11d(imin+1) = rwork(iu1cs+imin-1)*b11d(imin+1) -
641 $ rwork(iu1sn+imin-1)*b11e(imin)
642 b11e(imin) = temp
643 IF( imax .GT. imin+1 ) THEN
644 b11bulge = rwork(iu1sn+imin-1)*b11e(imin+1)
645 b11e(imin+1) = rwork(iu1cs+imin-1)*b11e(imin+1)
646 END IF
647 temp = rwork(iu1cs+imin-1)*b12d(imin) +
648 $ rwork(iu1sn+imin-1)*b12e(imin)
649 b12e(imin) = rwork(iu1cs+imin-1)*b12e(imin) -
650 $ rwork(iu1sn+imin-1)*b12d(imin)
651 b12d(imin) = temp
652 b12bulge = rwork(iu1sn+imin-1)*b12d(imin+1)
653 b12d(imin+1) = rwork(iu1cs+imin-1)*b12d(imin+1)
654 temp = rwork(iu2cs+imin-1)*b21e(imin) +
655 $ rwork(iu2sn+imin-1)*b21d(imin+1)
656 b21d(imin+1) = rwork(iu2cs+imin-1)*b21d(imin+1) -
657 $ rwork(iu2sn+imin-1)*b21e(imin)
658 b21e(imin) = temp
659 IF( imax .GT. imin+1 ) THEN
660 b21bulge = rwork(iu2sn+imin-1)*b21e(imin+1)
661 b21e(imin+1) = rwork(iu2cs+imin-1)*b21e(imin+1)
662 END IF
663 temp = rwork(iu2cs+imin-1)*b22d(imin) +
664 $ rwork(iu2sn+imin-1)*b22e(imin)
665 b22e(imin) = rwork(iu2cs+imin-1)*b22e(imin) -
666 $ rwork(iu2sn+imin-1)*b22d(imin)
667 b22d(imin) = temp
668 b22bulge = rwork(iu2sn+imin-1)*b22d(imin+1)
669 b22d(imin+1) = rwork(iu2cs+imin-1)*b22d(imin+1)
670*
671* Inner loop: chase bulges from B11(IMIN,IMIN+2),
672* B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
673* bottom-right
674*
675 DO i = imin+1, imax-1
676*
677* Compute PHI(I-1)
678*
679 x1 = sin(theta(i-1))*b11e(i-1) + cos(theta(i-1))*b21e(i-1)
680 x2 = sin(theta(i-1))*b11bulge + cos(theta(i-1))*b21bulge
681 y1 = sin(theta(i-1))*b12d(i-1) + cos(theta(i-1))*b22d(i-1)
682 y2 = sin(theta(i-1))*b12bulge + cos(theta(i-1))*b22bulge
683*
684 phi(i-1) = atan2( sqrt(x1**2+x2**2), sqrt(y1**2+y2**2) )
685*
686* Determine if there are bulges to chase or if a new direct
687* summand has been reached
688*
689 restart11 = b11e(i-1)**2 + b11bulge**2 .LE. thresh**2
690 restart21 = b21e(i-1)**2 + b21bulge**2 .LE. thresh**2
691 restart12 = b12d(i-1)**2 + b12bulge**2 .LE. thresh**2
692 restart22 = b22d(i-1)**2 + b22bulge**2 .LE. thresh**2
693*
694* If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
695* B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
696* chasing by applying the original shift again.
697*
698 IF( .NOT. restart11 .AND. .NOT. restart21 ) THEN
699 CALL dlartgp( x2, x1, rwork(iv1tsn+i-1),
700 $ rwork(iv1tcs+i-1), r )
701 ELSE IF( .NOT. restart11 .AND. restart21 ) THEN
702 CALL dlartgp( b11bulge, b11e(i-1), rwork(iv1tsn+i-1),
703 $ rwork(iv1tcs+i-1), r )
704 ELSE IF( restart11 .AND. .NOT. restart21 ) THEN
705 CALL dlartgp( b21bulge, b21e(i-1), rwork(iv1tsn+i-1),
706 $ rwork(iv1tcs+i-1), r )
707 ELSE IF( mu .LE. nu ) THEN
708 CALL dlartgs( b11d(i), b11e(i), mu, rwork(iv1tcs+i-1),
709 $ rwork(iv1tsn+i-1) )
710 ELSE
711 CALL dlartgs( b21d(i), b21e(i), nu, rwork(iv1tcs+i-1),
712 $ rwork(iv1tsn+i-1) )
713 END IF
714 rwork(iv1tcs+i-1) = -rwork(iv1tcs+i-1)
715 rwork(iv1tsn+i-1) = -rwork(iv1tsn+i-1)
716 IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
717 CALL dlartgp( y2, y1, rwork(iv2tsn+i-1-1),
718 $ rwork(iv2tcs+i-1-1), r )
719 ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
720 CALL dlartgp( b12bulge, b12d(i-1), rwork(iv2tsn+i-1-1),
721 $ rwork(iv2tcs+i-1-1), r )
722 ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
723 CALL dlartgp( b22bulge, b22d(i-1), rwork(iv2tsn+i-1-1),
724 $ rwork(iv2tcs+i-1-1), r )
725 ELSE IF( nu .LT. mu ) THEN
726 CALL dlartgs( b12e(i-1), b12d(i), nu,
727 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
728 ELSE
729 CALL dlartgs( b22e(i-1), b22d(i), mu,
730 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
731 END IF
732*
733 temp = rwork(iv1tcs+i-1)*b11d(i) + rwork(iv1tsn+i-1)*b11e(i)
734 b11e(i) = rwork(iv1tcs+i-1)*b11e(i) -
735 $ rwork(iv1tsn+i-1)*b11d(i)
736 b11d(i) = temp
737 b11bulge = rwork(iv1tsn+i-1)*b11d(i+1)
738 b11d(i+1) = rwork(iv1tcs+i-1)*b11d(i+1)
739 temp = rwork(iv1tcs+i-1)*b21d(i) + rwork(iv1tsn+i-1)*b21e(i)
740 b21e(i) = rwork(iv1tcs+i-1)*b21e(i) -
741 $ rwork(iv1tsn+i-1)*b21d(i)
742 b21d(i) = temp
743 b21bulge = rwork(iv1tsn+i-1)*b21d(i+1)
744 b21d(i+1) = rwork(iv1tcs+i-1)*b21d(i+1)
745 temp = rwork(iv2tcs+i-1-1)*b12e(i-1) +
746 $ rwork(iv2tsn+i-1-1)*b12d(i)
747 b12d(i) = rwork(iv2tcs+i-1-1)*b12d(i) -
748 $ rwork(iv2tsn+i-1-1)*b12e(i-1)
749 b12e(i-1) = temp
750 b12bulge = rwork(iv2tsn+i-1-1)*b12e(i)
751 b12e(i) = rwork(iv2tcs+i-1-1)*b12e(i)
752 temp = rwork(iv2tcs+i-1-1)*b22e(i-1) +
753 $ rwork(iv2tsn+i-1-1)*b22d(i)
754 b22d(i) = rwork(iv2tcs+i-1-1)*b22d(i) -
755 $ rwork(iv2tsn+i-1-1)*b22e(i-1)
756 b22e(i-1) = temp
757 b22bulge = rwork(iv2tsn+i-1-1)*b22e(i)
758 b22e(i) = rwork(iv2tcs+i-1-1)*b22e(i)
759*
760* Compute THETA(I)
761*
762 x1 = cos(phi(i-1))*b11d(i) + sin(phi(i-1))*b12e(i-1)
763 x2 = cos(phi(i-1))*b11bulge + sin(phi(i-1))*b12bulge
764 y1 = cos(phi(i-1))*b21d(i) + sin(phi(i-1))*b22e(i-1)
765 y2 = cos(phi(i-1))*b21bulge + sin(phi(i-1))*b22bulge
766*
767 theta(i) = atan2( sqrt(y1**2+y2**2), sqrt(x1**2+x2**2) )
768*
769* Determine if there are bulges to chase or if a new direct
770* summand has been reached
771*
772 restart11 = b11d(i)**2 + b11bulge**2 .LE. thresh**2
773 restart12 = b12e(i-1)**2 + b12bulge**2 .LE. thresh**2
774 restart21 = b21d(i)**2 + b21bulge**2 .LE. thresh**2
775 restart22 = b22e(i-1)**2 + b22bulge**2 .LE. thresh**2
776*
777* If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
778* B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
779* chasing by applying the original shift again.
780*
781 IF( .NOT. restart11 .AND. .NOT. restart12 ) THEN
782 CALL dlartgp( x2, x1, rwork(iu1sn+i-1), rwork(iu1cs+i-1),
783 $ r )
784 ELSE IF( .NOT. restart11 .AND. restart12 ) THEN
785 CALL dlartgp( b11bulge, b11d(i), rwork(iu1sn+i-1),
786 $ rwork(iu1cs+i-1), r )
787 ELSE IF( restart11 .AND. .NOT. restart12 ) THEN
788 CALL dlartgp( b12bulge, b12e(i-1), rwork(iu1sn+i-1),
789 $ rwork(iu1cs+i-1), r )
790 ELSE IF( mu .LE. nu ) THEN
791 CALL dlartgs( b11e(i), b11d(i+1), mu, rwork(iu1cs+i-1),
792 $ rwork(iu1sn+i-1) )
793 ELSE
794 CALL dlartgs( b12d(i), b12e(i), nu, rwork(iu1cs+i-1),
795 $ rwork(iu1sn+i-1) )
796 END IF
797 IF( .NOT. restart21 .AND. .NOT. restart22 ) THEN
798 CALL dlartgp( y2, y1, rwork(iu2sn+i-1), rwork(iu2cs+i-1),
799 $ r )
800 ELSE IF( .NOT. restart21 .AND. restart22 ) THEN
801 CALL dlartgp( b21bulge, b21d(i), rwork(iu2sn+i-1),
802 $ rwork(iu2cs+i-1), r )
803 ELSE IF( restart21 .AND. .NOT. restart22 ) THEN
804 CALL dlartgp( b22bulge, b22e(i-1), rwork(iu2sn+i-1),
805 $ rwork(iu2cs+i-1), r )
806 ELSE IF( nu .LT. mu ) THEN
807 CALL dlartgs( b21e(i), b21e(i+1), nu, rwork(iu2cs+i-1),
808 $ rwork(iu2sn+i-1) )
809 ELSE
810 CALL dlartgs( b22d(i), b22e(i), mu, rwork(iu2cs+i-1),
811 $ rwork(iu2sn+i-1) )
812 END IF
813 rwork(iu2cs+i-1) = -rwork(iu2cs+i-1)
814 rwork(iu2sn+i-1) = -rwork(iu2sn+i-1)
815*
816 temp = rwork(iu1cs+i-1)*b11e(i) + rwork(iu1sn+i-1)*b11d(i+1)
817 b11d(i+1) = rwork(iu1cs+i-1)*b11d(i+1) -
818 $ rwork(iu1sn+i-1)*b11e(i)
819 b11e(i) = temp
820 IF( i .LT. imax - 1 ) THEN
821 b11bulge = rwork(iu1sn+i-1)*b11e(i+1)
822 b11e(i+1) = rwork(iu1cs+i-1)*b11e(i+1)
823 END IF
824 temp = rwork(iu2cs+i-1)*b21e(i) + rwork(iu2sn+i-1)*b21d(i+1)
825 b21d(i+1) = rwork(iu2cs+i-1)*b21d(i+1) -
826 $ rwork(iu2sn+i-1)*b21e(i)
827 b21e(i) = temp
828 IF( i .LT. imax - 1 ) THEN
829 b21bulge = rwork(iu2sn+i-1)*b21e(i+1)
830 b21e(i+1) = rwork(iu2cs+i-1)*b21e(i+1)
831 END IF
832 temp = rwork(iu1cs+i-1)*b12d(i) + rwork(iu1sn+i-1)*b12e(i)
833 b12e(i) = rwork(iu1cs+i-1)*b12e(i) -
834 $ rwork(iu1sn+i-1)*b12d(i)
835 b12d(i) = temp
836 b12bulge = rwork(iu1sn+i-1)*b12d(i+1)
837 b12d(i+1) = rwork(iu1cs+i-1)*b12d(i+1)
838 temp = rwork(iu2cs+i-1)*b22d(i) + rwork(iu2sn+i-1)*b22e(i)
839 b22e(i) = rwork(iu2cs+i-1)*b22e(i) -
840 $ rwork(iu2sn+i-1)*b22d(i)
841 b22d(i) = temp
842 b22bulge = rwork(iu2sn+i-1)*b22d(i+1)
843 b22d(i+1) = rwork(iu2cs+i-1)*b22d(i+1)
844*
845 END DO
846*
847* Compute PHI(IMAX-1)
848*
849 x1 = sin(theta(imax-1))*b11e(imax-1) +
850 $ cos(theta(imax-1))*b21e(imax-1)
851 y1 = sin(theta(imax-1))*b12d(imax-1) +
852 $ cos(theta(imax-1))*b22d(imax-1)
853 y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
854*
855 phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
856*
857* Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
858*
859 restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
860 restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
861*
862 IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
863 CALL dlartgp( y2, y1, rwork(iv2tsn+imax-1-1),
864 $ rwork(iv2tcs+imax-1-1), r )
865 ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
866 CALL dlartgp( b12bulge, b12d(imax-1),
867 $ rwork(iv2tsn+imax-1-1),
868 $ rwork(iv2tcs+imax-1-1), r )
869 ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
870 CALL dlartgp( b22bulge, b22d(imax-1),
871 $ rwork(iv2tsn+imax-1-1),
872 $ rwork(iv2tcs+imax-1-1), r )
873 ELSE IF( nu .LT. mu ) THEN
874 CALL dlartgs( b12e(imax-1), b12d(imax), nu,
875 $ rwork(iv2tcs+imax-1-1),
876 $ rwork(iv2tsn+imax-1-1) )
877 ELSE
878 CALL dlartgs( b22e(imax-1), b22d(imax), mu,
879 $ rwork(iv2tcs+imax-1-1),
880 $ rwork(iv2tsn+imax-1-1) )
881 END IF
882*
883 temp = rwork(iv2tcs+imax-1-1)*b12e(imax-1) +
884 $ rwork(iv2tsn+imax-1-1)*b12d(imax)
885 b12d(imax) = rwork(iv2tcs+imax-1-1)*b12d(imax) -
886 $ rwork(iv2tsn+imax-1-1)*b12e(imax-1)
887 b12e(imax-1) = temp
888 temp = rwork(iv2tcs+imax-1-1)*b22e(imax-1) +
889 $ rwork(iv2tsn+imax-1-1)*b22d(imax)
890 b22d(imax) = rwork(iv2tcs+imax-1-1)*b22d(imax) -
891 $ rwork(iv2tsn+imax-1-1)*b22e(imax-1)
892 b22e(imax-1) = temp
893*
894* Update singular vectors
895*
896 IF( wantu1 ) THEN
897 IF( colmajor ) THEN
898 CALL zlasr( 'R', 'V', 'F', p, imax-imin+1,
899 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
900 $ u1(1,imin), ldu1 )
901 ELSE
902 CALL zlasr( 'L', 'V', 'F', imax-imin+1, p,
903 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
904 $ u1(imin,1), ldu1 )
905 END IF
906 END IF
907 IF( wantu2 ) THEN
908 IF( colmajor ) THEN
909 CALL zlasr( 'R', 'V', 'F', m-p, imax-imin+1,
910 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
911 $ u2(1,imin), ldu2 )
912 ELSE
913 CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-p,
914 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
915 $ u2(imin,1), ldu2 )
916 END IF
917 END IF
918 IF( wantv1t ) THEN
919 IF( colmajor ) THEN
920 CALL zlasr( 'L', 'V', 'f', IMAX-IMIN+1, Q,
921 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1),
922 $ V1T(IMIN,1), LDV1T )
923 ELSE
924 CALL ZLASR( 'r', 'v', 'f', Q, IMAX-IMIN+1,
925 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1),
926 $ V1T(1,IMIN), LDV1T )
927 END IF
928 END IF
929 IF( WANTV2T ) THEN
930 IF( COLMAJOR ) THEN
931 CALL ZLASR( 'l', 'v', 'f', IMAX-IMIN+1, M-Q,
932 $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1),
933 $ V2T(IMIN,1), LDV2T )
934 ELSE
935 CALL ZLASR( 'r', 'v', 'f', M-Q, IMAX-IMIN+1,
936 $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1),
937 $ V2T(1,IMIN), LDV2T )
938 END IF
939 END IF
940*
941* Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
942*
943.GT. IF( B11E(IMAX-1)+B21E(IMAX-1) 0 ) THEN
944 B11D(IMAX) = -B11D(IMAX)
945 B21D(IMAX) = -B21D(IMAX)
946 IF( WANTV1T ) THEN
947 IF( COLMAJOR ) THEN
948 CALL ZSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T )
949 ELSE
950 CALL ZSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 )
951 END IF
952 END IF
953 END IF
954*
955* Compute THETA(IMAX)
956*
957 X1 = COS(PHI(IMAX-1))*B11D(IMAX) +
958 $ SIN(PHI(IMAX-1))*B12E(IMAX-1)
959 Y1 = COS(PHI(IMAX-1))*B21D(IMAX) +
960 $ SIN(PHI(IMAX-1))*B22E(IMAX-1)
961*
962 THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) )
963*
964* Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
965* and B22(IMAX,IMAX-1)
966*
967.LT. IF( B11D(IMAX)+B12E(IMAX-1) 0 ) THEN
968 B12D(IMAX) = -B12D(IMAX)
969 IF( WANTU1 ) THEN
970 IF( COLMAJOR ) THEN
971 CALL ZSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 )
972 ELSE
973 CALL ZSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 )
974 END IF
975 END IF
976 END IF
977.GT. IF( B21D(IMAX)+B22E(IMAX-1) 0 ) THEN
978 B22D(IMAX) = -B22D(IMAX)
979 IF( WANTU2 ) THEN
980 IF( COLMAJOR ) THEN
981 CALL ZSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 )
982 ELSE
983 CALL ZSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 )
984 END IF
985 END IF
986 END IF
987*
988* Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
989*
990.LT. IF( B12D(IMAX)+B22D(IMAX) 0 ) THEN
991 IF( WANTV2T ) THEN
992 IF( COLMAJOR ) THEN
993 CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T )
994 ELSE
995 CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 )
996 END IF
997 END IF
998 END IF
999*
1000* Test for negligible sines or cosines
1001*
1002 DO I = IMIN, IMAX
1003.LT. IF( THETA(I) THRESH ) THEN
1004 THETA(I) = ZERO
1005.GT. ELSE IF( THETA(I) PIOVER2-THRESH ) THEN
1006 THETA(I) = PIOVER2
1007 END IF
1008 END DO
1009 DO I = IMIN, IMAX-1
1010.LT. IF( PHI(I) THRESH ) THEN
1011 PHI(I) = ZERO
1012.GT. ELSE IF( PHI(I) PIOVER2-THRESH ) THEN
1013 PHI(I) = PIOVER2
1014 END IF
1015 END DO
1016*
1017* Deflate
1018*
1019.GT. IF (IMAX 1) THEN
1020.EQ. DO WHILE( PHI(IMAX-1) ZERO )
1021 IMAX = IMAX - 1
1022.LE. IF (IMAX 1) EXIT
1023 END DO
1024 END IF
1025.GT. IF( IMIN IMAX - 1 )
1026 $ IMIN = IMAX - 1
1027.GT. IF (IMIN 1) THEN
1028.NE. DO WHILE (PHI(IMIN-1) ZERO)
1029 IMIN = IMIN - 1
1030.LE. IF (IMIN 1) EXIT
1031 END DO
1032 END IF
1033*
1034* Repeat main iteration loop
1035*
1036 END DO
1037*
1038* Postprocessing: order THETA from least to greatest
1039*
1040 DO I = 1, Q
1041*
1042 MINI = I
1043 THETAMIN = THETA(I)
1044 DO J = I+1, Q
1045.LT. IF( THETA(J) THETAMIN ) THEN
1046 MINI = J
1047 THETAMIN = THETA(J)
1048 END IF
1049 END DO
1050*
1051.NE. IF( MINI I ) THEN
1052 THETA(MINI) = THETA(I)
1053 THETA(I) = THETAMIN
1054 IF( COLMAJOR ) THEN
1055 IF( WANTU1 )
1056 $ CALL ZSWAP( P, U1(1,I), 1, U1(1,MINI), 1 )
1057 IF( WANTU2 )
1058 $ CALL ZSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 )
1059 IF( WANTV1T )
1060 $ CALL ZSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T )
1061 IF( WANTV2T )
1062 $ CALL ZSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1),
1063 $ LDV2T )
1064 ELSE
1065 IF( WANTU1 )
1066 $ CALL ZSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 )
1067 IF( WANTU2 )
1068 $ CALL ZSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 )
1069 IF( WANTV1T )
1070 $ CALL ZSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 )
1071 IF( WANTV2T )
1072 $ CALL ZSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 )
1073 END IF
1074 END IF
1075*
1076 END DO
1077*
1078 RETURN
1079*
1080* End of ZBBCSD
1081*
1082 END
1083
dlas2
subroutine dlas2(f, g, h, ssmin, ssmax)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition dlas2.f:107
dlartgp
subroutine dlartgp(f, g, cs, sn, r)
DLARTGP generates a plane rotation so that the diagonal is nonnegative.
Definition dlartgp.f:95
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:90
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
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:200
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:332
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21