← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/calelm/wsort.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/calelm/wsort.f
 
1
      subroutine wsort(countr,counti,n,index,test)
 
2
c
 
3
c!purpose
 
4
c     wsort sort double precision array according to rule specified by
 
5
c      external test
 
6
c     maintaining an index array 
 
7
c
 
8
c!calling sequence
 
9
c     subroutine wsort(countr,counti,n,index,test)
 
10
c     integer n,index(n)
 
11
c     double precision count(n)
 
12
c     integer test
 
13
c     external test
 
14
c
 
15
c     count(r,i)   : array to be sorted
 
16
c     n       :size of count and index
 
17
c     index   : array containing on return index of sorted array
 
18
c     test    : external integer function which define formal order for
 
19
c               records 
 
20
c               test(r1,i1,r2,i2)
 
21
c               where
 
22
c               r1,i1 are real and imag part of first complex number
 
23
c               r2,i2 are real and imag part of second complex number
 
24
c               returns
 
25
c                1 :if 1 is greater than 2
 
26
c               -1 :if 1 is less than 2
 
27
c                0 :if 1 is equal to 2
 
28
c
 
29
c!method
 
30
c     quick sort metjod is used
 
31
c!restriction
 
32
c     n must be less than 2**(50/2) ! due to lengh of work space mark
 
33
c!
 
34
c     Copyright INRIA
 
35
      dimension mark(50),index(n)
 
36
      double precision countr(n),counti(n),avr,avi,xr,xi
 
37
      integer test
 
38
      external test
 
39
c
 
40
c  set index array to original order .
 
41
      do 10 i=1,n
 
42
      index(i)=i
 
43
   10 continue
 
44
c  check that a trivial case has not been entered .
 
45
      if(n.eq.1)goto 200
 
46
      if(n.ge.1)go to 30
 
47
      goto 200
 
48
c  'm' is the length of segment which is short enough to enter
 
49
c  the final sorting routine. it may be easily changed.
 
50
   30 m=12
 
51
c  set up initial values.
 
52
      la=2
 
53
      is=1
 
54
      if=n
 
55
      do 190 mloop=1,n
 
56
c  if segment is short enough sort with final sorting routine .
 
57
      ifka=if-is
 
58
      if((ifka+1).gt.m)goto 70
 
59
c********* final sorting ***
 
60
c  ( a simple bubble sort )
 
61
      is1=is+1
 
62
      do 60 j=is1,if
 
63
      i=j
 
64
   40 it=test(countr(i-1),counti(i-1),countr(i),counti(i))
 
65
      if(it.eq.1) goto 60
 
66
      if(it.eq.-1) goto 50
 
67
      if(index(i-1).lt.index(i))goto 60
 
68
   50 avr=countr(i-1)
 
69
      avi=counti(i-1)
 
70
      countr(i-1)=countr(i)
 
71
      counti(i-1)=counti(i)
 
72
      countr(i)=avr
 
73
      counti(i)=avi
 
74
      int=index(i-1)
 
75
      index(i-1)=index(i)
 
76
      index(i)=int
 
77
      i=i-1
 
78
      if(i.gt.is)goto  40
 
79
   60 continue
 
80
      la=la-2
 
81
      goto 170
 
82
c             *******  quicksort  ********
 
83
c  select the number in the central position in the segment as
 
84
c  the test number.replace it with the number from the segment's
 
85
c  highest address.
 
86
   70 iy=(is+if)/2
 
87
      xr=countr(iy)
 
88
      xi=counti(iy)
 
89
      intest=index(iy)
 
90
      countr(iy)=countr(if)
 
91
      counti(iy)=counti(if)
 
92
      index(iy)=index(if)
 
93
c  the markers 'i' and 'ifk' are used for the beginning and end
 
94
c  of the section not so far tested against the present value
 
95
c  of x .
 
96
      k=1
 
97
      ifk=if
 
98
c  we alternate between the outer loop that increases i and the
 
99
c  inner loop that reduces ifk, moving numbers and indices as
 
100
c  necessary, until they meet .
 
101
      do 110 i=is,if
 
102
         it=test(xr,xi,countr(i),counti(i))
 
103
         if(it.lt.0) goto 110
 
104
         if(it.gt.0) goto 80
 
105
         if(intest.gt.index(i))goto 110
 
106
 80      if(i.ge.ifk)goto 120
 
107
         countr(ifk)=countr(i)
 
108
         counti(ifk)=counti(i)
 
109
         index(ifk)=index(i)
 
110
         k1=k
 
111
         do 100 k=k1,ifka
 
112
            ifk=if-k
 
113
            it=test(countr(ifk),counti(ifk),xr,xi)
 
114
            if(it.lt.0) goto 100
 
115
            if(it.gt.0) goto 90
 
116
            if(intest.le.index(ifk))goto 100
 
117
 90         if(i.ge.ifk)goto 130
 
118
            countr(i)=countr(ifk)
 
119
            counti(i)=counti(ifk)
 
120
            index(i)=index(ifk)
 
121
            go to 110
 
122
 100     continue
 
123
         goto 120
 
124
 110  continue
 
125
c  return the test number to the position marked by the marker
 
126
c  which did not move last. it divides the initial segment into
 
127
c  2 parts. any element in the first part is less than or equal
 
128
c  to any element in the second part, and they may now be sorted
 
129
c  independently .
 
130
  120 countr(ifk)=xr
 
131
      counti(ifk)=xi
 
132
      index(ifk)=intest
 
133
      ip=ifk
 
134
      goto 140
 
135
  130 countr(i)=xr
 
136
      counti(i)=xi
 
137
      index(i)=intest
 
138
      ip=i
 
139
c  store the longer subdivision in workspace.
 
140
  140 if((ip-is).gt.(if-ip))goto 150
 
141
      mark(la)=if
 
142
      mark(la-1)=ip+1
 
143
      if=ip-1
 
144
      goto 160
 
145
  150 mark(la)=ip-1
 
146
      mark(la-1)=is
 
147
      is=ip+1
 
148
c  find the length of the shorter subdivision.
 
149
  160 lngth=if-is
 
150
      if(lngth.le.0)goto 180
 
151
c  if it contains more than one element supply it with workspace .
 
152
      la=la+2
 
153
      goto 190
 
154
  170 if(la.le.0)goto 200
 
155
c  obtain the address of the shortest segment awaiting quicksort
 
156
  180 if=mark(la)
 
157
      is=mark(la-1)
 
158
  190 continue
 
159
  200 return
 
160
      end
 
161
      integer function rptest(r1,i1,r2,i2)
 
162
      double precision r1,i1,r2,i2
 
163
      if(r1.gt.r2) then
 
164
         rptest=1
 
165
      elseif(r1.lt.r2) then
 
166
         rptest=-1
 
167
      else
 
168
         rptest=0
 
169
      endif
 
170
      end
 
171
      integer function modtest(r1,i1,r2,i2)
 
172
      double precision r1,i1,r2,i2,a1,a2
 
173
      a1=r1**2+i1**2
 
174
      a2=r2**2+i2**2
 
175
      if(a1.gt.a2) then
 
176
         modtest=1
 
177
      elseif(a1.lt.a2) then
 
178
         modtest=-1
 
179
      else
 
180
         modtest=0
 
181
      endif
 
182
      end