#include "wwm_functions.h"
!-----------------------------------------------------------------------
! subroutine from blas1.f90
!-----------------------------------------------------------------------
function dnrm2 ( n, dx, incx)
!Mathieu here zero was not defined and implicit none not specified ...
! major bug!
use datapool, only : rkind, THR, zero, one
real(rkind) :: dnrm2
integer, intent(in) :: incx
integer next
real(rkind) dx(*), cutlo, cuthi, hitest, sumd, xmax
integer n, nn, i, j
!
! euclidean norm of the n-vector stored in dx() with storage
! increment incx .
! if n .le. 0 return with result = 0.
! if n .ge. 1 then incx must be .ge. 1
!
! c.l.lawson, 1978 jan 08
!
! four phase method using two built-in constants that are
! hopefully applicable to all machines.
! cutlo = maximum of dsqrt(u/eps) over all known machines.
! cuthi = minimum of dsqrt(v) over all known machines.
! where
! eps = smallest no. such that eps + 1. .gt. 1.
! u = smallest positive no. (underflow limit)
! v = largest no. (overflow limit)
!
! brief outline of algorithm..
!
! phase 1 scans zero components.
! move to phase 2 when a component is nonzero and .le. cutlo
! move to phase 3 when a component is .gt. cutlo
! move to phase 4 when a component is .ge. cuthi/m
! where m = n for x() real and m = 2*n for complex.
!
! values for cutlo and cuthi..
! from the environmental parameters listed in the imsl converter
! document the limiting values are as follows..
! cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are
! univac and dec at 2**(-103)
! thus cutlo = 2**(-51) = 4.44089e-16
! cuthi, s.p. v = 2**127 for univac, honeywell, and dec.
! thus cuthi = 2**(63.5) = 1.30438e19
! cutlo, d.p. u/eps = 2**(-67) for honeywell and dec.
! thus cutlo = 2**(-33.5) = 8.23181d-11
! cuthi, d.p. same as s.p. cuthi = 1.30438d19
! data cutlo, cuthi / 8.232d-11, 1.304d19 /
! data cutlo, cuthi / 4.441e-16, 1.304e19 /
data cutlo, cuthi / 8.232d-11, 1.304d19 /
!
if(n .gt. 0) go to 10
dnrm2 = zero
go to 300
!
10 assign 30 to next
sumd = zero
nn = n * incx
! begin main loop
i = 1
20 go to next,(30, 50, 70, 110)
30 if( MyABS(dx(i)) .gt. cutlo) go to 85
assign 50 to next
xmax = zero
!
! phase 1. sumd is zero
!
50 if( abs(dx(i)) .lt. THR) go to 200
if( MyABS(dx(i)) .gt. cutlo) go to 85
!
! prepare for phase 2.
assign 70 to next
go to 105
!
! prepare for phase 4.
!
100 i = j
assign 110 to next
sumd = (sumd / dx(i)) / dx(i)
105 xmax = MyABS(dx(i))
go to 115
!
! phase 2. sumd is small.
! scale to avoid destructive underflow.
!
70 if( MyABS(dx(i)) .gt. cutlo ) go to 75
!
! common code for phases 2 and 4.
! in phase 4 sumd is large. scale to avoid overflow.
!
110 if( MyABS(dx(i)) .le. xmax ) go to 115
sumd = one + sumd * (xmax / dx(i))**2
xmax = MyABS(dx(i))
go to 200
!
115 sumd = sumd + (dx(i)/xmax)**2
go to 200
!
!
! prepare for phase 3.
!
75 sumd = (sumd * xmax) * xmax
!
!
! for real or d.p. set hitest = cuthi/n
! for complex set hitest = cuthi/(2*n)
!
85 hitest = cuthi/float( n )
!
! phase 3. sumd is mid-range. no scaling.
!
do 95 j =i,nn,incx
if(MyABS(dx(j)) .ge. hitest) go to 100
95 sumd = sumd + dx(j)**2
dnrm2 = MySQRT( sumd )
go to 300
!
200 continue
i = i + incx
if ( i .le. nn ) go to 20
!
! end of main loop.
!
! compute square root and adjust for scaling.
!
dnrm2 = xmax * MySQRT(sumd)
300 continue
return
end
!-------------------------------------------------------------------------
#ifndef SCHISM
real(rkind) function ddot(n,dx,incx,dy,incy)
use datapool, only : rkind, ZERO
!
! forms the dot product of two vectors.
! uses unrolled loops for increments equal to one.
! jack dongarra, linpack, 3/11/78.
!
real(rkind) dx(*),dy(*),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
!
ddot = ZERO
dtemp = ZERO
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
!
! code for unequal increments or equal increments
! not equal to 1
!
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
10 continue
ddot = dtemp
return
!
! code for both increments equal to 1
!
!
! clean-up loop
!
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dx(i)*dy(i)
30 continue
if( n .lt. 5 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,5
dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + &
& dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
50 continue
60 ddot = dtemp
return
end
#endif
!----------------------------------------------------------------------
subroutine daxpy(n,da,dx,incx,dy,incy)
use datapool, only : rkind, THR
!
! constant times a vector plus a vector.
! uses unrolled loops for increments equal to one.
! jack dongarra, linpack, 3/11/78.
!
real(rkind) dx(1),dy(1),da
integer i,incx,incy,ix,iy,m,mp1,n
!
if(n.le.0)return
if (abs(da) .lt. THR) return
if(incx.eq.1.and.incy.eq.1)go to 20
!
! code for unequal increments or equal increments
! not equal to 1
!
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dy(iy) = dy(iy) + da*dx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
!
! code for both increments equal to 1
!
!
! clean-up loop
!
20 m = mod(n,4)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dy(i) = dy(i) + da*dx(i)
30 continue
if( n .lt. 4 ) return
40 mp1 = m + 1
do 50 i = mp1,n,4
dy(i) = dy(i) + da*dx(i)
dy(i + 1) = dy(i + 1) + da*dx(i + 1)
dy(i + 2) = dy(i + 2) + da*dx(i + 2)
dy(i + 3) = dy(i + 3) + da*dx(i + 3)
50 continue
return
end