Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65229464
dlasv2.f
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Jun 2, 01:18
Size
8 KB
Mime Type
text/html
Expires
Tue, Jun 4, 01:18 (2 d)
Engine
blob
Format
Raw Data
Handle
17966275
Attached To
rLAMMPS lammps
dlasv2.f
View Options
*> \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.
*
*> \date September 2012
*
*> \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.2) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* September 2012
*
* .. 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
Event Timeline
Log In to Comment