1 |
C ddot.f |
2 |
C is freely available from netlib. It is not subject to any GNU License |
3 |
C set by the authors of the ASCEND math programming system. |
4 |
C $Date: 1996/04/30 18:14:07 $ $Revision: 1.1.1.1 $ |
5 |
C |
6 |
double precision function ddot(n,dx,incx,dy,incy) |
7 |
c |
8 |
c forms the dot product of two vectors. |
9 |
c uses unrolled loops for increments equal to one. |
10 |
c jack dongarra, linpack, 3/11/78. |
11 |
c |
12 |
double precision dx(*),dy(*),dtemp |
13 |
integer i,incx,incy,ix,iy,m,mp1,n |
14 |
c |
15 |
ddot = 0.0d0 |
16 |
dtemp = 0.0d0 |
17 |
if(n.le.0)return |
18 |
if(incx.eq.1.and.incy.eq.1)go to 20 |
19 |
c |
20 |
c code for unequal increments or equal increments |
21 |
c not equal to 1 |
22 |
c |
23 |
ix = 1 |
24 |
iy = 1 |
25 |
if(incx.lt.0)ix = (-n+1)*incx + 1 |
26 |
if(incy.lt.0)iy = (-n+1)*incy + 1 |
27 |
do 10 i = 1,n |
28 |
dtemp = dtemp + dx(ix)*dy(iy) |
29 |
ix = ix + incx |
30 |
iy = iy + incy |
31 |
10 continue |
32 |
ddot = dtemp |
33 |
return |
34 |
c |
35 |
c code for both increments equal to 1 |
36 |
c |
37 |
c |
38 |
c clean-up loop |
39 |
c |
40 |
20 m = mod(n,5) |
41 |
if( m .eq. 0 ) go to 40 |
42 |
do 30 i = 1,m |
43 |
dtemp = dtemp + dx(i)*dy(i) |
44 |
30 continue |
45 |
if( n .lt. 5 ) go to 60 |
46 |
40 mp1 = m + 1 |
47 |
do 50 i = mp1,n,5 |
48 |
dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + |
49 |
* dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) |
50 |
50 continue |
51 |
60 ddot = dtemp |
52 |
return |
53 |
end |
54 |
|
55 |
|