dlasv2.f

Go to the documentation of this file.

*> \brief \b DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLASV2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasv2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasv2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasv2.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
*
*       .. Scalar Arguments ..
*       DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLASV2 computes the singular value decomposition of a 2-by-2
*> triangular matrix
*>    [  F   G  ]
*>    [  0   H  ].
*> On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
*> smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
*> right singular vectors for abs(SSMAX), giving the decomposition
*>
*>    [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
*>    [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] F
*> \verbatim
*>          F is DOUBLE PRECISION
*>          The (1,1) element of the 2-by-2 matrix.
*> \endverbatim
*>
*> \param[in] G
*> \verbatim
*>          G is DOUBLE PRECISION
*>          The (1,2) element of the 2-by-2 matrix.
*> \endverbatim
*>
*> \param[in] H
*> \verbatim
*>          H is DOUBLE PRECISION
*>          The (2,2) element of the 2-by-2 matrix.
*> \endverbatim
*>
*> \param[out] SSMIN
*> \verbatim
*>          SSMIN is DOUBLE PRECISION
*>          abs(SSMIN) is the smaller singular value.
*> \endverbatim
*>
*> \param[out] SSMAX
*> \verbatim
*>          SSMAX is DOUBLE PRECISION
*>          abs(SSMAX) is the larger singular value.
*> \endverbatim
*>
*> \param[out] SNL
*> \verbatim
*>          SNL is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[out] CSL
*> \verbatim
*>          CSL is DOUBLE PRECISION
*>          The vector (CSL, SNL) is a unit left singular vector for the
*>          singular value abs(SSMAX).
*> \endverbatim
*>
*> \param[out] SNR
*> \verbatim
*>          SNR is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[out] CSR
*> \verbatim
*>          CSR is DOUBLE PRECISION
*>          The vector (CSR, SNR) is a unit right singular vector for the
*>          singular value abs(SSMAX).
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup OTHERauxiliary
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  Any input parameter may be aliased with any output parameter.
*>
*>  Barring over/underflow and assuming a guard digit in subtraction, all
*>  output quantities are correct to within a few units in the last
*>  place (ulps).
*>
*>  In IEEE arithmetic, the code works correctly if one matrix element is
*>  infinite.
*>
*>  Overflow will not occur unless the largest singular value itself
*>  overflows or is within a few ulps of overflow. (On machines with
*>  partial overflow, like the Cray, overflow may occur if the largest
*>  singular value is within a factor of 2 of overflow.)
*>
*>  Underflow is harmless if underflow is gradual. Otherwise, results
*>  may correspond to a matrix modified by perturbations of size near
*>  the underflow threshold.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE dlasv2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      parameter( zero = 0.0d0 )
      DOUBLE PRECISION   HALF
      parameter( half = 0.5d0 )
      DOUBLE PRECISION   ONE
      parameter( one = 1.0d0 )
      DOUBLE PRECISION   TWO
      parameter( two = 2.0d0 )
      DOUBLE PRECISION   FOUR
      parameter( four = 4.0d0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            GASMAL, SWAP
      INTEGER            PMAX
      DOUBLE PRECISION   A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
     $                   MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, sign, sqrt
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           dlamch
*     ..
*     .. Executable Statements ..
*
      ft = f
      fa = abs( ft )
      ht = h
      ha = abs( h )
*
*     PMAX points to the maximum absolute element of matrix
*       PMAX = 1 if F largest in absolute values
*       PMAX = 2 if G largest in absolute values
*       PMAX = 3 if H largest in absolute values
*
      pmax = 1
      swap = ( ha.GT.fa )
      IF( swap ) THEN
         pmax = 3
         temp = ft
         ft = ht
         ht = temp
         temp = fa
         fa = ha
         ha = temp
*
*        Now FA .ge. HA
*
      END IF
      gt = g
      ga = abs( gt )
      IF( ga.EQ.zero ) THEN
*
*        Diagonal matrix
*
         ssmin = ha
         ssmax = fa
         clt = one
         crt = one
         slt = zero
         srt = zero
      ELSE
         gasmal = .true.
         IF( ga.GT.fa ) THEN
            pmax = 2
            IF( ( fa / ga ).LT.dlamch( 'eps' ) ) THEN
*
*              Case of very large GA
*
               GASMAL = .FALSE.
               SSMAX = GA
.GT.               IF( HAONE ) THEN
                  SSMIN = FA / ( GA / HA )
               ELSE
                  SSMIN = ( FA / GA )*HA
               END IF
               CLT = ONE
               SLT = HT / GT
               SRT = ONE
               CRT = FT / GT
            END IF
         END IF
         IF( GASMAL ) THEN
*
*           Normal case
*
            D = FA - HA
.EQ.            IF( DFA ) THEN
*
*              Copes with infinite F or H
*
               L = ONE
            ELSE
               L = D / FA
            END IF
*
*           Note that 0 .le. L .le. 1
*
            M = GT / FT
*
*           Note that abs(M) .le. 1/macheps
*
            T = TWO - L
*
*           Note that T .ge. 1
*
            MM = M*M
            TT = T*T
            S = SQRT( TT+MM )
*
*           Note that 1 .le. S .le. 1 + 1/macheps
*
.EQ.            IF( LZERO ) THEN
               R = ABS( M )
            ELSE
               R = SQRT( L*L+MM )
            END IF
*
*           Note that 0 .le. R .le. 1 + 1/macheps
*
            A = HALF*( S+R )
*
*           Note that 1 .le. A .le. 1 + abs(M)
*
            SSMIN = HA / A
            SSMAX = FA*A
.EQ.            IF( MMZERO ) THEN
*
*              Note that M is very tiny
*
.EQ.               IF( LZERO ) THEN
                  T = SIGN( TWO, FT )*SIGN( ONE, GT )
               ELSE
                  T = GT / SIGN( D, FT ) + M / T
               END IF
            ELSE
               T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
            END IF
            L = SQRT( T*T+FOUR )
            CRT = TWO / L
            SRT = T / L
            CLT = ( CRT+SRT*M ) / A
            SLT = ( HT / FT )*SRT / A
         END IF
      END IF
      IF( SWAP ) THEN
         CSL = SRT
         SNL = CRT
         CSR = SLT
         SNR = CLT
      ELSE
         CSL = CLT
         SNL = SLT
         CSR = CRT
         SNR = SRT
      END IF
*
*     Correct signs of SSMAX and SSMIN
*
.EQ.      IF( PMAX1 )
     $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
.EQ.      IF( PMAX2 )
     $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
.EQ.      IF( PMAX3 )
     $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
      SSMAX = SIGN( SSMAX, TSIGN )
      SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
      RETURN
*
*     End of DLASV2
*
      END