*> \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 *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \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 December 2016 * *> \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 (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * * .. 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