1: !> \brief \b DROTG
2: !
3: ! =========== DOCUMENTATION ===========
4: !
5: ! Online html documentation available at
6: ! http://www.netlib.org/lapack/explore-html/
7: !
8: ! Definition:
9: ! ===========
10: !
11: ! DROTG constructs a plane rotation
12: ! [ c s ] [ a ] = [ r ]
13: ! [ -s c ] [ b ] [ 0 ]
14: ! satisfying c**2 + s**2 = 1.
15: !
16: !> \par Purpose:
17: ! =============
18: !>
19: !> \verbatim
20: !>
21: !> The computation uses the formulas
22: !> sigma = sgn(a) if |a| > |b|
23: !> = sgn(b) if |b| >= |a|
24: !> r = sigma*sqrt( a**2 + b**2 )
25: !> c = 1; s = 0 if r = 0
26: !> c = a/r; s = b/r if r != 0
27: !> The subroutine also computes
28: !> z = s if |a| > |b|,
29: !> = 1/c if |b| >= |a| and c != 0
30: !> = 1 if c = 0
31: !> This allows c and s to be reconstructed from z as follows:
32: !> If z = 1, set c = 0, s = 1.
33: !> If |z| < 1, set c = sqrt(1 - z**2) and s = z.
34: !> If |z| > 1, set c = 1/z and s = sqrt( 1 - c**2).
35: !>
36: !> \endverbatim
37: !
38: ! Arguments:
39: ! ==========
40: !
41: !> \param[in,out] A
42: !> \verbatim
43: !> A is DOUBLE PRECISION
44: !> On entry, the scalar a.
45: !> On exit, the scalar r.
46: !> \endverbatim
47: !>
48: !> \param[in,out] B
49: !> \verbatim
50: !> B is DOUBLE PRECISION
51: !> On entry, the scalar b.
52: !> On exit, the scalar z.
53: !> \endverbatim
54: !>
55: !> \param[out] C
56: !> \verbatim
57: !> C is DOUBLE PRECISION
58: !> The scalar c.
59: !> \endverbatim
60: !>
61: !> \param[out] S
62: !> \verbatim
63: !> S is DOUBLE PRECISION
64: !> The scalar s.
65: !> \endverbatim
66: !
67: ! Authors:
68: ! ========
69: !
70: !> \author Edward Anderson, Lockheed Martin
71: !
72: !> \par Contributors:
73: ! ==================
74: !>
75: !> Weslley Pereira, University of Colorado Denver, USA
76: !
77: !> \ingroup single_blas_level1
78: !
79: !> \par Further Details:
80: ! =====================
81: !>
82: !> \verbatim
83: !>
84: !> Anderson E. (2017)
85: !> Algorithm 978: Safe Scaling in the Level 1 BLAS
86: !> ACM Trans Math Softw 44:1--28
87: !> https://doi.org/10.1145/3061665
88: !>
89: !> \endverbatim
90: !
91: ! =====================================================================
92: subroutine DROTG( a, b, c, s )
93: integer, parameter :: wp = kind(1.d0)
94: !
95: ! -- Reference BLAS level1 routine --
96: ! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
97: ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
98: !
99: ! .. Constants ..
100: real(wp), parameter :: zero = 0.0_wp
101: real(wp), parameter :: one = 1.0_wp
102: ! ..
103: ! .. Scaling constants ..
104: real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( &
105: minexponent(0._wp)-1, &
106: 1-maxexponent(0._wp) &
107: )
108: real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( &
109: 1-minexponent(0._wp), &
110: maxexponent(0._wp)-1 &
111: )
112: ! ..
113: ! .. Scalar Arguments ..
114: real(wp) :: a, b, c, s
115: ! ..
116: ! .. Local Scalars ..
117: real(wp) :: anorm, bnorm, scl, sigma, r, z
118: ! ..
119: anorm = abs(a)
120: bnorm = abs(b)
121: if( bnorm == zero ) then
122: c = one
123: s = zero
124: b = zero
125: else if( anorm == zero ) then
126: c = zero
127: s = one
128: a = b
129: b = one
130: else
131: scl = min( safmax, max( safmin, anorm, bnorm ) )
132: if( anorm > bnorm ) then
133: sigma = sign(one,a)
134: else
135: sigma = sign(one,b)
136: end if
137: r = sigma*( scl*sqrt((a/scl)**2 + (b/scl)**2) )
138: c = a/r
139: s = b/r
140: if( anorm > bnorm ) then
141: z = s
142: else if( c /= zero ) then
143: z = one/c
144: else
145: z = one
146: end if
147: a = r
148: b = z
149: end if
150: return
151: end subroutine
CVSweb interface <joel.bertrand@systella.fr>