LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ sbdsqr()

 subroutine sbdsqr ( character uplo, integer n, integer ncvt, integer nru, integer ncc, real, dimension( * ) d, real, dimension( * ) e, real, dimension( ldvt, * ) vt, integer ldvt, real, dimension( ldu, * ) u, integer ldu, real, dimension( ldc, * ) c, integer ldc, real, dimension( * ) work, integer info )

SBDSQR

Purpose:
``` SBDSQR computes the singular values and, optionally, the right and/or
left singular vectors from the singular value decomposition (SVD) of
a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
zero-shift QR algorithm.  The SVD of B has the form

B = Q * S * P**T

where S is the diagonal matrix of singular values, Q is an orthogonal
matrix of left singular vectors, and P is an orthogonal matrix of
right singular vectors.  If left singular vectors are requested, this
subroutine actually returns U*Q instead of Q, and, if right singular
vectors are requested, this subroutine returns P**T*VT instead of
P**T, for given real input matrices U and VT.  When U and VT are the
orthogonal matrices that reduce a general matrix A to bidiagonal
form:  A = U*B*VT, as computed by SGEBRD, then

A = (U*Q) * S * (P**T*VT)

is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
for a given real input matrix C.

See "Computing  Small Singular Values of Bidiagonal Matrices With
Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
no. 5, pp. 873-912, Sept 1990) and
"Accurate singular values and differential qd algorithms," by
B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
Department, University of California at Berkeley, July 1992
for a detailed description of the algorithm.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 = 'U': B is upper bidiagonal; = 'L': B is lower bidiagonal.``` [in] N ``` N is INTEGER The order of the matrix B. N >= 0.``` [in] NCVT ``` NCVT is INTEGER The number of columns of the matrix VT. NCVT >= 0.``` [in] NRU ``` NRU is INTEGER The number of rows of the matrix U. NRU >= 0.``` [in] NCC ``` NCC is INTEGER The number of columns of the matrix C. NCC >= 0.``` [in,out] D ``` D is REAL array, dimension (N) On entry, the n diagonal elements of the bidiagonal matrix B. On exit, if INFO=0, the singular values of B in decreasing order.``` [in,out] E ``` E is REAL array, dimension (N-1) On entry, the N-1 offdiagonal elements of the bidiagonal matrix B. On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E will contain the diagonal and superdiagonal elements of a bidiagonal matrix orthogonally equivalent to the one given as input.``` [in,out] VT ``` VT is REAL array, dimension (LDVT, NCVT) On entry, an N-by-NCVT matrix VT. On exit, VT is overwritten by P**T * VT. Not referenced if NCVT = 0.``` [in] LDVT ``` LDVT is INTEGER The leading dimension of the array VT. LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.``` [in,out] U ``` U is REAL array, dimension (LDU, N) On entry, an NRU-by-N matrix U. On exit, U is overwritten by U * Q. Not referenced if NRU = 0.``` [in] LDU ``` LDU is INTEGER The leading dimension of the array U. LDU >= max(1,NRU).``` [in,out] C ``` C is REAL array, dimension (LDC, NCC) On entry, an N-by-NCC matrix C. On exit, C is overwritten by Q**T * C. Not referenced if NCC = 0.``` [in] LDC ``` LDC is INTEGER The leading dimension of the array C. LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.``` [out] WORK ` WORK is REAL array, dimension (4*N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: If INFO = -i, the i-th argument had an illegal value > 0: if NCVT = NRU = NCC = 0, = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 30*N iterations (in inner while loop) = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks) else NCVT = NRU = NCC = 0, the algorithm did not converge; D and E contain the elements of a bidiagonal matrix which is orthogonally similar to the input matrix B; if INFO = i, i elements of E have not converged to zero.```
Internal Parameters:
```  TOLMUL  REAL, default = max(10,min(100,EPS**(-1/8)))
TOLMUL controls the convergence criterion of the QR loop.
If it is positive, TOLMUL*EPS is the desired relative
precision in the computed singular values.
If it is negative, abs(TOLMUL*EPS*sigma_max) is the
desired absolute accuracy in the computed singular
values (corresponds to relative accuracy
abs(TOLMUL*EPS) in the largest singular value.
abs(TOLMUL) should be between 1 and 1/EPS, and preferably
between 10 (for fast convergence) and .1/EPS
(for there to be some accuracy in the results).
Default is to lose at either one eighth or 2 of the
available decimal digits in each computed singular value
(whichever is smaller).

MAXITR  INTEGER, default = 6
MAXITR controls the maximum number of passes of the
algorithm through its inner loop. The algorithms stops
(and so fails to converge) if the number of passes
through the inner loop exceeds MAXITR*N**2.```
Note:
```  Bug report from Cezary Dendek.
On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is
removed since it can overflow pretty easily (for N larger or equal
than 18,919). We instead use MAXITDIVN = MAXITR*N.```

Definition at line 238 of file sbdsqr.f.

240*
241* -- LAPACK computational routine --
242* -- LAPACK is a software package provided by Univ. of Tennessee, --
243* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
244*
245* .. Scalar Arguments ..
246 CHARACTER UPLO
247 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
248* ..
249* .. Array Arguments ..
250 REAL C( LDC, * ), D( * ), E( * ), U( LDU, * ),
251 \$ VT( LDVT, * ), WORK( * )
252* ..
253*
254* =====================================================================
255*
256* .. Parameters ..
257 REAL ZERO
258 parameter( zero = 0.0e0 )
259 REAL ONE
260 parameter( one = 1.0e0 )
261 REAL NEGONE
262 parameter( negone = -1.0e0 )
263 REAL HNDRTH
264 parameter( hndrth = 0.01e0 )
265 REAL TEN
266 parameter( ten = 10.0e0 )
267 REAL HNDRD
268 parameter( hndrd = 100.0e0 )
269 REAL MEIGTH
270 parameter( meigth = -0.125e0 )
271 INTEGER MAXITR
272 parameter( maxitr = 6 )
273* ..
274* .. Local Scalars ..
275 LOGICAL LOWER, ROTATE
276 INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
277 \$ MAXITDIVN, NM1, NM12, NM13, OLDLL, OLDM
278 REAL ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
279 \$ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
280 \$ SINR, SLL, SMAX, SMIN, SMINOA,
281 \$ SN, THRESH, TOL, TOLMUL, UNFL
282* ..
283* .. External Functions ..
284 LOGICAL LSAME
285 REAL SLAMCH
286 EXTERNAL lsame, slamch
287* ..
288* .. External Subroutines ..
289 EXTERNAL slartg, slas2, slasq1, slasr, slasv2, srot,
290 \$ sscal, sswap, xerbla
291* ..
292* .. Intrinsic Functions ..
293 INTRINSIC abs, max, min, real, sign, sqrt
294* ..
295* .. Executable Statements ..
296*
297* Test the input parameters.
298*
299 info = 0
300 lower = lsame( uplo, 'L' )
301 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
302 info = -1
303 ELSE IF( n.LT.0 ) THEN
304 info = -2
305 ELSE IF( ncvt.LT.0 ) THEN
306 info = -3
307 ELSE IF( nru.LT.0 ) THEN
308 info = -4
309 ELSE IF( ncc.LT.0 ) THEN
310 info = -5
311 ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
312 \$ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
313 info = -9
314 ELSE IF( ldu.LT.max( 1, nru ) ) THEN
315 info = -11
316 ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
317 \$ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
318 info = -13
319 END IF
320 IF( info.NE.0 ) THEN
321 CALL xerbla( 'SBDSQR', -info )
322 RETURN
323 END IF
324 IF( n.EQ.0 )
325 \$ RETURN
326 IF( n.EQ.1 )
327 \$ GO TO 160
328*
329* ROTATE is true if any singular vectors desired, false otherwise
330*
331 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
332*
333* If no singular vectors desired, use qd algorithm
334*
335 IF( .NOT.rotate ) THEN
336 CALL slasq1( n, d, e, work, info )
337*
338* If INFO equals 2, dqds didn't finish, try to finish
339*
340 IF( info .NE. 2 ) RETURN
341 info = 0
342 END IF
343*
344 nm1 = n - 1
345 nm12 = nm1 + nm1
346 nm13 = nm12 + nm1
347 idir = 0
348*
349* Get machine constants
350*
351 eps = slamch( 'Epsilon' )
352 unfl = slamch( 'Safe minimum' )
353*
354* If matrix lower bidiagonal, rotate to be upper bidiagonal
355* by applying Givens rotations on the left
356*
357 IF( lower ) THEN
358 DO 10 i = 1, n - 1
359 CALL slartg( d( i ), e( i ), cs, sn, r )
360 d( i ) = r
361 e( i ) = sn*d( i+1 )
362 d( i+1 ) = cs*d( i+1 )
363 work( i ) = cs
364 work( nm1+i ) = sn
365 10 CONTINUE
366*
367* Update singular vectors if desired
368*
369 IF( nru.GT.0 )
370 \$ CALL slasr( 'R', 'V', 'F', nru, n, work( 1 ), work( n ), u,
371 \$ ldu )
372 IF( ncc.GT.0 )
373 \$ CALL slasr( 'L', 'V', 'F', n, ncc, work( 1 ), work( n ), c,
374 \$ ldc )
375 END IF
376*
377* Compute singular values to relative accuracy TOL
378* (By setting TOL to be negative, algorithm will compute
379* singular values to absolute accuracy ABS(TOL)*norm(input matrix))
380*
381 tolmul = max( ten, min( hndrd, eps**meigth ) )
382 tol = tolmul*eps
383*
384* Compute approximate maximum, minimum singular values
385*
386 smax = zero
387 DO 20 i = 1, n
388 smax = max( smax, abs( d( i ) ) )
389 20 CONTINUE
390 DO 30 i = 1, n - 1
391 smax = max( smax, abs( e( i ) ) )
392 30 CONTINUE
393 smin = zero
394 IF( tol.GE.zero ) THEN
395*
396* Relative accuracy desired
397*
398 sminoa = abs( d( 1 ) )
399 IF( sminoa.EQ.zero )
400 \$ GO TO 50
401 mu = sminoa
402 DO 40 i = 2, n
403 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
404 sminoa = min( sminoa, mu )
405 IF( sminoa.EQ.zero )
406 \$ GO TO 50
407 40 CONTINUE
408 50 CONTINUE
409 sminoa = sminoa / sqrt( real( n ) )
410 thresh = max( tol*sminoa, maxitr*(n*(n*unfl)) )
411 ELSE
412*
413* Absolute accuracy desired
414*
415 thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
416 END IF
417*
418* Prepare for main iteration loop for the singular values
419* (MAXIT is the maximum number of passes through the inner
420* loop permitted before nonconvergence signalled.)
421*
422 maxitdivn = maxitr*n
423 iterdivn = 0
424 iter = -1
425 oldll = -1
426 oldm = -1
427*
428* M points to last element of unconverged part of matrix
429*
430 m = n
431*
432* Begin main iteration loop
433*
434 60 CONTINUE
435*
436* Check for convergence or exceeding iteration count
437*
438 IF( m.LE.1 )
439 \$ GO TO 160
440*
441 IF( iter.GE.n ) THEN
442 iter = iter - n
443 iterdivn = iterdivn + 1
444 IF( iterdivn.GE.maxitdivn )
445 \$ GO TO 200
446 END IF
447*
448* Find diagonal block of matrix to work on
449*
450 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
451 \$ d( m ) = zero
452 smax = abs( d( m ) )
453 DO 70 lll = 1, m - 1
454 ll = m - lll
455 abss = abs( d( ll ) )
456 abse = abs( e( ll ) )
457 IF( tol.LT.zero .AND. abss.LE.thresh )
458 \$ d( ll ) = zero
459 IF( abse.LE.thresh )
460 \$ GO TO 80
461 smax = max( smax, abss, abse )
462 70 CONTINUE
463 ll = 0
464 GO TO 90
465 80 CONTINUE
466 e( ll ) = zero
467*
468* Matrix splits since E(LL) = 0
469*
470 IF( ll.EQ.m-1 ) THEN
471*
473*
474 m = m - 1
475 GO TO 60
476 END IF
477 90 CONTINUE
478 ll = ll + 1
479*
480* E(LL) through E(M-1) are nonzero, E(LL-1) is zero
481*
482 IF( ll.EQ.m-1 ) THEN
483*
484* 2 by 2 block, handle separately
485*
486 CALL slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
487 \$ cosr, sinl, cosl )
488 d( m-1 ) = sigmx
489 e( m-1 ) = zero
490 d( m ) = sigmn
491*
492* Compute singular vectors, if desired
493*
494 IF( ncvt.GT.0 )
495 \$ CALL srot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt, cosr,
496 \$ sinr )
497 IF( nru.GT.0 )
498 \$ CALL srot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
499 IF( ncc.GT.0 )
500 \$ CALL srot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
501 \$ sinl )
502 m = m - 2
503 GO TO 60
504 END IF
505*
506* If working on new submatrix, choose shift direction
507* (from larger end diagonal element towards smaller)
508*
509 IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
510 IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
511*
512* Chase bulge from top (big end) to bottom (small end)
513*
514 idir = 1
515 ELSE
516*
517* Chase bulge from bottom (big end) to top (small end)
518*
519 idir = 2
520 END IF
521 END IF
522*
523* Apply convergence tests
524*
525 IF( idir.EQ.1 ) THEN
526*
527* Run convergence test in forward direction
528* First apply standard test to bottom of matrix
529*
530 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
531 \$ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
532 e( m-1 ) = zero
533 GO TO 60
534 END IF
535*
536 IF( tol.GE.zero ) THEN
537*
538* If relative accuracy desired,
539* apply convergence criterion forward
540*
541 mu = abs( d( ll ) )
542 smin = mu
543 DO 100 lll = ll, m - 1
544 IF( abs( e( lll ) ).LE.tol*mu ) THEN
545 e( lll ) = zero
546 GO TO 60
547 END IF
548 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
549 smin = min( smin, mu )
550 100 CONTINUE
551 END IF
552*
553 ELSE
554*
555* Run convergence test in backward direction
556* First apply standard test to top of matrix
557*
558 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
559 \$ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
560 e( ll ) = zero
561 GO TO 60
562 END IF
563*
564 IF( tol.GE.zero ) THEN
565*
566* If relative accuracy desired,
567* apply convergence criterion backward
568*
569 mu = abs( d( m ) )
570 smin = mu
571 DO 110 lll = m - 1, ll, -1
572 IF( abs( e( lll ) ).LE.tol*mu ) THEN
573 e( lll ) = zero
574 GO TO 60
575 END IF
576 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
577 smin = min( smin, mu )
578 110 CONTINUE
579 END IF
580 END IF
581 oldll = ll
582 oldm = m
583*
584* Compute shift. First, test if shifting would ruin relative
585* accuracy, and if so set the shift to zero.
586*
587 IF( tol.GE.zero .AND. n*tol*( smin / smax ).LE.
588 \$ max( eps, hndrth*tol ) ) THEN
589*
590* Use a zero shift to avoid loss of relative accuracy
591*
592 shift = zero
593 ELSE
594*
595* Compute the shift from 2-by-2 block at end of matrix
596*
597 IF( idir.EQ.1 ) THEN
598 sll = abs( d( ll ) )
599 CALL slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
600 ELSE
601 sll = abs( d( m ) )
602 CALL slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
603 END IF
604*
605* Test if shift negligible, and if so set to zero
606*
607 IF( sll.GT.zero ) THEN
608 IF( ( shift / sll )**2.LT.eps )
609 \$ shift = zero
610 END IF
611 END IF
612*
613* Increment iteration count
614*
615 iter = iter + m - ll
616*
617* If SHIFT = 0, do simplified QR iteration
618*
619 IF( shift.EQ.zero ) THEN
620 IF( idir.EQ.1 ) THEN
621*
622* Chase bulge from top to bottom
623* Save cosines and sines for later singular vector updates
624*
625 cs = one
626 oldcs = one
627 DO 120 i = ll, m - 1
628 CALL slartg( d( i )*cs, e( i ), cs, sn, r )
629 IF( i.GT.ll )
630 \$ e( i-1 ) = oldsn*r
631 CALL slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
632 work( i-ll+1 ) = cs
633 work( i-ll+1+nm1 ) = sn
634 work( i-ll+1+nm12 ) = oldcs
635 work( i-ll+1+nm13 ) = oldsn
636 120 CONTINUE
637 h = d( m )*cs
638 d( m ) = h*oldcs
639 e( m-1 ) = h*oldsn
640*
641* Update singular vectors
642*
643 IF( ncvt.GT.0 )
644 \$ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
645 \$ work( n ), vt( ll, 1 ), ldvt )
646 IF( nru.GT.0 )
647 \$ CALL slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
648 \$ work( nm13+1 ), u( 1, ll ), ldu )
649 IF( ncc.GT.0 )
650 \$ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
651 \$ work( nm13+1 ), c( ll, 1 ), ldc )
652*
653* Test convergence
654*
655 IF( abs( e( m-1 ) ).LE.thresh )
656 \$ e( m-1 ) = zero
657*
658 ELSE
659*
660* Chase bulge from bottom to top
661* Save cosines and sines for later singular vector updates
662*
663 cs = one
664 oldcs = one
665 DO 130 i = m, ll + 1, -1
666 CALL slartg( d( i )*cs, e( i-1 ), cs, sn, r )
667 IF( i.LT.m )
668 \$ e( i ) = oldsn*r
669 CALL slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
670 work( i-ll ) = cs
671 work( i-ll+nm1 ) = -sn
672 work( i-ll+nm12 ) = oldcs
673 work( i-ll+nm13 ) = -oldsn
674 130 CONTINUE
675 h = d( ll )*cs
676 d( ll ) = h*oldcs
677 e( ll ) = h*oldsn
678*
679* Update singular vectors
680*
681 IF( ncvt.GT.0 )
682 \$ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
683 \$ work( nm13+1 ), vt( ll, 1 ), ldvt )
684 IF( nru.GT.0 )
685 \$ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
686 \$ work( n ), u( 1, ll ), ldu )
687 IF( ncc.GT.0 )
688 \$ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
689 \$ work( n ), c( ll, 1 ), ldc )
690*
691* Test convergence
692*
693 IF( abs( e( ll ) ).LE.thresh )
694 \$ e( ll ) = zero
695 END IF
696 ELSE
697*
698* Use nonzero shift
699*
700 IF( idir.EQ.1 ) THEN
701*
702* Chase bulge from top to bottom
703* Save cosines and sines for later singular vector updates
704*
705 f = ( abs( d( ll ) )-shift )*
706 \$ ( sign( one, d( ll ) )+shift / d( ll ) )
707 g = e( ll )
708 DO 140 i = ll, m - 1
709 CALL slartg( f, g, cosr, sinr, r )
710 IF( i.GT.ll )
711 \$ e( i-1 ) = r
712 f = cosr*d( i ) + sinr*e( i )
713 e( i ) = cosr*e( i ) - sinr*d( i )
714 g = sinr*d( i+1 )
715 d( i+1 ) = cosr*d( i+1 )
716 CALL slartg( f, g, cosl, sinl, r )
717 d( i ) = r
718 f = cosl*e( i ) + sinl*d( i+1 )
719 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
720 IF( i.LT.m-1 ) THEN
721 g = sinl*e( i+1 )
722 e( i+1 ) = cosl*e( i+1 )
723 END IF
724 work( i-ll+1 ) = cosr
725 work( i-ll+1+nm1 ) = sinr
726 work( i-ll+1+nm12 ) = cosl
727 work( i-ll+1+nm13 ) = sinl
728 140 CONTINUE
729 e( m-1 ) = f
730*
731* Update singular vectors
732*
733 IF( ncvt.GT.0 )
734 \$ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
735 \$ work( n ), vt( ll, 1 ), ldvt )
736 IF( nru.GT.0 )
737 \$ CALL slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
738 \$ work( nm13+1 ), u( 1, ll ), ldu )
739 IF( ncc.GT.0 )
740 \$ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
741 \$ work( nm13+1 ), c( ll, 1 ), ldc )
742*
743* Test convergence
744*
745 IF( abs( e( m-1 ) ).LE.thresh )
746 \$ e( m-1 ) = zero
747*
748 ELSE
749*
750* Chase bulge from bottom to top
751* Save cosines and sines for later singular vector updates
752*
753 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
754 \$ d( m ) )
755 g = e( m-1 )
756 DO 150 i = m, ll + 1, -1
757 CALL slartg( f, g, cosr, sinr, r )
758 IF( i.LT.m )
759 \$ e( i ) = r
760 f = cosr*d( i ) + sinr*e( i-1 )
761 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
762 g = sinr*d( i-1 )
763 d( i-1 ) = cosr*d( i-1 )
764 CALL slartg( f, g, cosl, sinl, r )
765 d( i ) = r
766 f = cosl*e( i-1 ) + sinl*d( i-1 )
767 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
768 IF( i.GT.ll+1 ) THEN
769 g = sinl*e( i-2 )
770 e( i-2 ) = cosl*e( i-2 )
771 END IF
772 work( i-ll ) = cosr
773 work( i-ll+nm1 ) = -sinr
774 work( i-ll+nm12 ) = cosl
775 work( i-ll+nm13 ) = -sinl
776 150 CONTINUE
777 e( ll ) = f
778*
779* Test convergence
780*
781 IF( abs( e( ll ) ).LE.thresh )
782 \$ e( ll ) = zero
783*
784* Update singular vectors if desired
785*
786 IF( ncvt.GT.0 )
787 \$ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
788 \$ work( nm13+1 ), vt( ll, 1 ), ldvt )
789 IF( nru.GT.0 )
790 \$ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
791 \$ work( n ), u( 1, ll ), ldu )
792 IF( ncc.GT.0 )
793 \$ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
794 \$ work( n ), c( ll, 1 ), ldc )
795 END IF
796 END IF
797*
798* QR iteration finished, go back and check convergence
799*
800 GO TO 60
801*
802* All singular values converged, so make them positive
803*
804 160 CONTINUE
805 DO 170 i = 1, n
806 IF( d( i ).LT.zero ) THEN
807 d( i ) = -d( i )
808*
809* Change sign of singular vectors, if desired
810*
811 IF( ncvt.GT.0 )
812 \$ CALL sscal( ncvt, negone, vt( i, 1 ), ldvt )
813 END IF
814 170 CONTINUE
815*
816* Sort the singular values into decreasing order (insertion sort on
817* singular values, but only one transposition per singular vector)
818*
819 DO 190 i = 1, n - 1
820*
821* Scan for smallest D(I)
822*
823 isub = 1
824 smin = d( 1 )
825 DO 180 j = 2, n + 1 - i
826 IF( d( j ).LE.smin ) THEN
827 isub = j
828 smin = d( j )
829 END IF
830 180 CONTINUE
831 IF( isub.NE.n+1-i ) THEN
832*
833* Swap singular values and vectors
834*
835 d( isub ) = d( n+1-i )
836 d( n+1-i ) = smin
837 IF( ncvt.GT.0 )
838 \$ CALL sswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
839 \$ ldvt )
840 IF( nru.GT.0 )
841 \$ CALL sswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
842 IF( ncc.GT.0 )
843 \$ CALL sswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
844 END IF
845 190 CONTINUE
846 GO TO 220
847*
848* Maximum number of iterations exceeded, failure to converge
849*
850 200 CONTINUE
851 info = 0
852 DO 210 i = 1, n - 1
853 IF( e( i ).NE.zero )
854 \$ info = info + 1
855 210 CONTINUE
856 220 CONTINUE
857 RETURN
858*
859* End of SBDSQR
860*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
subroutine slas2(f, g, h, ssmin, ssmax)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition slas2.f:105
subroutine slasq1(n, d, e, work, info)
SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
Definition slasq1.f:108
subroutine slasr(side, pivot, direct, m, n, c, s, a, lda)
SLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition slasr.f:199
subroutine slasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition slasv2.f:136
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: