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