V3FIT
ddsum.f
1  function ddsum(n, din)
2  use stel_kinds
3  implicit none
4  integer :: n
5  real(rprec), intent(in) :: din(n)
6  real(rprec) :: ddsum
7 !
8 ! accurate summation by sorting and
9 ! adding in increasing order
10 !
11  logical is_all_negative, is_all_positive
12  integer istart
13  real(rprec) :: possum,negsum
14  real(rprec), parameter :: zero = 0
15 
16 
17  integer i, info
18  integer, allocatable :: iperm(:)
19  real(rprec), allocatable :: dinabs(:)
20  real(rprec) dsum
21 
22 !SIMPLE FUNCTIONALITY ddsum = SUM(din)
23 ! return
24 
25  allocate (iperm(n), dinabs(n))
26 
27 ! ------------------------
28 ! use shell sort to sort increasing abs value
29 ! ------------------------
30  dinabs = abs(din)
31  call dshell(n, dinabs, iperm )
32 
33  dsum = zero
34  do i=1,n
35  dsum = dsum + din(iperm(i))
36  enddo
37 
38  ddsum = dsum
39 
40  deallocate (iperm, dinabs)
41 
42  return
43  end function ddsum
44 
45 
46  subroutine dshell(n,dx,iperm)
47  use stel_kinds
48  integer, intent(in) :: n
49  integer, intent(out) :: iperm(n)
50  real(rprec), intent(in) :: dx(n)
51  integer :: i
52  real(rprec) :: dtemp,dxert
53 
54  integer j,k,m,maxmn
55  integer itemp
56 
57  do i=1,n
58  iperm(i) = i
59  enddo
60 
61  m = n/2
62  do while (m.gt.0)
63  maxmn = n - m
64 
65 
66  do 23004 j = 1,maxmn
67  do 23006 k = j,1,-m
68  dxert = dx(iperm(k+m))
69  if (dxert .ge. dx(iperm(k))) exit
70 
71  itemp = iperm(k+m)
72  iperm(k+m) = iperm(k)
73  iperm(k) = itemp
74 
75 
76 c dtemp = dx(k+m)
77 c dx(k+m) = dx(k)
78 c dx(k) = dtemp
79 
80 
81 23006 continue
82 C end do k
83 
84 23004 continue
85 C end do j
86 
87  m = m/2
88 
89  end do
90 C end while
91 
92  end subroutine dshell
93 
94  subroutine dshell2(n,arrin,iperm)
95  use stel_kinds
96  implicit none
97  integer, intent(in) :: n
98  integer, intent(out) :: iperm(n)
99  real(rprec), intent(in) :: arrin(n)
100  integer :: i, j, ir, m, indx
101  real(rprec) :: dtemp
102 
103  do i=1,n
104  iperm(i) = i
105  enddo
106 
107  m = n/2+1
108  ir = n
109 
110  do while (m .gt. 0)
111  if (m .gt. 1) then
112  m = m-1
113  indx = iperm(m)
114  dtemp = arrin(indx)
115  else
116  indx = iperm(ir)
117  dtemp = arrin(indx)
118  iperm(ir) = iperm(1)
119  ir = ir-1
120  if (ir .le. 1) then
121 ! if (ir .eq. 1) then
122  iperm(1) = indx
123  exit
124  end if
125  end if
126 
127  i = m
128  j = m+m
129  do while (j .le. ir)
130  if (j .lt. ir) then
131  if (arrin(iperm(j)) .lt. arrin(iperm(j+1))) j=j+1
132  endif
133  if (dtemp .lt. arrin(iperm(j))) then
134  iperm(i) = iperm(j)
135  i = j
136  j = j+j
137  else
138  j = ir+1
139  endif
140  end do
141 
142  iperm(i) = indx
143 
144  end do
145 
146  end subroutine dshell2