ADMB Documentation  -a65f1c97
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ivsort2.cpp
Go to the documentation of this file.
1 /*
2  * $Id$
3  * Author: David Fournier
4  *
5  * Copyright (c) 2009 ADMB Foundation
6  */
7 #include <fvar.hpp>
8 
9 namespace admb_deprecated {
10 
20 ivector sort(const ivector& v, int NSTACK)
21 {
22  const int M=7;
23  const int FM=7875;
24  const int FA=211;
25  const int FC=1663;
26 
27  int n=v.size();
28  ivector arr(v.indexmin(),v.indexmax());
29  arr=v;
30  arr.shift(1);
31  int l=1,jstack=0,j,ir,iq,i;
32  ivector istack(1,NSTACK+1);
33  long int fx=0L;
34  int a;
35 
36  ir=n;
37  for (;;)
38  {
39  if (ir-l < M)
40  {
41  for (j=l+1;j<=ir;j++)
42  {
43  a=arr[j];
44  for (i=j-1;i>0 && arr[i]>a;i--) arr[i+1]=arr[i];
45  arr[i+1]=a;
46  }
47  if (jstack == 0)
48  {
49  arr.shift(v.indexmin());
50  return arr;
51  }
52  ir=istack[jstack--];
53  l=istack[jstack--];
54  }
55  else
56  {
57  i=l;
58  j=ir;
59  fx=(fx*FA+FC) % FM;
60  iq=l+((ir-l+1)*fx)/FM;
61  a=arr[iq];
62  arr[iq]=arr[l];
63  for (;;)
64  {
65  while (j > 0 && a < arr[j]) j--;
66  if (j <= i)
67  {
68  arr[i]=a;
69  break;
70  }
71  arr[i++]=arr[j];
72  while (i <= n && a > arr[i] ) i++;
73  if (j <= i)
74  {
75  arr[(i=j)]=a;
76  break;
77  }
78  arr[j--]=arr[i];
79  }
80  if (ir-i >= i-l)
81  {
82  istack[++jstack]=i+1;
83  istack[++jstack]=ir;
84  ir=i-1;
85  } else
86  {
87  istack[++jstack]=l;
88  istack[++jstack]=i-1;
89  l=i+1;
90  }
91  if (jstack > NSTACK)
92  {
93  cerr << "Need to increase the stack in sort(const ivector&)\n";
94  ad_exit(1);
95  }
96  }
97  }
98 }
99 
110 ivector sort(const ivector& _v, const ivector& _index, int NSTACK)
111 {
112  const int M=7;
113  const int FM=7875;
114  const int FA=211;
115  const int FC=1663;
116  int iii;
117  ivector& index = (ivector&) _index;
118  ivector& v = (ivector&) _v;
119 
120  if (v.size() != index.size())
121  {
122  cerr << " Incompatible array sizes in vector v and ivector index\n"
123  << " in ivector sort(const ivector& v, const ivector& index)\n";
124  ad_exit(1);
125  }
126  int n=v.size();
127  int One=1;
128  int nne=v.indexmin();
129  index.fill_seqadd(nne,One);
130  ivector arr(v.indexmin(),v.indexmax());
131  arr=v;
132  arr.shift(1);
133  index.shift(1);
134  int l=1,jstack=0,j,ir,iq,i;
135  ivector istack(1,NSTACK+1);
136  long int fx=0L;
137  int a;
138 
139  ir=n;
140  int i1;
141  int tmp;
142  for (;;)
143  {
144  if (ir-l < M)
145  {
146  for (j=l+1;j<=ir;j++)
147  {
148  a=arr[j];
149  iii=index[j];
150  for (i=j-1;i>0;i--)
151  {
152  tmp=arr[i];
153  i1=i+1;
154  if (tmp <=a) break;
155  arr[i1]=tmp;
156  index[i1]=index[i];
157  }
158  arr[i+1]=a;
159  index[i+1]=iii;
160  }
161  if (jstack == 0)
162  {
163  arr.shift(v.indexmin());
164  index.shift(v.indexmin());
165  return arr;
166  }
167  ir=istack[jstack--];
168  l=istack[jstack--];
169  }
170  else
171  {
172  i=l;
173  j=ir;
174  fx=(fx*FA+FC) % FM;
175  iq=l+((ir-l+1)*fx)/FM;
176  a=arr[iq];
177  iii=index[iq];
178  arr[iq]=arr[l];
179  index[iq]=index[l];
180  for (;;)
181  {
182  while (j > 0 && a < arr[j]) j--;
183  if (j <= i)
184  {
185  arr[i]=a;
186  index[i]=iii;
187  break;
188  }
189  arr[i]=arr[j];
190  index[i++]=index[j];
191  while (i <= n && a > arr[i]) i++;
192  if (j <= i)
193  {
194  arr[(i=j)]=a;
195  index[i]=iii;
196  break;
197  }
198  arr[j]=arr[i];
199  index[j--]=index[i];
200  }
201  if (ir-i >= i-l)
202  {
203  istack[++jstack]=i+1;
204  istack[++jstack]=ir;
205  ir=i-1;
206  } else
207  {
208  istack[++jstack]=l;
209  istack[++jstack]=i-1;
210  l=i+1;
211  }
212  if (jstack > NSTACK)
213  {
214  cerr << "Need to increase the stack in sort(const ivector&)\n";
215  ad_exit(1);
216  }
217  }
218  }
219 }
220 
221 }
dvector sort(const dvector &v, int NSTACK)
Quicksort.
Definition: dvsort2.cpp:20
unsigned int size() const
Definition: ivector.h:109
exitptr ad_exit
Definition: gradstrc.cpp:53
void fill_seqadd(int, int)
Fills ivector elements with values starting from base and incremented by offset.
Definition: cranfill.cpp:73
Array of integers(int) with indexes from index_min to indexmax.
Definition: ivector.h:50
Author: David Fournier Copyright (c) 2008-2012 Regents of the University of California.
int indexmin() const
Definition: ivector.h:99
#define M
Definition: rngen.cpp:57
int indexmax() const
Definition: ivector.h:104
ivector & shift(int min)
Author: David Fournier Copyright (c) 2008-2012 Regents of the University of California.
Definition: ivec.cpp:17