#include "ESMF_LapackBlas.inc" !> \brief \b DLASV2 ! ! =========== 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. ! !> \date November 2011 ! !> \ingroup auxOTHERauxiliary ! !> \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 (version 3.4.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! November 2011 ! ! .. 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 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