slasq1.f

Go to the documentation of this file.

*> \brief \b SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download SLASQ1 + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq1.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq1.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq1.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE SLASQ1( N, D, E, WORK, INFO )
* 
*       .. Scalar Arguments ..
*       INTEGER            INFO, N
*       ..
*       .. Array Arguments ..
*       REAL               D( * ), E( * ), WORK( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SLASQ1 computes the singular values of a real N-by-N bidiagonal
*> matrix with diagonal D and off-diagonal E. The singular values
*> are computed to high relative accuracy, in the absence of
*> denormalization, underflow and overflow. The algorithm was first
*> presented in
*>
*> "Accurate singular values and differential qd algorithms" by K. V.
*> Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
*> 1994,
*>
*> and the present implementation is described in "An implementation of
*> the dqds Algorithm (Positive Case)", LAPACK Working Note.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>        The number of rows and columns in the matrix. N >= 0.
*> \endverbatim
*>
*> \param[in,out] D
*> \verbatim
*>          D is REAL array, dimension (N)
*>        On entry, D contains the diagonal elements of the
*>        bidiagonal matrix whose SVD is desired. On normal exit,
*>        D contains the singular values in decreasing order.
*> \endverbatim
*>
*> \param[in,out] E
*> \verbatim
*>          E is REAL array, dimension (N)
*>        On entry, elements E(1:N-1) contain the off-diagonal elements
*>        of the bidiagonal matrix whose SVD is desired.
*>        On exit, E is overwritten.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (4*N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>        = 0: successful exit
*>        < 0: if INFO = -i, the i-th argument had an illegal value
*>        > 0: the algorithm failed
*>             = 1, a split was marked by a positive value in E
*>             = 2, current block of Z not diagonalized after 100*N
*>                  iterations (in inner while loop)  On exit D and E
*>                  represent a matrix with the same singular values
*>                  which the calling subroutine could use to finish the
*>                  computation, or even feed back into SLASQ1
*>             = 3, termination criterion of outer while loop not met 
*>                  (program created more than N unreduced blocks)
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2015
*
*> \ingroup auxOTHERcomputational
*
*  =====================================================================
      SUBROUTINE slasq1( N, D, E, WORK, INFO )
*
*  -- LAPACK computational routine (version 3.6.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2015
*
*     .. Scalar Arguments ..
      INTEGER            INFO, N
*     ..
*     .. Array Arguments ..
      REAL               D( * ), E( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      parameter                ( zero = 0.0e0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IINFO
      REAL               EPS, SCALE, SAFMIN, SIGMN, SIGMX
*     ..
*     .. External Subroutines ..
      EXTERNAL           scopy, slas2, slascl, slasq2, slasrt, xerbla
*     ..
*     .. External Functions ..
      REAL               SLAMCH
      EXTERNAL           slamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, sqrt
*     ..
*     .. Executable Statements ..
*
      info = 0
      IF( n.LT.0 ) THEN
         info = -1
         CALL xerbla( 'SLASQ1', -info )
         RETURN
      ELSE IF( n.EQ.0 ) THEN
         RETURN
      ELSE IF( n.EQ.1 ) THEN
         d( 1 ) = abs( d( 1 ) )
         RETURN
      ELSE IF( n.EQ.2 ) THEN
         CALL slas2( d( 1 ), e( 1 ), d( 2 ), sigmn, sigmx )
         d( 1 ) = sigmx
         d( 2 ) = sigmn
         RETURN
      END IF
*
*     Estimate the largest singular value.
*
      sigmx = zero
      DO 10 i = 1, n - 1
         d( i ) = abs( d( i ) )
         sigmx = max( sigmx, abs( e( i ) ) )
   10 CONTINUE
      d( n ) = abs( d( n ) )
*
*     Early return if SIGMX is zero (matrix is already diagonal).
*
      IF( sigmx.EQ.zero ) THEN
         CALL slasrt( 'D', n, d, iinfo )
         RETURN
      END IF
*
      DO 20 i = 1, n
         sigmx = max( sigmx, d( i ) )
   20 CONTINUE
*
*     Copy D and E into WORK (in the Z format) and scale (squaring the
*     input data makes scaling by a power of the radix pointless).
*
      eps = slamch( 'Precision' )
      safmin = slamch( 'Safe minimum' )
      scale = sqrt( eps / safmin )
      CALL scopy( n, d, 1, work( 1 ), 2 )
      CALL scopy( n-1, e, 1, work( 2 ), 2 )
      CALL slascl( 'G', 0, 0, sigmx, scale, 2*n-1, 1, work, 2*n-1,
     $             iinfo )
*         
*     Compute the q's and e's.
*
      DO 30 i = 1, 2*n - 1
         work( i ) = work( i )**2
   30 CONTINUE
      work( 2*n ) = zero
*
      CALL slasq2( n, work, info )
*
      IF( info.EQ.0 ) THEN
         DO 40 i = 1, n
            d( i ) = sqrt( work( i ) )
   40    CONTINUE
         CALL slascl( 'G', 0, 0, scale, sigmx, n, 1, d, n, iinfo )
      ELSE IF( info.EQ.2 ) THEN
*
*     Maximum number of iterations exceeded.  Move data from WORK
*     into D and E so the calling subroutine can try to finish
*
         DO i = 1, n
            d( i ) = sqrt( work( 2*i-1 ) )
            e( i ) = sqrt( work( 2*i ) )
         END DO
         CALL slascl( 'G', 0, 0, scale, sigmx, n, 1, d, n, iinfo )
         CALL slascl( 'G', 0, 0, scale, sigmx, n, 1, e, n, iinfo )
      END IF
*
      RETURN
*
*     End of SLASQ1
*
      END