Chimie Théorique » OpenPath

      SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )

*

*  -- LAPACK auxiliary routine (version 3.2) --

*  -- LAPACK is a software package provided by Univ. of Tennessee,    --

*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--

*     November 2006

*

*     .. Scalar Arguments ..

      DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN

*     ..

*

*  Purpose

*  =======

*

*  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 ].

*

*  Arguments

*  =========

*

*  F       (input) DOUBLE PRECISION

*          The (1,1) element of the 2-by-2 matrix.

*

*  G       (input) DOUBLE PRECISION

*          The (1,2) element of the 2-by-2 matrix.

*

*  H       (input) DOUBLE PRECISION

*          The (2,2) element of the 2-by-2 matrix.

*

*  SSMIN   (output) DOUBLE PRECISION

*          abs(SSMIN) is the smaller singular value.

*

*  SSMAX   (output) DOUBLE PRECISION

*          abs(SSMAX) is the larger singular value.

*

*  SNL     (output) DOUBLE PRECISION

*  CSL     (output) DOUBLE PRECISION

*          The vector (CSL, SNL) is a unit left singular vector for the

*          singular value abs(SSMAX).

*

*  SNR     (output) DOUBLE PRECISION

*  CSR     (output) DOUBLE PRECISION

*          The vector (CSR, SNR) is a unit right singular vector for the

*          singular value abs(SSMAX).

*

*  Further Details

*  ===============

*

*  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.

*

* =====================================================================

*

*     .. 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

               IF( HA.GT.ONE ) 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

            IF( D.EQ.FA ) 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

*

            IF( L.EQ.ZERO ) 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

            IF( MM.EQ.ZERO ) THEN

*

*              Note that M is very tiny

*

               IF( L.EQ.ZERO ) 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

*

      IF( PMAX.EQ.1 )

     $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )

      IF( PMAX.EQ.2 )

     $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )

      IF( PMAX.EQ.3 )

     $   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