1!> \brief \b ZLASSQ updates a sum of squares represented in scaled form.
2!
3!  =========== DOCUMENTATION ===========
4!
5! Online html documentation available at
6!            http://www.netlib.org/lapack/explore-html/
7!
8!> \htmlonly
9!> Download ZLASSQ + dependencies
10!> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlassq.f90">
11!> [TGZ]</a>
12!> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlassq.f90">
13!> [ZIP]</a>
14!> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlassq.f90">
15!> [TXT]</a>
16!> \endhtmlonly
17!
18!  Definition:
19!  ===========
20!
21!       SUBROUTINE ZLASSQ( N, X, INCX, SCALE, SUMSQ )
22!
23!       .. Scalar Arguments ..
24!       INTEGER            INCX, N
25!       DOUBLE PRECISION   SCALE, SUMSQ
26!       ..
27!       .. Array Arguments ..
28!       DOUBLE COMPLEX     X( * )
29!       ..
30!
31!
32!> \par Purpose:
33!  =============
34!>
35!> \verbatim
36!>
37!> ZLASSQ  returns the values  scl  and  smsq  such that
38!>
39!>    ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,
40!>
41!> where  x( i ) = X( 1 + ( i - 1 )*INCX ). The value of  sumsq  is
42!> assumed to be non-negative and  scl  returns the value
43!>
44!>    scl = max( scale, abs( x( i ) ) ).
45!>
46!> scale and sumsq must be supplied in SCALE and SUMSQ and
47!> scl and smsq are overwritten on SCALE and SUMSQ respectively.
48!>
49!> \endverbatim
50!
51!  Arguments:
52!  ==========
53!
54!> \param[in] N
55!> \verbatim
56!>          N is INTEGER
57!>          The number of elements to be used from the vector x.
58!> \endverbatim
59!>
60!> \param[in] X
61!> \verbatim
62!>          X is DOUBLE COMPLEX array, dimension (1+(N-1)*abs(INCX))
63!>          The vector for which a scaled sum of squares is computed.
64!>             x( i )  = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
65!> \endverbatim
66!>
67!> \param[in] INCX
68!> \verbatim
69!>          INCX is INTEGER
70!>          The increment between successive values of the vector x.
71!>          If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
72!>          If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
73!>          If INCX = 0, x isn't a vector so there is no need to call
74!>          this subroutine.  If you call it anyway, it will count x(1)
75!>          in the vector norm N times.
76!> \endverbatim
77!>
78!> \param[in,out] SCALE
79!> \verbatim
80!>          SCALE is DOUBLE PRECISION
81!>          On entry, the value  scale  in the equation above.
82!>          On exit, SCALE is overwritten with  scl , the scaling factor
83!>          for the sum of squares.
84!> \endverbatim
85!>
86!> \param[in,out] SUMSQ
87!> \verbatim
88!>          SUMSQ is DOUBLE PRECISION
89!>          On entry, the value  sumsq  in the equation above.
90!>          On exit, SUMSQ is overwritten with  smsq , the basic sum of
91!>          squares from which  scl  has been factored out.
92!> \endverbatim
93!
94!  Authors:
95!  ========
96!
97!> \author Edward Anderson, Lockheed Martin
98!
99!> \par Contributors:
100!  ==================
101!>
102!> Weslley Pereira, University of Colorado Denver, USA
103!> Nick Papior, Technical University of Denmark, DK
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!>  Blue, James L. (1978)
116!>  A Portable Fortran Program to Find the Euclidean Norm of a Vector
117!>  ACM Trans Math Softw 4:15--23
118!>  https://doi.org/10.1145/355769.355771
119!>
120!> \endverbatim
121!
122!> \ingroup OTHERauxiliary
123!
124!  =====================================================================
125subroutine ZLASSQ( n, x, incx, scl, sumsq )
126   use LA_CONSTANTS, &
127      only: wp=>dp, zero=>dzero, one=>done, &
128            sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml
129   use LA_XISNAN
130!
131!  -- LAPACK auxiliary routine --
132!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
133!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
134!
135!  .. Scalar Arguments ..
136   integer :: incx, n
137   real(wp) :: scl, sumsq
138!  ..
139!  .. Array Arguments ..
140   complex(wp) :: x(*)
141!  ..
142!  .. Local Scalars ..
143   integer :: i, ix
144   logical :: notbig
145   real(wp) :: abig, amed, asml, ax, ymax, ymin
146!  ..
147!
148!  Quick return if possible
149!
150   if( LA_ISNAN(scl) .or. LA_ISNAN(sumsq) ) return
151   if( sumsq == zero ) scl = one
152   if( scl == zero ) then
153      scl = one
154      sumsq = zero
155   end if
156   if (n <= 0) then
157      return
158   end if
159!
160!  Compute the sum of squares in 3 accumulators:
161!     abig -- sums of squares scaled down to avoid overflow
162!     asml -- sums of squares scaled up to avoid underflow
163!     amed -- sums of squares that do not require scaling
164!  The thresholds and multipliers are
165!     tbig -- values bigger than this are scaled down by sbig
166!     tsml -- values smaller than this are scaled up by ssml
167!
168   notbig = .true.
169   asml = zero
170   amed = zero
171   abig = zero
172   ix = 1
173   if( incx < 0 ) ix = 1 - (n-1)*incx
174   do i = 1, n
175      ax = abs(real(x(ix)))
176      if (ax > tbig) then
177         abig = abig + (ax*sbig)**2
178         notbig = .false.
179      else if (ax < tsml) then
180         if (notbig) asml = asml + (ax*ssml)**2
181      else
182         amed = amed + ax**2
183      end if
184      ax = abs(aimag(x(ix)))
185      if (ax > tbig) then
186         abig = abig + (ax*sbig)**2
187         notbig = .false.
188      else if (ax < tsml) then
189         if (notbig) asml = asml + (ax*ssml)**2
190      else
191         amed = amed + ax**2
192      end if
193      ix = ix + incx
194   end do
195!
196!  Put the existing sum of squares into one of the accumulators
197!
198   if( sumsq > zero ) then
199      ax = scl*sqrt( sumsq )
200      if (ax > tbig) then
201         abig = abig + (ax*sbig)**2
202         notbig = .false.
203      else if (ax < tsml) then
204         if (notbig) asml = asml + (ax*ssml)**2
205      else
206         amed = amed + ax**2
207      end if
208   end if
209!
210!  Combine abig and amed or amed and asml if more than one
211!  accumulator was used.
212!
213   if (abig > zero) then
214!
215!     Combine abig and amed if abig > 0.
216!
217      if (amed > zero .or. LA_ISNAN(amed)) then
218         abig = abig + (amed*sbig)*sbig
219      end if
220      scl = one / sbig
221      sumsq = abig
222   else if (asml > zero) then
223!
224!     Combine amed and asml if asml > 0.
225!
226      if (amed > zero .or. LA_ISNAN(amed)) then
227         amed = sqrt(amed)
228         asml = sqrt(asml) / ssml
229         if (asml > amed) then
230            ymin = amed
231            ymax = asml
232         else
233            ymin = asml
234            ymax = amed
235         end if
236         scl = one
237         sumsq = ymax**2*( one + (ymin/ymax)**2 )
238      else
239         scl = one / ssml
240         sumsq = asml
241      end if
242   else
243!
244!     Otherwise all values are mid-range or zero
245!
246      scl = one
247      sumsq = amed
248   end if
249   return
250end subroutine
251