1!> \brief \b ZLARTG generates a plane rotation with real cosine and complex sine.
2!
3!  =========== DOCUMENTATION ===========
4!
5! Online html documentation available at
6!            http://www.netlib.org/lapack/explore-html/
7!
8!  Definition:
9!  ===========
10!
11!       SUBROUTINE ZLARTG( F, G, C, S, R )
12!
13!       .. Scalar Arguments ..
14!       REAL(wp)           C
15!       COMPLEX(wp)        F, G, R, S
16!       ..
17!
18!> \par Purpose:
19!  =============
20!>
21!> \verbatim
22!>
23!> ZLARTG generates a plane rotation so that
24!>
25!>    [  C         S  ] . [ F ]  =  [ R ]
26!>    [ -conjg(S)  C  ]   [ G ]     [ 0 ]
27!>
28!> where C is real and C**2 + |S|**2 = 1.
29!>
30!> The mathematical formulas used for C and S are
31!>
32!>    sgn(x) = {  x / |x|,   x != 0
33!>             {  1,         x = 0
34!>
35!>    R = sgn(F) * sqrt(|F|**2 + |G|**2)
36!>
37!>    C = |F| / sqrt(|F|**2 + |G|**2)
38!>
39!>    S = sgn(F) * conjg(G) / sqrt(|F|**2 + |G|**2)
40!>
41!> When F and G are real, the formulas simplify to C = F/R and
42!> S = G/R, and the returned values of C, S, and R should be
43!> identical to those returned by DLARTG.
44!>
45!> The algorithm used to compute these quantities incorporates scaling
46!> to avoid overflow or underflow in computing the square root of the
47!> sum of squares.
48!>
49!> This is a faster version of the BLAS1 routine ZROTG, except for
50!> the following differences:
51!>    F and G are unchanged on return.
52!>    If G=0, then C=1 and S=0.
53!>    If F=0, then C=0 and S is chosen so that R is real.
54!>
55!> Below, wp=>dp stands for double precision from LA_CONSTANTS module.
56!> \endverbatim
57!
58!  Arguments:
59!  ==========
60!
61!> \param[in] F
62!> \verbatim
63!>          F is COMPLEX(wp)
64!>          The first component of vector to be rotated.
65!> \endverbatim
66!>
67!> \param[in] G
68!> \verbatim
69!>          G is COMPLEX(wp)
70!>          The second component of vector to be rotated.
71!> \endverbatim
72!>
73!> \param[out] C
74!> \verbatim
75!>          C is REAL(wp)
76!>          The cosine of the rotation.
77!> \endverbatim
78!>
79!> \param[out] S
80!> \verbatim
81!>          S is COMPLEX(wp)
82!>          The sine of the rotation.
83!> \endverbatim
84!>
85!> \param[out] R
86!> \verbatim
87!>          R is COMPLEX(wp)
88!>          The nonzero component of the rotated vector.
89!> \endverbatim
90!
91!  Authors:
92!  ========
93!
94!> \author Edward Anderson, Lockheed Martin
95!
96!> \date August 2016
97!
98!> \ingroup OTHERauxiliary
99!
100!> \par Contributors:
101!  ==================
102!>
103!> Weslley Pereira, University of Colorado Denver, USA
104!
105!> \par Further Details:
106!  =====================
107!>
108!> \verbatim
109!>
110!>  Anderson E. (2017)
111!>  Algorithm 978: Safe Scaling in the Level 1 BLAS
112!>  ACM Trans Math Softw 44:1--28
113!>  https://doi.org/10.1145/3061665
114!>
115!> \endverbatim
116!
117subroutine ZLARTG( f, g, c, s, r )
118   use LA_CONSTANTS, &
119   only: wp=>dp, zero=>dzero, one=>done, two=>dtwo, czero=>zzero, &
120         rtmin=>drtmin, rtmax=>drtmax, safmin=>dsafmin, safmax=>dsafmax
121!
122!  -- LAPACK auxiliary routine (version 3.10.0) --
123!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
124!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
125!     February 2021
126!
127!  .. Scalar Arguments ..
128   real(wp)           c
129   complex(wp)        f, g, r, s
130!  ..
131!  .. Local Scalars ..
132   real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w
133   complex(wp) :: fs, gs, t
134!  ..
135!  .. Intrinsic Functions ..
136   intrinsic :: abs, aimag, conjg, max, min, real, sqrt
137!  ..
138!  .. Statement Functions ..
139   real(wp) :: ABSSQ
140!  ..
141!  .. Statement Function definitions ..
142   ABSSQ( t ) = real( t )**2 + aimag( t )**2
143!  ..
144!  .. Executable Statements ..
145!
146   if( g == czero ) then
147      c = one
148      s = czero
149      r = f
150   else if( f == czero ) then
151      c = zero
152      g1 = max( abs(real(g)), abs(aimag(g)) )
153      if( g1 > rtmin .and. g1 < rtmax ) then
154!
155!        Use unscaled algorithm
156!
157         g2 = ABSSQ( g )
158         d = sqrt( g2 )
159         s = conjg( g ) / d
160         r = d
161      else
162!
163!        Use scaled algorithm
164!
165         u = min( safmax, max( safmin, g1 ) )
166         uu = one / u
167         gs = g*uu
168         g2 = ABSSQ( gs )
169         d = sqrt( g2 )
170         s = conjg( gs ) / d
171         r = d*u
172      end if
173   else
174      f1 = max( abs(real(f)), abs(aimag(f)) )
175      g1 = max( abs(real(g)), abs(aimag(g)) )
176      if( f1 > rtmin .and. f1 < rtmax .and. &
177          g1 > rtmin .and. g1 < rtmax ) then
178!
179!        Use unscaled algorithm
180!
181         f2 = ABSSQ( f )
182         g2 = ABSSQ( g )
183         h2 = f2 + g2
184         if( f2 > rtmin .and. h2 < rtmax ) then
185            d = sqrt( f2*h2 )
186         else
187            d = sqrt( f2 )*sqrt( h2 )
188         end if
189         p = 1 / d
190         c = f2*p
191         s = conjg( g )*( f*p )
192         r = f*( h2*p )
193      else
194!
195!        Use scaled algorithm
196!
197         u = min( safmax, max( safmin, f1, g1 ) )
198         uu = one / u
199         gs = g*uu
200         g2 = ABSSQ( gs )
201         if( f1*uu < rtmin ) then
202!
203!           f is not well-scaled when scaled by g1.
204!           Use a different scaling for f.
205!
206            v = min( safmax, max( safmin, f1 ) )
207            vv = one / v
208            w = v * uu
209            fs = f*vv
210            f2 = ABSSQ( fs )
211            h2 = f2*w**2 + g2
212         else
213!
214!           Otherwise use the same scaling for f and g.
215!
216            w = one
217            fs = f*uu
218            f2 = ABSSQ( fs )
219            h2 = f2 + g2
220         end if
221         if( f2 > rtmin .and. h2 < rtmax ) then
222            d = sqrt( f2*h2 )
223         else
224            d = sqrt( f2 )*sqrt( h2 )
225         end if
226         p = 1 / d
227         c = ( f2*p )*w
228         s = conjg( gs )*( fs*p )
229         r = ( fs*( h2*p ) )*u
230      end if
231   end if
232   return
233end subroutine
234