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