OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
real
Reference BLAS » Level1

Functions

subroutine crotg (a, b, c, s)
 CROTG
real(wp) function dnrm2 (n, x, incx)
 DNRM2
subroutine drotg (a, b, c, s)
 DROTG
real(wp) function dznrm2 (n, x, incx)
 DZNRM2
real function sasum (n, sx, incx)
 SASUM
subroutine saxpy (n, sa, sx, incx, sy, incy)
 SAXPY
real function scabs1 (z)
 SCABS1
real function scasum (n, cx, incx)
 SCASUM
real(wp) function scnrm2 (n, x, incx)
 SCNRM2
subroutine scopy (n, sx, incx, sy, incy)
 SCOPY
real function sdot (n, sx, incx, sy, incy)
 SDOT
real function sdsdot (n, sb, sx, incx, sy, incy)
 SDSDOT
real(wp) function snrm2 (n, x, incx)
 SNRM2
subroutine srot (n, sx, incx, sy, incy, c, s)
 SROT
subroutine srotg (a, b, c, s)
 SROTG
subroutine srotm (n, sx, incx, sy, incy, sparam)
 SROTM
subroutine srotmg (sd1, sd2, sx1, sy1, sparam)
 SROTMG
subroutine sscal (n, sa, sx, incx)
 SSCAL
subroutine sswap (n, sx, incx, sy, incy)
 SSWAP
subroutine zrotg (a, b, c, s)
 ZROTG

Detailed Description

This is the group of real LEVEL 1 BLAS routines.

Function Documentation

◆ crotg()

subroutine crotg ( complex(wp) a,
complex(wp) b,
real(wp) c,
complex(wp) s )

CROTG

Purpose:
!>
!> The computation uses the formulas
!>    |x| = sqrt( Re(x)**2 + Im(x)**2 )
!>    sgn(x) = x / |x|  if x /= 0
!>           = 1        if x  = 0
!>    c = |a| / sqrt(|a|**2 + |b|**2)
!>    s = sgn(a) * conjg(b) / sqrt(|a|**2 + |b|**2)
!> When a and b are real and r /= 0, the formulas simplify to
!>    r = sgn(a)*sqrt(|a|**2 + |b|**2)
!>    c = a / r
!>    s = b / r
!> the same as in SROTG when |a| > |b|.  When |b| >= |a|, the
!> sign of c and s will be different from those computed by SROTG
!> if the signs of a and b are not the same.
!>
!>
Parameters
[in,out]A
!>          A is COMPLEX
!>          On entry, the scalar a.
!>          On exit, the scalar r.
!>
[in]B
!>          B is COMPLEX
!>          The scalar b.
!>
[out]C
!>          C is REAL
!>          The scalar c.
!>
[out]S
!>          S is COMPLEX
!>          The scalar s.
!>
Author
Edward Anderson, Lockheed Martin
Contributors:
Weslley Pereira, University of Colorado Denver, USA
Further Details:
!>
!>  Anderson E. (2017)
!>  Algorithm 978: Safe Scaling in the Level 1 BLAS
!>  ACM Trans Math Softw 44:1--28
!>  https://doi.org/10.1145/3061665
!>
!>

Definition at line 90 of file crotg.f90.

91 integer, parameter :: wp = kind(1.e0)
92!
93! -- Reference BLAS level1 routine --
94! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
95! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
96!
97! .. Constants ..
98 real(wp), parameter :: zero = 0.0_wp
99 real(wp), parameter :: one = 1.0_wp
100 complex(wp), parameter :: czero = 0.0_wp
101! ..
102! .. Scaling constants ..
103 real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( &
104 minexponent(0._wp)-1, &
105 1-maxexponent(0._wp) &
106 )
107 real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( &
108 1-minexponent(0._wp), &
109 maxexponent(0._wp)-1 &
110 )
111 real(wp), parameter :: rtmin = sqrt( real(radix(0._wp),wp)**max( &
112 minexponent(0._wp)-1, &
113 1-maxexponent(0._wp) &
114 ) / epsilon(0._wp) )
115 real(wp), parameter :: rtmax = sqrt( real(radix(0._wp),wp)**max( &
116 1-minexponent(0._wp), &
117 maxexponent(0._wp)-1 &
118 ) * epsilon(0._wp) )
119! ..
120! .. Scalar Arguments ..
121 real(wp) :: c
122 complex(wp) :: a, b, s
123! ..
124! .. Local Scalars ..
125 real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w
126 complex(wp) :: f, fs, g, gs, r, t
127! ..
128! .. Intrinsic Functions ..
129 intrinsic :: abs, aimag, conjg, max, min, real, sqrt
130! ..
131! .. Statement Functions ..
132 real(wp) :: ABSSQ
133! ..
134! .. Statement Function definitions ..
135 abssq( t ) = real( t )**2 + aimag( t )**2
136! ..
137! .. Executable Statements ..
138!
139 f = a
140 g = b
141 if( g == czero ) then
142 c = one
143 s = czero
144 r = f
145 else if( f == czero ) then
146 c = zero
147 g1 = max( abs(real(g)), abs(aimag(g)) )
148 if( g1 > rtmin .and. g1 < rtmax ) then
149!
150! Use unscaled algorithm
151!
152 g2 = abssq( g )
153 d = sqrt( g2 )
154 s = conjg( g ) / d
155 r = d
156 else
157!
158! Use scaled algorithm
159!
160 u = min( safmax, max( safmin, g1 ) )
161 uu = one / u
162 gs = g*uu
163 g2 = abssq( gs )
164 d = sqrt( g2 )
165 s = conjg( gs ) / d
166 r = d*u
167 end if
168 else
169 f1 = max( abs(real(f)), abs(aimag(f)) )
170 g1 = max( abs(real(g)), abs(aimag(g)) )
171 if( f1 > rtmin .and. f1 < rtmax .and. &
172 g1 > rtmin .and. g1 < rtmax ) then
173!
174! Use unscaled algorithm
175!
176 f2 = abssq( f )
177 g2 = abssq( g )
178 h2 = f2 + g2
179 if( f2 > rtmin .and. h2 < rtmax ) then
180 d = sqrt( f2*h2 )
181 else
182 d = sqrt( f2 )*sqrt( h2 )
183 end if
184 p = 1 / d
185 c = f2*p
186 s = conjg( g )*( f*p )
187 r = f*( h2*p )
188 else
189!
190! Use scaled algorithm
191!
192 u = min( safmax, max( safmin, f1, g1 ) )
193 uu = one / u
194 gs = g*uu
195 g2 = abssq( gs )
196 if( f1*uu < rtmin ) then
197!
198! f is not well-scaled when scaled by g1.
199! Use a different scaling for f.
200!
201 v = min( safmax, max( safmin, f1 ) )
202 vv = one / v
203 w = v * uu
204 fs = f*vv
205 f2 = abssq( fs )
206 h2 = f2*w**2 + g2
207 else
208!
209! Otherwise use the same scaling for f and g.
210!
211 w = one
212 fs = f*uu
213 f2 = abssq( fs )
214 h2 = f2 + g2
215 end if
216 if( f2 > rtmin .and. h2 < rtmax ) then
217 d = sqrt( f2*h2 )
218 else
219 d = sqrt( f2 )*sqrt( h2 )
220 end if
221 p = 1 / d
222 c = ( f2*p )*w
223 s = conjg( gs )*( fs*p )
224 r = ( fs*( h2*p ) )*u
225 end if
226 end if
227 a = r
228 return
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ dnrm2()

real(wp) function dnrm2 ( integer n,
real(wp), dimension(*) x,
integer incx )

DNRM2

Purpose:
!>
!> DNRM2 returns the euclidean norm of a vector via the function
!> name, so that
!>
!>    DNRM2 := sqrt( x'*x )
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in]X
!>          X is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER, storage spacing between elements of X
!>          If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
!>          If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
!>          If INCX = 0, x isn't a vector so there is no need to call
!>          this subroutine.  If you call it anyway, it will count x(1)
!>          in the vector norm N times.
!>
Author
Edward Anderson, Lockheed Martin
Date
August 2016
Contributors:
Weslley Pereira, University of Colorado Denver, USA
Further Details:
!>
!>  Anderson E. (2017)
!>  Algorithm 978: Safe Scaling in the Level 1 BLAS
!>  ACM Trans Math Softw 44:1--28
!>  https://doi.org/10.1145/3061665
!>
!>  Blue, James L. (1978)
!>  A Portable Fortran Program to Find the Euclidean Norm of a Vector
!>  ACM Trans Math Softw 4:15--23
!>  https://doi.org/10.1145/355769.355771
!>
!>

Definition at line 88 of file dnrm2.f90.

89 integer, parameter :: wp = kind(1.d0)
90 real(wp) :: DNRM2
91!
92! -- Reference BLAS level1 routine (version 3.9.1) --
93! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
94! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
95! March 2021
96!
97! .. Constants ..
98 real(wp), parameter :: zero = 0.0_wp
99 real(wp), parameter :: one = 1.0_wp
100 real(wp), parameter :: maxN = huge(0.0_wp)
101! ..
102! .. Blue's scaling constants ..
103 real(wp), parameter :: tsml = real(radix(0._wp), wp)**ceiling( &
104 (minexponent(0._wp) - 1) * 0.5_wp)
105 real(wp), parameter :: tbig = real(radix(0._wp), wp)**floor( &
106 (maxexponent(0._wp) - digits(0._wp) + 1) * 0.5_wp)
107 real(wp), parameter :: ssml = real(radix(0._wp), wp)**( - floor( &
108 (minexponent(0._wp) - digits(0._wp)) * 0.5_wp))
109 real(wp), parameter :: sbig = real(radix(0._wp), wp)**( - ceiling( &
110 (maxexponent(0._wp) + digits(0._wp) - 1) * 0.5_wp))
111! ..
112! .. Scalar Arguments ..
113 integer :: incx, n
114! ..
115! .. Array Arguments ..
116 real(wp) :: x(*)
117! ..
118! .. Local Scalars ..
119 integer :: i, ix
120 logical :: notbig
121 real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
122!
123! Quick return if possible
124!
125 dnrm2 = zero
126 if( n <= 0 ) return
127!
128 scl = one
129 sumsq = zero
130!
131! Compute the sum of squares in 3 accumulators:
132! abig -- sums of squares scaled down to avoid overflow
133! asml -- sums of squares scaled up to avoid underflow
134! amed -- sums of squares that do not require scaling
135! The thresholds and multipliers are
136! tbig -- values bigger than this are scaled down by sbig
137! tsml -- values smaller than this are scaled up by ssml
138!
139 notbig = .true.
140 asml = zero
141 amed = zero
142 abig = zero
143 ix = 1
144 if( incx < 0 ) ix = 1 - (n-1)*incx
145 do i = 1, n
146 ax = abs(x(ix))
147 if (ax > tbig) then
148 abig = abig + (ax*sbig)**2
149 notbig = .false.
150 else if (ax < tsml) then
151 if (notbig) asml = asml + (ax*ssml)**2
152 else
153 amed = amed + ax**2
154 end if
155 ix = ix + incx
156 end do
157!
158! Combine abig and amed or amed and asml if more than one
159! accumulator was used.
160!
161 if (abig > zero) then
162!
163! Combine abig and amed if abig > 0.
164!
165 if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
166 abig = abig + (amed*sbig)*sbig
167 end if
168 scl = one / sbig
169 sumsq = abig
170 else if (asml > zero) then
171!
172! Combine amed and asml if asml > 0.
173!
174 if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
175 amed = sqrt(amed)
176 asml = sqrt(asml) / ssml
177 if (asml > amed) then
178 ymin = amed
179 ymax = asml
180 else
181 ymin = asml
182 ymax = amed
183 end if
184 scl = one
185 sumsq = ymax**2*( one + (ymin/ymax)**2 )
186 else
187 scl = one / ssml
188 sumsq = asml
189 end if
190 else
191!
192! Otherwise all values are mid-range
193!
194 scl = one
195 sumsq = amed
196 end if
197 dnrm2 = scl*sqrt( sumsq )
198 return
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
ymax
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
Definition law100_upd.F:272

◆ drotg()

subroutine drotg ( real(wp) a,
real(wp) b,
real(wp) c,
real(wp) s )

DROTG

Purpose:
!>
!> The computation uses the formulas
!>    sigma = sgn(a)    if |a| >  |b|
!>          = sgn(b)    if |b| >= |a|
!>    r = sigma*sqrt( a**2 + b**2 )
!>    c = 1; s = 0      if r = 0
!>    c = a/r; s = b/r  if r != 0
!> The subroutine also computes
!>    z = s    if |a| > |b|,
!>      = 1/c  if |b| >= |a| and c != 0
!>      = 1    if c = 0
!> This allows c and s to be reconstructed from z as follows:
!>    If z = 1, set c = 0, s = 1.
!>    If |z| < 1, set c = sqrt(1 - z**2) and s = z.
!>    If |z| > 1, set c = 1/z and s = sqrt( 1 - c**2).
!>
!>
Parameters
[in,out]A
!>          A is DOUBLE PRECISION
!>          On entry, the scalar a.
!>          On exit, the scalar r.
!>
[in,out]B
!>          B is DOUBLE PRECISION
!>          On entry, the scalar b.
!>          On exit, the scalar z.
!>
[out]C
!>          C is DOUBLE PRECISION
!>          The scalar c.
!>
[out]S
!>          S is DOUBLE PRECISION
!>          The scalar s.
!>
Author
Edward Anderson, Lockheed Martin
Contributors:
Weslley Pereira, University of Colorado Denver, USA
Further Details:
!>
!>  Anderson E. (2017)
!>  Algorithm 978: Safe Scaling in the Level 1 BLAS
!>  ACM Trans Math Softw 44:1--28
!>  https://doi.org/10.1145/3061665
!>
!>

Definition at line 92 of file drotg.f90.

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

◆ dznrm2()

real(wp) function dznrm2 ( integer n,
complex(wp), dimension(*) x,
integer incx )

DZNRM2

Purpose:
!>
!> DZNRM2 returns the euclidean norm of a vector via the function
!> name, so that
!>
!>    DZNRM2 := sqrt( x**H*x )
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in]X
!>          X is COMPLEX*16 array, dimension (N)
!>         complex vector with N elements
!>
[in]INCX
!>          INCX is INTEGER, storage spacing between elements of X
!>          If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
!>          If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
!>          If INCX = 0, x isn't a vector so there is no need to call
!>          this subroutine.  If you call it anyway, it will count x(1)
!>          in the vector norm N times.
!>
Author
Edward Anderson, Lockheed Martin
Date
August 2016
Contributors:
Weslley Pereira, University of Colorado Denver, USA
Further Details:
!>
!>  Anderson E. (2017)
!>  Algorithm 978: Safe Scaling in the Level 1 BLAS
!>  ACM Trans Math Softw 44:1--28
!>  https://doi.org/10.1145/3061665
!>
!>  Blue, James L. (1978)
!>  A Portable Fortran Program to Find the Euclidean Norm of a Vector
!>  ACM Trans Math Softw 4:15--23
!>  https://doi.org/10.1145/355769.355771
!>
!>

Definition at line 89 of file dznrm2.f90.

90 integer, parameter :: wp = kind(1.d0)
91 real(wp) :: DZNRM2
92!
93! -- Reference BLAS level1 routine (version 3.9.1) --
94! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
95! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
96! March 2021
97!
98! .. Constants ..
99 real(wp), parameter :: zero = 0.0_wp
100 real(wp), parameter :: one = 1.0_wp
101 real(wp), parameter :: maxN = huge(0.0_wp)
102! ..
103! .. Blue's scaling constants ..
104 real(wp), parameter :: tsml = real(radix(0._wp), wp)**ceiling( &
105 (minexponent(0._wp) - 1) * 0.5_wp)
106 real(wp), parameter :: tbig = real(radix(0._wp), wp)**floor( &
107 (maxexponent(0._wp) - digits(0._wp) + 1) * 0.5_wp)
108 real(wp), parameter :: ssml = real(radix(0._wp), wp)**( - floor( &
109 (minexponent(0._wp) - digits(0._wp)) * 0.5_wp))
110 real(wp), parameter :: sbig = real(radix(0._wp), wp)**( - ceiling( &
111 (maxexponent(0._wp) + digits(0._wp) - 1) * 0.5_wp))
112! ..
113! .. Scalar Arguments ..
114 integer :: incx, n
115! ..
116! .. Array Arguments ..
117 complex(wp) :: x(*)
118! ..
119! .. Local Scalars ..
120 integer :: i, ix
121 logical :: notbig
122 real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
123!
124! Quick return if possible
125!
126 dznrm2 = zero
127 if( n <= 0 ) return
128!
129 scl = one
130 sumsq = zero
131!
132! Compute the sum of squares in 3 accumulators:
133! abig -- sums of squares scaled down to avoid overflow
134! asml -- sums of squares scaled up to avoid underflow
135! amed -- sums of squares that do not require scaling
136! The thresholds and multipliers are
137! tbig -- values bigger than this are scaled down by sbig
138! tsml -- values smaller than this are scaled up by ssml
139!
140 notbig = .true.
141 asml = zero
142 amed = zero
143 abig = zero
144 ix = 1
145 if( incx < 0 ) ix = 1 - (n-1)*incx
146 do i = 1, n
147 ax = abs(real(x(ix)))
148 if (ax > tbig) then
149 abig = abig + (ax*sbig)**2
150 notbig = .false.
151 else if (ax < tsml) then
152 if (notbig) asml = asml + (ax*ssml)**2
153 else
154 amed = amed + ax**2
155 end if
156 ax = abs(aimag(x(ix)))
157 if (ax > tbig) then
158 abig = abig + (ax*sbig)**2
159 notbig = .false.
160 else if (ax < tsml) then
161 if (notbig) asml = asml + (ax*ssml)**2
162 else
163 amed = amed + ax**2
164 end if
165 ix = ix + incx
166 end do
167!
168! Combine abig and amed or amed and asml if more than one
169! accumulator was used.
170!
171 if (abig > zero) then
172!
173! Combine abig and amed if abig > 0.
174!
175 if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
176 abig = abig + (amed*sbig)*sbig
177 end if
178 scl = one / sbig
179 sumsq = abig
180 else if (asml > zero) then
181!
182! Combine amed and asml if asml > 0.
183!
184 if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
185 amed = sqrt(amed)
186 asml = sqrt(asml) / ssml
187 if (asml > amed) then
188 ymin = amed
189 ymax = asml
190 else
191 ymin = asml
192 ymax = amed
193 end if
194 scl = one
195 sumsq = ymax**2*( one + (ymin/ymax)**2 )
196 else
197 scl = one / ssml
198 sumsq = asml
199 end if
200 else
201!
202! Otherwise all values are mid-range
203!
204 scl = one
205 sumsq = amed
206 end if
207 dznrm2 = scl*sqrt( sumsq )
208 return
dznrm2
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition dznrm2.f90:90

◆ sasum()

real function sasum ( integer n,
real, dimension(*) sx,
integer incx )

SASUM

Purpose:
!>
!>    SASUM takes the sum of the absolute values.
!>    uses unrolled loops for increment equal to one.
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in]SX
!>          SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER
!>         storage spacing between elements of SX
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>     jack dongarra, linpack, 3/11/78.
!>     modified 3/93 to return if incx .le. 0.
!>     modified 12/3/93, array(1) declarations changed to array(*)
!>

Definition at line 71 of file sasum.f.

72*
73* -- Reference BLAS level1 routine --
74* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
75* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
76*
77* .. Scalar Arguments ..
78 INTEGER INCX,N
79* ..
80* .. Array Arguments ..
81 REAL SX(*)
82* ..
83*
84* =====================================================================
85*
86* .. Local Scalars ..
87 REAL STEMP
88 INTEGER I,M,MP1,NINCX
89* ..
90* .. Intrinsic Functions ..
91 INTRINSIC abs,mod
92* ..
93 sasum = 0.0e0
94 stemp = 0.0e0
95 IF (n.LE.0 .OR. incx.LE.0) RETURN
96 IF (incx.EQ.1) THEN
97* code for increment equal to 1
98*
99*
100* clean-up loop
101*
102 m = mod(n,6)
103 IF (m.NE.0) THEN
104 DO i = 1,m
105 stemp = stemp + abs(sx(i))
106 END DO
107 IF (n.LT.6) THEN
108 sasum = stemp
109 RETURN
110 END IF
111 END IF
112 mp1 = m + 1
113 DO i = mp1,n,6
114 stemp = stemp + abs(sx(i)) + abs(sx(i+1)) +
115 $ abs(sx(i+2)) + abs(sx(i+3)) +
116 $ abs(sx(i+4)) + abs(sx(i+5))
117 END DO
118 ELSE
119*
120* code for increment not equal to 1
121*
122 nincx = n*incx
123 DO i = 1,nincx,incx
124 stemp = stemp + abs(sx(i))
125 END DO
126 END IF
127 sasum = stemp
128 RETURN
129*
130* End of SASUM
131*
sasum
real function sasum(n, sx, incx)
SASUM
Definition sasum.f:72

◆ saxpy()

subroutine saxpy ( integer n,
real sa,
real, dimension(*) sx,
integer incx,
real, dimension(*) sy,
integer incy )

SAXPY

Purpose:
!>
!>    SAXPY constant times a vector plus a vector.
!>    uses unrolled loops for increments equal to one.
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in]SA
!>          SA is REAL
!>           On entry, SA specifies the scalar alpha.
!>
[in]SX
!>          SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER
!>         storage spacing between elements of SX
!>
[in,out]SY
!>          SY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
!>
[in]INCY
!>          INCY is INTEGER
!>         storage spacing between elements of SY
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>     jack dongarra, linpack, 3/11/78.
!>     modified 12/3/93, array(1) declarations changed to array(*)
!>

Definition at line 88 of file saxpy.f.

89*
90* -- Reference BLAS level1 routine --
91* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
92* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
93*
94* .. Scalar Arguments ..
95 REAL SA
96 INTEGER INCX,INCY,N
97* ..
98* .. Array Arguments ..
99 REAL SX(*),SY(*)
100* ..
101*
102* =====================================================================
103*
104* .. Local Scalars ..
105 INTEGER I,IX,IY,M,MP1
106* ..
107* .. Intrinsic Functions ..
108 INTRINSIC mod
109* ..
110 IF (n.LE.0) RETURN
111 IF (sa.EQ.0.0) RETURN
112 IF (incx.EQ.1 .AND. incy.EQ.1) THEN
113*
114* code for both increments equal to 1
115*
116*
117* clean-up loop
118*
119 m = mod(n,4)
120 IF (m.NE.0) THEN
121 DO i = 1,m
122 sy(i) = sy(i) + sa*sx(i)
123 END DO
124 END IF
125 IF (n.LT.4) RETURN
126 mp1 = m + 1
127 DO i = mp1,n,4
128 sy(i) = sy(i) + sa*sx(i)
129 sy(i+1) = sy(i+1) + sa*sx(i+1)
130 sy(i+2) = sy(i+2) + sa*sx(i+2)
131 sy(i+3) = sy(i+3) + sa*sx(i+3)
132 END DO
133 ELSE
134*
135* code for unequal increments or equal increments
136* not equal to 1
137*
138 ix = 1
139 iy = 1
140 IF (incx.LT.0) ix = (-n+1)*incx + 1
141 IF (incy.LT.0) iy = (-n+1)*incy + 1
142 DO i = 1,n
143 sy(iy) = sy(iy) + sa*sx(ix)
144 ix = ix + incx
145 iy = iy + incy
146 END DO
147 END IF
148 RETURN
149*
150* End of SAXPY
151*

◆ scabs1()

real function scabs1 ( complex z)

SCABS1

Purpose:
!>
!> SCABS1 computes |Re(.)| + |Im(.)| of a complex number
!>
Parameters
[in]Z
!>          Z is COMPLEX
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 45 of file scabs1.f.

46*
47* -- Reference BLAS level1 routine --
48* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
49* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
50*
51* .. Scalar Arguments ..
52 COMPLEX Z
53* ..
54*
55* =====================================================================
56*
57* .. Intrinsic Functions ..
58 INTRINSIC abs,aimag,real
59* ..
60 scabs1 = abs(real(z)) + abs(aimag(z))
61 RETURN
62*
63* End of SCABS1
64*
scabs1
real function scabs1(z)
SCABS1
Definition scabs1.f:46

◆ scasum()

real function scasum ( integer n,
complex, dimension(*) cx,
integer incx )

SCASUM

Purpose:
!>
!>    SCASUM takes the sum of the (|Re(.)| + |Im(.)|)'s of a complex vector and
!>    returns a single precision result.
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in,out]CX
!>          CX is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER
!>         storage spacing between elements of SX
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>     jack dongarra, linpack, 3/11/78.
!>     modified 3/93 to return if incx .le. 0.
!>     modified 12/3/93, array(1) declarations changed to array(*)
!>

Definition at line 71 of file scasum.f.

72*
73* -- Reference BLAS level1 routine --
74* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
75* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
76*
77* .. Scalar Arguments ..
78 INTEGER INCX,N
79* ..
80* .. Array Arguments ..
81 COMPLEX CX(*)
82* ..
83*
84* =====================================================================
85*
86* .. Local Scalars ..
87 REAL STEMP
88 INTEGER I,NINCX
89* ..
90* .. Intrinsic Functions ..
91 INTRINSIC abs,aimag,real
92* ..
93 scasum = 0.0e0
94 stemp = 0.0e0
95 IF (n.LE.0 .OR. incx.LE.0) RETURN
96 IF (incx.EQ.1) THEN
97*
98* code for increment equal to 1
99*
100 DO i = 1,n
101 stemp = stemp + abs(real(cx(i))) + abs(aimag(cx(i)))
102 END DO
103 ELSE
104*
105* code for increment not equal to 1
106*
107 nincx = n*incx
108 DO i = 1,nincx,incx
109 stemp = stemp + abs(real(cx(i))) + abs(aimag(cx(i)))
110 END DO
111 END IF
112 scasum = stemp
113 RETURN
114*
115* End of SCASUM
116*
scasum
real function scasum(n, cx, incx)
SCASUM
Definition scasum.f:72

◆ scnrm2()

real(wp) function scnrm2 ( integer n,
complex(wp), dimension(*) x,
integer incx )

SCNRM2

Purpose:
!>
!> SCNRM2 returns the euclidean norm of a vector via the function
!> name, so that
!>
!>    SCNRM2 := sqrt( x**H*x )
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in]X
!>          X is COMPLEX array, dimension (N)
!>         complex vector with N elements
!>
[in]INCX
!>          INCX is INTEGER, storage spacing between elements of X
!>          If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
!>          If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
!>          If INCX = 0, x isn't a vector so there is no need to call
!>          this subroutine.  If you call it anyway, it will count x(1)
!>          in the vector norm N times.
!>
Author
Edward Anderson, Lockheed Martin
Date
August 2016
Contributors:
Weslley Pereira, University of Colorado Denver, USA
Further Details:
!>
!>  Anderson E. (2017)
!>  Algorithm 978: Safe Scaling in the Level 1 BLAS
!>  ACM Trans Math Softw 44:1--28
!>  https://doi.org/10.1145/3061665
!>
!>  Blue, James L. (1978)
!>  A Portable Fortran Program to Find the Euclidean Norm of a Vector
!>  ACM Trans Math Softw 4:15--23
!>  https://doi.org/10.1145/355769.355771
!>
!>

Definition at line 89 of file scnrm2.f90.

90 integer, parameter :: wp = kind(1.e0)
91 real(wp) :: SCNRM2
92!
93! -- Reference BLAS level1 routine (version 3.9.1) --
94! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
95! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
96! March 2021
97!
98! .. Constants ..
99 real(wp), parameter :: zero = 0.0_wp
100 real(wp), parameter :: one = 1.0_wp
101 real(wp), parameter :: maxN = huge(0.0_wp)
102! ..
103! .. Blue's scaling constants ..
104 real(wp), parameter :: tsml = real(radix(0._wp), wp)**ceiling( &
105 (minexponent(0._wp) - 1) * 0.5_wp)
106 real(wp), parameter :: tbig = real(radix(0._wp), wp)**floor( &
107 (maxexponent(0._wp) - digits(0._wp) + 1) * 0.5_wp)
108 real(wp), parameter :: ssml = real(radix(0._wp), wp)**( - floor( &
109 (minexponent(0._wp) - digits(0._wp)) * 0.5_wp))
110 real(wp), parameter :: sbig = real(radix(0._wp), wp)**( - ceiling( &
111 (maxexponent(0._wp) + digits(0._wp) - 1) * 0.5_wp))
112! ..
113! .. Scalar Arguments ..
114 integer :: incx, n
115! ..
116! .. Array Arguments ..
117 complex(wp) :: x(*)
118! ..
119! .. Local Scalars ..
120 integer :: i, ix
121 logical :: notbig
122 real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
123!
124! Quick return if possible
125!
126 scnrm2 = zero
127 if( n <= 0 ) return
128!
129 scl = one
130 sumsq = zero
131!
132! Compute the sum of squares in 3 accumulators:
133! abig -- sums of squares scaled down to avoid overflow
134! asml -- sums of squares scaled up to avoid underflow
135! amed -- sums of squares that do not require scaling
136! The thresholds and multipliers are
137! tbig -- values bigger than this are scaled down by sbig
138! tsml -- values smaller than this are scaled up by ssml
139!
140 notbig = .true.
141 asml = zero
142 amed = zero
143 abig = zero
144 ix = 1
145 if( incx < 0 ) ix = 1 - (n-1)*incx
146 do i = 1, n
147 ax = abs(real(x(ix)))
148 if (ax > tbig) then
149 abig = abig + (ax*sbig)**2
150 notbig = .false.
151 else if (ax < tsml) then
152 if (notbig) asml = asml + (ax*ssml)**2
153 else
154 amed = amed + ax**2
155 end if
156 ax = abs(aimag(x(ix)))
157 if (ax > tbig) then
158 abig = abig + (ax*sbig)**2
159 notbig = .false.
160 else if (ax < tsml) then
161 if (notbig) asml = asml + (ax*ssml)**2
162 else
163 amed = amed + ax**2
164 end if
165 ix = ix + incx
166 end do
167!
168! Combine abig and amed or amed and asml if more than one
169! accumulator was used.
170!
171 if (abig > zero) then
172!
173! Combine abig and amed if abig > 0.
174!
175 if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
176 abig = abig + (amed*sbig)*sbig
177 end if
178 scl = one / sbig
179 sumsq = abig
180 else if (asml > zero) then
181!
182! Combine amed and asml if asml > 0.
183!
184 if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
185 amed = sqrt(amed)
186 asml = sqrt(asml) / ssml
187 if (asml > amed) then
188 ymin = amed
189 ymax = asml
190 else
191 ymin = asml
192 ymax = amed
193 end if
194 scl = one
195 sumsq = ymax**2*( one + (ymin/ymax)**2 )
196 else
197 scl = one / ssml
198 sumsq = asml
199 end if
200 else
201!
202! Otherwise all values are mid-range
203!
204 scl = one
205 sumsq = amed
206 end if
207 scnrm2 = scl*sqrt( sumsq )
208 return
scnrm2
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition scnrm2.f90:90

◆ scopy()

subroutine scopy ( integer n,
real, dimension(*) sx,
integer incx,
real, dimension(*) sy,
integer incy )

SCOPY

Purpose:
!>
!>    SCOPY copies a vector, x, to a vector, y.
!>    uses unrolled loops for increments equal to 1.
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in]SX
!>          SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER
!>         storage spacing between elements of SX
!>
[out]SY
!>          SY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
!>
[in]INCY
!>          INCY is INTEGER
!>         storage spacing between elements of SY
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>     jack dongarra, linpack, 3/11/78.
!>     modified 12/3/93, array(1) declarations changed to array(*)
!>

Definition at line 81 of file scopy.f.

82*
83* -- Reference BLAS level1 routine --
84* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
85* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
86*
87* .. Scalar Arguments ..
88 INTEGER INCX,INCY,N
89* ..
90* .. Array Arguments ..
91 REAL SX(*),SY(*)
92* ..
93*
94* =====================================================================
95*
96* .. Local Scalars ..
97 INTEGER I,IX,IY,M,MP1
98* ..
99* .. Intrinsic Functions ..
100 INTRINSIC mod
101* ..
102 IF (n.LE.0) RETURN
103 IF (incx.EQ.1 .AND. incy.EQ.1) THEN
104*
105* code for both increments equal to 1
106*
107*
108* clean-up loop
109*
110 m = mod(n,7)
111 IF (m.NE.0) THEN
112 DO i = 1,m
113 sy(i) = sx(i)
114 END DO
115 IF (n.LT.7) RETURN
116 END IF
117 mp1 = m + 1
118 DO i = mp1,n,7
119 sy(i) = sx(i)
120 sy(i+1) = sx(i+1)
121 sy(i+2) = sx(i+2)
122 sy(i+3) = sx(i+3)
123 sy(i+4) = sx(i+4)
124 sy(i+5) = sx(i+5)
125 sy(i+6) = sx(i+6)
126 END DO
127 ELSE
128*
129* code for unequal increments or equal increments
130* not equal to 1
131*
132 ix = 1
133 iy = 1
134 IF (incx.LT.0) ix = (-n+1)*incx + 1
135 IF (incy.LT.0) iy = (-n+1)*incy + 1
136 DO i = 1,n
137 sy(iy) = sx(ix)
138 ix = ix + incx
139 iy = iy + incy
140 END DO
141 END IF
142 RETURN
143*
144* End of SCOPY
145*

◆ sdot()

real function sdot ( integer n,
real, dimension(*) sx,
integer incx,
real, dimension(*) sy,
integer incy )

SDOT

Purpose:
!>
!>    SDOT forms the dot product of two vectors.
!>    uses unrolled loops for increments equal to one.
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in]SX
!>          SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER
!>         storage spacing between elements of SX
!>
[in]SY
!>          SY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
!>
[in]INCY
!>          INCY is INTEGER
!>         storage spacing between elements of SY
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>     jack dongarra, linpack, 3/11/78.
!>     modified 12/3/93, array(1) declarations changed to array(*)
!>

Definition at line 81 of file sdot.f.

82*
83* -- Reference BLAS level1 routine --
84* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
85* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
86*
87* .. Scalar Arguments ..
88 INTEGER INCX,INCY,N
89* ..
90* .. Array Arguments ..
91 REAL SX(*),SY(*)
92* ..
93*
94* =====================================================================
95*
96* .. Local Scalars ..
97 REAL STEMP
98 INTEGER I,IX,IY,M,MP1
99* ..
100* .. Intrinsic Functions ..
101 INTRINSIC mod
102* ..
103 stemp = 0.0e0
104 sdot = 0.0e0
105 IF (n.LE.0) RETURN
106 IF (incx.EQ.1 .AND. incy.EQ.1) THEN
107*
108* code for both increments equal to 1
109*
110*
111* clean-up loop
112*
113 m = mod(n,5)
114 IF (m.NE.0) THEN
115 DO i = 1,m
116 stemp = stemp + sx(i)*sy(i)
117 END DO
118 IF (n.LT.5) THEN
119 sdot=stemp
120 RETURN
121 END IF
122 END IF
123 mp1 = m + 1
124 DO i = mp1,n,5
125 stemp = stemp + sx(i)*sy(i) + sx(i+1)*sy(i+1) +
126 $ sx(i+2)*sy(i+2) + sx(i+3)*sy(i+3) + sx(i+4)*sy(i+4)
127 END DO
128 ELSE
129*
130* code for unequal increments or equal increments
131* not equal to 1
132*
133 ix = 1
134 iy = 1
135 IF (incx.LT.0) ix = (-n+1)*incx + 1
136 IF (incy.LT.0) iy = (-n+1)*incy + 1
137 DO i = 1,n
138 stemp = stemp + sx(ix)*sy(iy)
139 ix = ix + incx
140 iy = iy + incy
141 END DO
142 END IF
143 sdot = stemp
144 RETURN
145*
146* End of SDOT
147*
sdot
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82

◆ sdsdot()

real function sdsdot ( integer n,
real sb,
real, dimension(*) sx,
integer incx,
real, dimension(*) sy,
integer incy )

SDSDOT

Purpose:
!>
!>   Compute the inner product of two vectors with extended
!>   precision accumulation.
!>
!>   Returns S.P. result with dot product accumulated in D.P.
!>   SDSDOT = SB + sum for I = 0 to N-1 of SX(LX+I*INCX)*SY(LY+I*INCY),
!>   where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
!>   defined in a similar way using INCY.
!>
Parameters
[in]N
!>          N is INTEGER
!>          number of elements in input vector(s)
!>
[in]SB
!>          SB is REAL
!>          single precision scalar to be added to inner product
!>
[in]SX
!>          SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>          single precision vector with N elements
!>
[in]INCX
!>          INCX is INTEGER
!>          storage spacing between elements of SX
!>
[in]SY
!>          SY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>          single precision vector with N elements
!>
[in]INCY
!>          INCY is INTEGER
!>          storage spacing between elements of SY
!>
Author
Lawson, C. L., (JPL), Hanson, R. J., (SNLA),
Kincaid, D. R., (U. of Texas), Krogh, F. T., (JPL)
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>    REFERENCES
!>
!>    C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
!>    Krogh, Basic linear algebra subprograms for Fortran
!>    usage, Algorithm No. 539, Transactions on Mathematical
!>    Software 5, 3 (September 1979), pp. 308-323.
!>
!>    REVISION HISTORY  (YYMMDD)
!>
!>    791001  DATE WRITTEN
!>    890531  Changed all specific intrinsics to generic.  (WRB)
!>    890831  Modified array declarations.  (WRB)
!>    890831  REVISION DATE from Version 3.2
!>    891214  Prologue converted to Version 4.0 format.  (BAB)
!>    920310  Corrected definition of LX in DESCRIPTION.  (WRB)
!>    920501  Reformatted the REFERENCES section.  (WRB)
!>    070118  Reformat to LAPACK coding style
!>

Definition at line 112 of file sdsdot.f.

113*
114* -- Reference BLAS level1 routine --
115* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
116* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
117*
118* .. Scalar Arguments ..
119 REAL SB
120 INTEGER INCX,INCY,N
121* ..
122* .. Array Arguments ..
123 REAL SX(*),SY(*)
124* .. Local Scalars ..
125 DOUBLE PRECISION DSDOT
126 INTEGER I,KX,KY,NS
127* ..
128* .. Intrinsic Functions ..
129 INTRINSIC dble
130* ..
131 dsdot = sb
132 IF (n.LE.0) THEN
133 sdsdot = dsdot
134 RETURN
135 END IF
136 IF (incx.EQ.incy .AND. incx.GT.0) THEN
137*
138* Code for equal and positive increments.
139*
140 ns = n*incx
141 DO i = 1,ns,incx
142 dsdot = dsdot + dble(sx(i))*dble(sy(i))
143 END DO
144 ELSE
145*
146* Code for unequal or nonpositive increments.
147*
148 kx = 1
149 ky = 1
150 IF (incx.LT.0) kx = 1 + (1-n)*incx
151 IF (incy.LT.0) ky = 1 + (1-n)*incy
152 DO i = 1,n
153 dsdot = dsdot + dble(sx(kx))*dble(sy(ky))
154 kx = kx + incx
155 ky = ky + incy
156 END DO
157 END IF
158 sdsdot = dsdot
159 RETURN
160*
161* End of SDSDOT
162*
dsdot
double precision function dsdot(n, sx, incx, sy, incy)
DSDOT
Definition dsdot.f:119
sdsdot
real function sdsdot(n, sb, sx, incx, sy, incy)
SDSDOT
Definition sdsdot.f:113

◆ snrm2()

real(wp) function snrm2 ( integer n,
real(wp), dimension(*) x,
integer incx )

SNRM2

Purpose:
!>
!> SNRM2 returns the euclidean norm of a vector via the function
!> name, so that
!>
!>    SNRM2 := sqrt( x'*x ).
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in]X
!>          X is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER, storage spacing between elements of X
!>          If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
!>          If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
!>          If INCX = 0, x isn't a vector so there is no need to call
!>          this subroutine.  If you call it anyway, it will count x(1)
!>          in the vector norm N times.
!>
Author
Edward Anderson, Lockheed Martin
Date
August 2016
Contributors:
Weslley Pereira, University of Colorado Denver, USA
Further Details:
!>
!>  Anderson E. (2017)
!>  Algorithm 978: Safe Scaling in the Level 1 BLAS
!>  ACM Trans Math Softw 44:1--28
!>  https://doi.org/10.1145/3061665
!>
!>  Blue, James L. (1978)
!>  A Portable Fortran Program to Find the Euclidean Norm of a Vector
!>  ACM Trans Math Softw 4:15--23
!>  https://doi.org/10.1145/355769.355771
!>
!>

Definition at line 88 of file snrm2.f90.

89 integer, parameter :: wp = kind(1.e0)
90 real(wp) :: SNRM2
91!
92! -- Reference BLAS level1 routine (version 3.9.1) --
93! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
94! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
95! March 2021
96!
97! .. Constants ..
98 real(wp), parameter :: zero = 0.0_wp
99 real(wp), parameter :: one = 1.0_wp
100 real(wp), parameter :: maxN = huge(0.0_wp)
101! ..
102! .. Blue's scaling constants ..
103 real(wp), parameter :: tsml = real(radix(0._wp), wp)**ceiling( &
104 (minexponent(0._wp) - 1) * 0.5_wp)
105 real(wp), parameter :: tbig = real(radix(0._wp), wp)**floor( &
106 (maxexponent(0._wp) - digits(0._wp) + 1) * 0.5_wp)
107 real(wp), parameter :: ssml = real(radix(0._wp), wp)**( - floor( &
108 (minexponent(0._wp) - digits(0._wp)) * 0.5_wp))
109 real(wp), parameter :: sbig = real(radix(0._wp), wp)**( - ceiling( &
110 (maxexponent(0._wp) + digits(0._wp) - 1) * 0.5_wp))
111! ..
112! .. Scalar Arguments ..
113 integer :: incx, n
114! ..
115! .. Array Arguments ..
116 real(wp) :: x(*)
117! ..
118! .. Local Scalars ..
119 integer :: i, ix
120 logical :: notbig
121 real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
122!
123! Quick return if possible
124!
125 snrm2 = zero
126 if( n <= 0 ) return
127!
128 scl = one
129 sumsq = zero
130!
131! Compute the sum of squares in 3 accumulators:
132! abig -- sums of squares scaled down to avoid overflow
133! asml -- sums of squares scaled up to avoid underflow
134! amed -- sums of squares that do not require scaling
135! The thresholds and multipliers are
136! tbig -- values bigger than this are scaled down by sbig
137! tsml -- values smaller than this are scaled up by ssml
138!
139 notbig = .true.
140 asml = zero
141 amed = zero
142 abig = zero
143 ix = 1
144 if( incx < 0 ) ix = 1 - (n-1)*incx
145 do i = 1, n
146 ax = abs(x(ix))
147 if (ax > tbig) then
148 abig = abig + (ax*sbig)**2
149 notbig = .false.
150 else if (ax < tsml) then
151 if (notbig) asml = asml + (ax*ssml)**2
152 else
153 amed = amed + ax**2
154 end if
155 ix = ix + incx
156 end do
157!
158! Combine abig and amed or amed and asml if more than one
159! accumulator was used.
160!
161 if (abig > zero) then
162!
163! Combine abig and amed if abig > 0.
164!
165 if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
166 abig = abig + (amed*sbig)*sbig
167 end if
168 scl = one / sbig
169 sumsq = abig
170 else if (asml > zero) then
171!
172! Combine amed and asml if asml > 0.
173!
174 if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
175 amed = sqrt(amed)
176 asml = sqrt(asml) / ssml
177 if (asml > amed) then
178 ymin = amed
179 ymax = asml
180 else
181 ymin = asml
182 ymax = amed
183 end if
184 scl = one
185 sumsq = ymax**2*( one + (ymin/ymax)**2 )
186 else
187 scl = one / ssml
188 sumsq = asml
189 end if
190 else
191!
192! Otherwise all values are mid-range
193!
194 scl = one
195 sumsq = amed
196 end if
197 snrm2 = scl*sqrt( sumsq )
198 return
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89

◆ srot()

subroutine srot ( integer n,
real, dimension(*) sx,
integer incx,
real, dimension(*) sy,
integer incy,
real c,
real s )

SROT

Purpose:
!>
!>    applies a plane rotation.
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in,out]SX
!>          SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER
!>         storage spacing between elements of SX
!>
[in,out]SY
!>          SY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
!>
[in]INCY
!>          INCY is INTEGER
!>         storage spacing between elements of SY
!>
[in]C
!>          C is REAL
!>
[in]S
!>          S is REAL
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>     jack dongarra, linpack, 3/11/78.
!>     modified 12/3/93, array(1) declarations changed to array(*)
!>

Definition at line 91 of file srot.f.

92*
93* -- Reference BLAS level1 routine --
94* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
95* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
96*
97* .. Scalar Arguments ..
98 REAL C,S
99 INTEGER INCX,INCY,N
100* ..
101* .. Array Arguments ..
102 REAL SX(*),SY(*)
103* ..
104*
105* =====================================================================
106*
107* .. Local Scalars ..
108 REAL STEMP
109 INTEGER I,IX,IY
110* ..
111 IF (n.LE.0) RETURN
112 IF (incx.EQ.1 .AND. incy.EQ.1) THEN
113*
114* code for both increments equal to 1
115*
116 DO i = 1,n
117 stemp = c*sx(i) + s*sy(i)
118 sy(i) = c*sy(i) - s*sx(i)
119 sx(i) = stemp
120 END DO
121 ELSE
122*
123* code for unequal increments or equal increments not equal
124* to 1
125*
126 ix = 1
127 iy = 1
128 IF (incx.LT.0) ix = (-n+1)*incx + 1
129 IF (incy.LT.0) iy = (-n+1)*incy + 1
130 DO i = 1,n
131 stemp = c*sx(ix) + s*sy(iy)
132 sy(iy) = c*sy(iy) - s*sx(ix)
133 sx(ix) = stemp
134 ix = ix + incx
135 iy = iy + incy
136 END DO
137 END IF
138 RETURN
139*
140* End of SROT
141*

◆ srotg()

subroutine srotg ( real(wp) a,
real(wp) b,
real(wp) c,
real(wp) s )

SROTG

Purpose:
!>
!> The computation uses the formulas
!>    sigma = sgn(a)    if |a| >  |b|
!>          = sgn(b)    if |b| >= |a|
!>    r = sigma*sqrt( a**2 + b**2 )
!>    c = 1; s = 0      if r = 0
!>    c = a/r; s = b/r  if r != 0
!> The subroutine also computes
!>    z = s    if |a| > |b|,
!>      = 1/c  if |b| >= |a| and c != 0
!>      = 1    if c = 0
!> This allows c and s to be reconstructed from z as follows:
!>    If z = 1, set c = 0, s = 1.
!>    If |z| < 1, set c = sqrt(1 - z**2) and s = z.
!>    If |z| > 1, set c = 1/z and s = sqrt( 1 - c**2).
!>
!>
Parameters
[in,out]A
!>          A is REAL
!>          On entry, the scalar a.
!>          On exit, the scalar r.
!>
[in,out]B
!>          B is REAL
!>          On entry, the scalar b.
!>          On exit, the scalar z.
!>
[out]C
!>          C is REAL
!>          The scalar c.
!>
[out]S
!>          S is REAL
!>          The scalar s.
!>
Author
Edward Anderson, Lockheed Martin
Contributors:
Weslley Pereira, University of Colorado Denver, USA
Further Details:
!>
!>  Anderson E. (2017)
!>  Algorithm 978: Safe Scaling in the Level 1 BLAS
!>  ACM Trans Math Softw 44:1--28
!>  https://doi.org/10.1145/3061665
!>
!>

Definition at line 92 of file srotg.f90.

93 integer, parameter :: wp = kind(1.e0)
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

◆ srotm()

subroutine srotm ( integer n,
real, dimension(*) sx,
integer incx,
real, dimension(*) sy,
integer incy,
real, dimension(5) sparam )

SROTM

Purpose:
!>
!>    APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX
!>
!>    (SX**T) , WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF SX ARE IN
!>    (SX**T)
!>
!>    SX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE
!>    LX = (-INCX)*N, AND SIMILARLY FOR SY USING USING LY AND INCY.
!>    WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
!>
!>    SFLAG=-1.E0     SFLAG=0.E0        SFLAG=1.E0     SFLAG=-2.E0
!>
!>      (SH11  SH12)    (1.E0  SH12)    (SH11  1.E0)    (1.E0  0.E0)
!>    H=(          )    (          )    (          )    (          )
!>      (SH21  SH22),   (SH21  1.E0),   (-1.E0 SH22),   (0.E0  1.E0).
!>    SEE  SROTMG FOR A DESCRIPTION OF DATA STORAGE IN SPARAM.
!>
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in,out]SX
!>          SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER
!>         storage spacing between elements of SX
!>
[in,out]SY
!>          SY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
!>
[in]INCY
!>          INCY is INTEGER
!>         storage spacing between elements of SY
!>
[in]SPARAM
!>          SPARAM is REAL array, dimension (5)
!>     SPARAM(1)=SFLAG
!>     SPARAM(2)=SH11
!>     SPARAM(3)=SH21
!>     SPARAM(4)=SH12
!>     SPARAM(5)=SH22
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 96 of file srotm.f.

97*
98* -- Reference BLAS level1 routine --
99* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
100* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
101*
102* .. Scalar Arguments ..
103 INTEGER INCX,INCY,N
104* ..
105* .. Array Arguments ..
106 REAL SPARAM(5),SX(*),SY(*)
107* ..
108*
109* =====================================================================
110*
111* .. Local Scalars ..
112 REAL SFLAG,SH11,SH12,SH21,SH22,TWO,W,Z,ZERO
113 INTEGER I,KX,KY,NSTEPS
114* ..
115* .. Data statements ..
116 DATA zero,two/0.e0,2.e0/
117* ..
118*
119 sflag = sparam(1)
120 IF (n.LE.0 .OR. (sflag+two.EQ.zero)) RETURN
121 IF (incx.EQ.incy.AND.incx.GT.0) THEN
122*
123 nsteps = n*incx
124 IF (sflag.LT.zero) THEN
125 sh11 = sparam(2)
126 sh12 = sparam(4)
127 sh21 = sparam(3)
128 sh22 = sparam(5)
129 DO i = 1,nsteps,incx
130 w = sx(i)
131 z = sy(i)
132 sx(i) = w*sh11 + z*sh12
133 sy(i) = w*sh21 + z*sh22
134 END DO
135 ELSE IF (sflag.EQ.zero) THEN
136 sh12 = sparam(4)
137 sh21 = sparam(3)
138 DO i = 1,nsteps,incx
139 w = sx(i)
140 z = sy(i)
141 sx(i) = w + z*sh12
142 sy(i) = w*sh21 + z
143 END DO
144 ELSE
145 sh11 = sparam(2)
146 sh22 = sparam(5)
147 DO i = 1,nsteps,incx
148 w = sx(i)
149 z = sy(i)
150 sx(i) = w*sh11 + z
151 sy(i) = -w + sh22*z
152 END DO
153 END IF
154 ELSE
155 kx = 1
156 ky = 1
157 IF (incx.LT.0) kx = 1 + (1-n)*incx
158 IF (incy.LT.0) ky = 1 + (1-n)*incy
159*
160 IF (sflag.LT.zero) THEN
161 sh11 = sparam(2)
162 sh12 = sparam(4)
163 sh21 = sparam(3)
164 sh22 = sparam(5)
165 DO i = 1,n
166 w = sx(kx)
167 z = sy(ky)
168 sx(kx) = w*sh11 + z*sh12
169 sy(ky) = w*sh21 + z*sh22
170 kx = kx + incx
171 ky = ky + incy
172 END DO
173 ELSE IF (sflag.EQ.zero) THEN
174 sh12 = sparam(4)
175 sh21 = sparam(3)
176 DO i = 1,n
177 w = sx(kx)
178 z = sy(ky)
179 sx(kx) = w + z*sh12
180 sy(ky) = w*sh21 + z
181 kx = kx + incx
182 ky = ky + incy
183 END DO
184 ELSE
185 sh11 = sparam(2)
186 sh22 = sparam(5)
187 DO i = 1,n
188 w = sx(kx)
189 z = sy(ky)
190 sx(kx) = w*sh11 + z
191 sy(ky) = -w + sh22*z
192 kx = kx + incx
193 ky = ky + incy
194 END DO
195 END IF
196 END IF
197 RETURN
198*
199* End of SROTM
200*

◆ srotmg()

subroutine srotmg ( real sd1,
real sd2,
real sx1,
real sy1,
real, dimension(5) sparam )

SROTMG

Purpose:
!>
!>    CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
!>    THE SECOND COMPONENT OF THE 2-VECTOR  (SQRT(SD1)*SX1,SQRT(SD2)*>    SY2)**T.
!>    WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
!>
!>    SFLAG=-1.E0     SFLAG=0.E0        SFLAG=1.E0     SFLAG=-2.E0
!>
!>      (SH11  SH12)    (1.E0  SH12)    (SH11  1.E0)    (1.E0  0.E0)
!>    H=(          )    (          )    (          )    (          )
!>      (SH21  SH22),   (SH21  1.E0),   (-1.E0 SH22),   (0.E0  1.E0).
!>    LOCATIONS 2-4 OF SPARAM CONTAIN SH11,SH21,SH12, AND SH22
!>    RESPECTIVELY. (VALUES OF 1.E0, -1.E0, OR 0.E0 IMPLIED BY THE
!>    VALUE OF SPARAM(1) ARE NOT STORED IN SPARAM.)
!>
!>    THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE
!>    INEXACT.  THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE
!>    OF SD1 AND SD2.  ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
!>
!>
Parameters
[in,out]SD1
!>          SD1 is REAL
!>
[in,out]SD2
!>          SD2 is REAL
!>
[in,out]SX1
!>          SX1 is REAL
!>
[in]SY1
!>          SY1 is REAL
!>
[out]SPARAM
!>          SPARAM is REAL array, dimension (5)
!>     SPARAM(1)=SFLAG
!>     SPARAM(2)=SH11
!>     SPARAM(3)=SH21
!>     SPARAM(4)=SH12
!>     SPARAM(5)=SH22
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 89 of file srotmg.f.

90*
91* -- Reference BLAS level1 routine --
92* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
93* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
94*
95* .. Scalar Arguments ..
96 REAL SD1,SD2,SX1,SY1
97* ..
98* .. Array Arguments ..
99 REAL SPARAM(5)
100* ..
101*
102* =====================================================================
103*
104* .. Local Scalars ..
105 REAL GAM,GAMSQ,ONE,RGAMSQ,SFLAG,SH11,SH12,SH21,SH22,SP1,SP2,SQ1,
106 $ SQ2,STEMP,SU,TWO,ZERO
107* ..
108* .. Intrinsic Functions ..
109 INTRINSIC abs
110* ..
111* .. Data statements ..
112*
113 DATA zero,one,two/0.e0,1.e0,2.e0/
114 DATA gam,gamsq,rgamsq/4096.e0,1.67772e7,5.96046e-8/
115* ..
116
117 IF (sd1.LT.zero) THEN
118* GO ZERO-H-D-AND-SX1..
119 sflag = -one
120 sh11 = zero
121 sh12 = zero
122 sh21 = zero
123 sh22 = zero
124*
125 sd1 = zero
126 sd2 = zero
127 sx1 = zero
128 ELSE
129* CASE-SD1-NONNEGATIVE
130 sp2 = sd2*sy1
131 IF (sp2.EQ.zero) THEN
132 sflag = -two
133 sparam(1) = sflag
134 RETURN
135 END IF
136* REGULAR-CASE..
137 sp1 = sd1*sx1
138 sq2 = sp2*sy1
139 sq1 = sp1*sx1
140*
141 IF (abs(sq1).GT.abs(sq2)) THEN
142 sh21 = -sy1/sx1
143 sh12 = sp2/sp1
144*
145 su = one - sh12*sh21
146*
147 IF (su.GT.zero) THEN
148 sflag = zero
149 sd1 = sd1/su
150 sd2 = sd2/su
151 sx1 = sx1*su
152 ELSE
153* This code path if here for safety. We do not expect this
154* condition to ever hold except in edge cases with rounding
155* errors. See DOI: 10.1145/355841.355847
156 sflag = -one
157 sh11 = zero
158 sh12 = zero
159 sh21 = zero
160 sh22 = zero
161*
162 sd1 = zero
163 sd2 = zero
164 sx1 = zero
165 END IF
166 ELSE
167
168 IF (sq2.LT.zero) THEN
169* GO ZERO-H-D-AND-SX1..
170 sflag = -one
171 sh11 = zero
172 sh12 = zero
173 sh21 = zero
174 sh22 = zero
175*
176 sd1 = zero
177 sd2 = zero
178 sx1 = zero
179 ELSE
180 sflag = one
181 sh11 = sp1/sp2
182 sh22 = sx1/sy1
183 su = one + sh11*sh22
184 stemp = sd2/su
185 sd2 = sd1/su
186 sd1 = stemp
187 sx1 = sy1*su
188 END IF
189 END IF
190
191* PROCEDURE..SCALE-CHECK
192 IF (sd1.NE.zero) THEN
193 DO WHILE ((sd1.LE.rgamsq) .OR. (sd1.GE.gamsq))
194 IF (sflag.EQ.zero) THEN
195 sh11 = one
196 sh22 = one
197 sflag = -one
198 ELSE
199 sh21 = -one
200 sh12 = one
201 sflag = -one
202 END IF
203 IF (sd1.LE.rgamsq) THEN
204 sd1 = sd1*gam**2
205 sx1 = sx1/gam
206 sh11 = sh11/gam
207 sh12 = sh12/gam
208 ELSE
209 sd1 = sd1/gam**2
210 sx1 = sx1*gam
211 sh11 = sh11*gam
212 sh12 = sh12*gam
213 END IF
214 ENDDO
215 END IF
216
217 IF (sd2.NE.zero) THEN
218 DO WHILE ( (abs(sd2).LE.rgamsq) .OR. (abs(sd2).GE.gamsq) )
219 IF (sflag.EQ.zero) THEN
220 sh11 = one
221 sh22 = one
222 sflag = -one
223 ELSE
224 sh21 = -one
225 sh12 = one
226 sflag = -one
227 END IF
228 IF (abs(sd2).LE.rgamsq) THEN
229 sd2 = sd2*gam**2
230 sh21 = sh21/gam
231 sh22 = sh22/gam
232 ELSE
233 sd2 = sd2/gam**2
234 sh21 = sh21*gam
235 sh22 = sh22*gam
236 END IF
237 END DO
238 END IF
239
240 END IF
241
242 IF (sflag.LT.zero) THEN
243 sparam(2) = sh11
244 sparam(3) = sh21
245 sparam(4) = sh12
246 sparam(5) = sh22
247 ELSE IF (sflag.EQ.zero) THEN
248 sparam(3) = sh21
249 sparam(4) = sh12
250 ELSE
251 sparam(2) = sh11
252 sparam(5) = sh22
253 END IF
254
255 sparam(1) = sflag
256 RETURN
257*
258* End of SROTMG
259*

◆ sscal()

subroutine sscal ( integer n,
real sa,
real, dimension(*) sx,
integer incx )

SSCAL

Purpose:
!>
!>    SSCAL scales a vector by a constant.
!>    uses unrolled loops for increment equal to 1.
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in]SA
!>          SA is REAL
!>           On entry, SA specifies the scalar alpha.
!>
[in,out]SX
!>          SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER
!>         storage spacing between elements of SX
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>     jack dongarra, linpack, 3/11/78.
!>     modified 3/93 to return if incx .le. 0.
!>     modified 12/3/93, array(1) declarations changed to array(*)
!>

Definition at line 78 of file sscal.f.

79*
80* -- Reference BLAS level1 routine --
81* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
82* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
83*
84* .. Scalar Arguments ..
85 REAL SA
86 INTEGER INCX,N
87* ..
88* .. Array Arguments ..
89 REAL SX(*)
90* ..
91*
92* =====================================================================
93*
94* .. Local Scalars ..
95 INTEGER I,M,MP1,NINCX
96* ..
97* .. Intrinsic Functions ..
98 INTRINSIC mod
99* ..
100 IF (n.LE.0 .OR. incx.LE.0) RETURN
101 IF (incx.EQ.1) THEN
102*
103* code for increment equal to 1
104*
105*
106* clean-up loop
107*
108 m = mod(n,5)
109 IF (m.NE.0) THEN
110 DO i = 1,m
111 sx(i) = sa*sx(i)
112 END DO
113 IF (n.LT.5) RETURN
114 END IF
115 mp1 = m + 1
116 DO i = mp1,n,5
117 sx(i) = sa*sx(i)
118 sx(i+1) = sa*sx(i+1)
119 sx(i+2) = sa*sx(i+2)
120 sx(i+3) = sa*sx(i+3)
121 sx(i+4) = sa*sx(i+4)
122 END DO
123 ELSE
124*
125* code for increment not equal to 1
126*
127 nincx = n*incx
128 DO i = 1,nincx,incx
129 sx(i) = sa*sx(i)
130 END DO
131 END IF
132 RETURN
133*
134* End of SSCAL
135*

◆ sswap()

subroutine sswap ( integer n,
real, dimension(*) sx,
integer incx,
real, dimension(*) sy,
integer incy )

SSWAP

Purpose:
!>
!>    SSWAP interchanges two vectors.
!>    uses unrolled loops for increments equal to 1.
!>
Parameters
[in]N
!>          N is INTEGER
!>         number of elements in input vector(s)
!>
[in,out]SX
!>          SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!>
[in]INCX
!>          INCX is INTEGER
!>         storage spacing between elements of SX
!>
[in,out]SY
!>          SY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
!>
[in]INCY
!>          INCY is INTEGER
!>         storage spacing between elements of SY
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>     jack dongarra, linpack, 3/11/78.
!>     modified 12/3/93, array(1) declarations changed to array(*)
!>

Definition at line 81 of file sswap.f.

82*
83* -- Reference BLAS level1 routine --
84* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
85* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
86*
87* .. Scalar Arguments ..
88 INTEGER INCX,INCY,N
89* ..
90* .. Array Arguments ..
91 REAL SX(*),SY(*)
92* ..
93*
94* =====================================================================
95*
96* .. Local Scalars ..
97 REAL STEMP
98 INTEGER I,IX,IY,M,MP1
99* ..
100* .. Intrinsic Functions ..
101 INTRINSIC mod
102* ..
103 IF (n.LE.0) RETURN
104 IF (incx.EQ.1 .AND. incy.EQ.1) THEN
105*
106* code for both increments equal to 1
107*
108*
109* clean-up loop
110*
111 m = mod(n,3)
112 IF (m.NE.0) THEN
113 DO i = 1,m
114 stemp = sx(i)
115 sx(i) = sy(i)
116 sy(i) = stemp
117 END DO
118 IF (n.LT.3) RETURN
119 END IF
120 mp1 = m + 1
121 DO i = mp1,n,3
122 stemp = sx(i)
123 sx(i) = sy(i)
124 sy(i) = stemp
125 stemp = sx(i+1)
126 sx(i+1) = sy(i+1)
127 sy(i+1) = stemp
128 stemp = sx(i+2)
129 sx(i+2) = sy(i+2)
130 sy(i+2) = stemp
131 END DO
132 ELSE
133*
134* code for unequal increments or equal increments not equal
135* to 1
136*
137 ix = 1
138 iy = 1
139 IF (incx.LT.0) ix = (-n+1)*incx + 1
140 IF (incy.LT.0) iy = (-n+1)*incy + 1
141 DO i = 1,n
142 stemp = sx(ix)
143 sx(ix) = sy(iy)
144 sy(iy) = stemp
145 ix = ix + incx
146 iy = iy + incy
147 END DO
148 END IF
149 RETURN
150*
151* End of SSWAP
152*

◆ zrotg()

subroutine zrotg ( complex(wp) a,
complex(wp) b,
real(wp) c,
complex(wp) s )

ZROTG

Purpose:
!>
!> The computation uses the formulas
!>    |x| = sqrt( Re(x)**2 + Im(x)**2 )
!>    sgn(x) = x / |x|  if x /= 0
!>           = 1        if x  = 0
!>    c = |a| / sqrt(|a|**2 + |b|**2)
!>    s = sgn(a) * conjg(b) / sqrt(|a|**2 + |b|**2)
!> When a and b are real and r /= 0, the formulas simplify to
!>    r = sgn(a)*sqrt(|a|**2 + |b|**2)
!>    c = a / r
!>    s = b / r
!> the same as in DROTG when |a| > |b|.  When |b| >= |a|, the
!> sign of c and s will be different from those computed by DROTG
!> if the signs of a and b are not the same.
!>
!>
Parameters
[in,out]A
!>          A is DOUBLE COMPLEX
!>          On entry, the scalar a.
!>          On exit, the scalar r.
!>
[in]B
!>          B is DOUBLE COMPLEX
!>          The scalar b.
!>
[out]C
!>          C is DOUBLE PRECISION
!>          The scalar c.
!>
[out]S
!>          S is DOUBLE COMPLEX
!>          The scalar s.
!>
Author
Edward Anderson, Lockheed Martin
Contributors:
Weslley Pereira, University of Colorado Denver, USA
Further Details:
!>
!>  Anderson E. (2017)
!>  Algorithm 978: Safe Scaling in the Level 1 BLAS
!>  ACM Trans Math Softw 44:1--28
!>  https://doi.org/10.1145/3061665
!>
!>

Definition at line 90 of file zrotg.f90.

91 integer, parameter :: wp = kind(1.d0)
92!
93! -- Reference BLAS level1 routine --
94! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
95! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
96!
97! .. Constants ..
98 real(wp), parameter :: zero = 0.0_wp
99 real(wp), parameter :: one = 1.0_wp
100 complex(wp), parameter :: czero = 0.0_wp
101! ..
102! .. Scaling constants ..
103 real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( &
104 minexponent(0._wp)-1, &
105 1-maxexponent(0._wp) &
106 )
107 real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( &
108 1-minexponent(0._wp), &
109 maxexponent(0._wp)-1 &
110 )
111 real(wp), parameter :: rtmin = sqrt( real(radix(0._wp),wp)**max( &
112 minexponent(0._wp)-1, &
113 1-maxexponent(0._wp) &
114 ) / epsilon(0._wp) )
115 real(wp), parameter :: rtmax = sqrt( real(radix(0._wp),wp)**max( &
116 1-minexponent(0._wp), &
117 maxexponent(0._wp)-1 &
118 ) * epsilon(0._wp) )
119! ..
120! .. Scalar Arguments ..
121 real(wp) :: c
122 complex(wp) :: a, b, s
123! ..
124! .. Local Scalars ..
125 real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w
126 complex(wp) :: f, fs, g, gs, r, t
127! ..
128! .. Intrinsic Functions ..
129 intrinsic :: abs, aimag, conjg, max, min, real, sqrt
130! ..
131! .. Statement Functions ..
132 real(wp) :: ABSSQ
133! ..
134! .. Statement Function definitions ..
135 abssq( t ) = real( t )**2 + aimag( t )**2
136! ..
137! .. Executable Statements ..
138!
139 f = a
140 g = b
141 if( g == czero ) then
142 c = one
143 s = czero
144 r = f
145 else if( f == czero ) then
146 c = zero
147 g1 = max( abs(real(g)), abs(aimag(g)) )
148 if( g1 > rtmin .and. g1 < rtmax ) then
149!
150! Use unscaled algorithm
151!
152 g2 = abssq( g )
153 d = sqrt( g2 )
154 s = conjg( g ) / d
155 r = d
156 else
157!
158! Use scaled algorithm
159!
160 u = min( safmax, max( safmin, g1 ) )
161 uu = one / u
162 gs = g*uu
163 g2 = abssq( gs )
164 d = sqrt( g2 )
165 s = conjg( gs ) / d
166 r = d*u
167 end if
168 else
169 f1 = max( abs(real(f)), abs(aimag(f)) )
170 g1 = max( abs(real(g)), abs(aimag(g)) )
171 if( f1 > rtmin .and. f1 < rtmax .and. &
172 g1 > rtmin .and. g1 < rtmax ) then
173!
174! Use unscaled algorithm
175!
176 f2 = abssq( f )
177 g2 = abssq( g )
178 h2 = f2 + g2
179 if( f2 > rtmin .and. h2 < rtmax ) then
180 d = sqrt( f2*h2 )
181 else
182 d = sqrt( f2 )*sqrt( h2 )
183 end if
184 p = 1 / d
185 c = f2*p
186 s = conjg( g )*( f*p )
187 r = f*( h2*p )
188 else
189!
190! Use scaled algorithm
191!
192 u = min( safmax, max( safmin, f1, g1 ) )
193 uu = one / u
194 gs = g*uu
195 g2 = abssq( gs )
196 if( f1*uu < rtmin ) then
197!
198! f is not well-scaled when scaled by g1.
199! Use a different scaling for f.
200!
201 v = min( safmax, max( safmin, f1 ) )
202 vv = one / v
203 w = v * uu
204 fs = f*vv
205 f2 = abssq( fs )
206 h2 = f2*w**2 + g2
207 else
208!
209! Otherwise use the same scaling for f and g.
210!
211 w = one
212 fs = f*uu
213 f2 = abssq( fs )
214 h2 = f2 + g2
215 end if
216 if( f2 > rtmin .and. h2 < rtmax ) then
217 d = sqrt( f2*h2 )
218 else
219 d = sqrt( f2 )*sqrt( h2 )
220 end if
221 p = 1 / d
222 c = ( f2*p )*w
223 s = conjg( gs )*( fs*p )
224 r = ( fs*( h2*p ) )*u
225 end if
226 end if
227 a = r
228 return