source: Sophya/trunk/SophyaLib/NTools/nbtrixx.cc@ 4051

Last change on this file since 4051 was 3685, checked in by cmv, 16 years ago
  • changement du code C-Fortran like addressing par du code C++ purement C-like addressing (NR C++ version)

cmv, 27/11/2009
File size: 6.5 KB
RevLine 
[3678]1#include "machdefs.h"
2#include <stdio.h>
3#include <stdlib.h>
4
5#include "nbtrixx.h"
6
7namespace SOPHYA {
8
[3685]9typedef int_4 Int; // 32 bit integer
10
[3678]11//-------------------------------------------------------------
12template <class T>
[3685]13void TabSort(uint_4 n, T* arr)
[3678]14{
[3685]15 static const Int M=7, NSTACK=64;
16 Int i,ir,j,k,jstack=-1,l=0;
17 T a;
18 Int* istack = new Int[NSTACK];
19 ir=n-1;
[3678]20 for (;;) {
21 if (ir-l < M) {
22 for (j=l+1;j<=ir;j++) {
23 a=arr[j];
24 for (i=j-1;i>=l;i--) {
25 if (arr[i] <= a) break;
26 arr[i+1]=arr[i];
27 }
28 arr[i+1]=a;
29 }
[3685]30 if (jstack < 0) break;
[3678]31 ir=istack[jstack--];
32 l=istack[jstack--];
33 } else {
34 k=(l+ir) >> 1;
[3685]35 _NR_SWAP(arr[k],arr[l+1]);
[3678]36 if (arr[l] > arr[ir]) {
[3685]37 _NR_SWAP(arr[l],arr[ir]);
[3678]38 }
39 if (arr[l+1] > arr[ir]) {
[3685]40 _NR_SWAP(arr[l+1],arr[ir]);
[3678]41 }
42 if (arr[l] > arr[l+1]) {
[3685]43 _NR_SWAP(arr[l],arr[l+1]);
[3678]44 }
45 i=l+1;
46 j=ir;
47 a=arr[l+1];
48 for (;;) {
49 do i++; while (arr[i] < a);
50 do j--; while (arr[j] > a);
51 if (j < i) break;
[3685]52 _NR_SWAP(arr[i],arr[j]);
[3678]53 }
54 arr[l+1]=arr[j];
55 arr[j]=a;
56 jstack += 2;
[3685]57 if (jstack >= NSTACK) throw("NSTACK too small in sort.");
[3678]58 if (ir-i+1 >= j-l) {
59 istack[jstack]=ir;
60 istack[jstack-1]=i;
61 ir=j-1;
62 } else {
63 istack[jstack]=j-1;
64 istack[jstack-1]=l;
65 l=i;
66 }
67 }
68 }
[3685]69 delete [] istack;
[3678]70}
71
72//-------------------------------------------------------------
[3685]73template <class T, class U>
74void TabSort2(uint_4 n, T *arr, U *brr)
[3678]75{
[3685]76 const Int M=7,NSTACK=64;
77 Int i,ir,j,k,jstack=-1,l=0;
78 T a;
79 U b;
80 Int* istack = new Int[NSTACK];
81 ir=n-1;
[3678]82 for (;;) {
83 if (ir-l < M) {
84 for (j=l+1;j<=ir;j++) {
85 a=arr[j];
86 b=brr[j];
87 for (i=j-1;i>=l;i--) {
88 if (arr[i] <= a) break;
89 arr[i+1]=arr[i];
90 brr[i+1]=brr[i];
91 }
92 arr[i+1]=a;
93 brr[i+1]=b;
94 }
[3685]95 if (jstack < 0) break;
96 ir=istack[jstack--];
97 l=istack[jstack--];
[3678]98 } else {
99 k=(l+ir) >> 1;
[3685]100 _NR_SWAP(arr[k],arr[l+1]);
101 _NR_SWAP(brr[k],brr[l+1]);
[3678]102 if (arr[l] > arr[ir]) {
[3685]103 _NR_SWAP(arr[l],arr[ir]);
104 _NR_SWAP(brr[l],brr[ir]);
[3678]105 }
106 if (arr[l+1] > arr[ir]) {
[3685]107 _NR_SWAP(arr[l+1],arr[ir]);
108 _NR_SWAP(brr[l+1],brr[ir]);
[3678]109 }
110 if (arr[l] > arr[l+1]) {
[3685]111 _NR_SWAP(arr[l],arr[l+1]);
112 _NR_SWAP(brr[l],brr[l+1]);
[3678]113 }
114 i=l+1;
115 j=ir;
116 a=arr[l+1];
117 b=brr[l+1];
118 for (;;) {
119 do i++; while (arr[i] < a);
120 do j--; while (arr[j] > a);
121 if (j < i) break;
[3685]122 _NR_SWAP(arr[i],arr[j]);
123 _NR_SWAP(brr[i],brr[j]);
[3678]124 }
125 arr[l+1]=arr[j];
126 arr[j]=a;
127 brr[l+1]=brr[j];
128 brr[j]=b;
129 jstack += 2;
[3685]130 if (jstack >= NSTACK) throw("NSTACK too small in sort2.");
[3678]131 if (ir-i+1 >= j-l) {
132 istack[jstack]=ir;
133 istack[jstack-1]=i;
134 ir=j-1;
135 } else {
136 istack[jstack]=j-1;
137 istack[jstack-1]=l;
138 l=i;
139 }
140 }
141 }
[3685]142 delete [] istack;
[3678]143}
144
[3685]145
[3678]146//-------------------------------------------------------------
147template <class T>
[3685]148void TabSortInd(uint_4 n, T *arr, uint_4 *indx)
[3678]149{
[3685]150 const Int M=7,NSTACK=64;
151 Int i,indxt,ir,j,k,jstack=-1,l=0;
[3678]152 T a;
[3685]153 Int* istack = new Int[NSTACK];
154 ir=n-1;
155 for (j=0;j<(Int)n;j++) indx[j]=j;
[3678]156 for (;;) {
157 if (ir-l < M) {
158 for (j=l+1;j<=ir;j++) {
159 indxt=indx[j];
160 a=arr[indxt];
161 for (i=j-1;i>=l;i--) {
162 if (arr[indx[i]] <= a) break;
163 indx[i+1]=indx[i];
164 }
165 indx[i+1]=indxt;
166 }
[3685]167 if (jstack < 0) break;
[3678]168 ir=istack[jstack--];
169 l=istack[jstack--];
170 } else {
171 k=(l+ir) >> 1;
[3685]172 _NR_SWAP(indx[k],indx[l+1]);
[3678]173 if (arr[indx[l]] > arr[indx[ir]]) {
[3685]174 _NR_SWAP(indx[l],indx[ir]);
[3678]175 }
176 if (arr[indx[l+1]] > arr[indx[ir]]) {
[3685]177 _NR_SWAP(indx[l+1],indx[ir]);
[3678]178 }
179 if (arr[indx[l]] > arr[indx[l+1]]) {
[3685]180 _NR_SWAP(indx[l],indx[l+1]);
[3678]181 }
182 i=l+1;
183 j=ir;
184 indxt=indx[l+1];
185 a=arr[indxt];
186 for (;;) {
187 do i++; while (arr[indx[i]] < a);
188 do j--; while (arr[indx[j]] > a);
189 if (j < i) break;
[3685]190 _NR_SWAP(indx[i],indx[j]);
[3678]191 }
192 indx[l+1]=indx[j];
193 indx[j]=indxt;
194 jstack += 2;
[3685]195 if (jstack >= NSTACK) throw("NSTACK too small in index.");
[3678]196 if (ir-i+1 >= j-l) {
197 istack[jstack]=ir;
198 istack[jstack-1]=i;
199 ir=j-1;
200 } else {
201 istack[jstack]=j-1;
202 istack[jstack-1]=l;
203 l=i;
204 }
205 }
206 }
[3685]207 delete [] istack;
[3678]208}
209
[3685]210
[3678]211//-------------------------------------------------------------
212#ifdef __CXX_PRAGMA_TEMPLATES__
213#pragma define_template TabSort<int_2>
214#pragma define_template TabSort<uint_2>
215#pragma define_template TabSort<int_4>
216#pragma define_template TabSort<uint_4>
217#pragma define_template TabSort<int_8>
218#pragma define_template TabSort<uint_8>
219#pragma define_template TabSort<r_4>
220#pragma define_template TabSort<r_8>
221
[3685]222#pragma define_template TabSort2<int_2,int_2>
223#pragma define_template TabSort2<uint_2,uint_2>
224#pragma define_template TabSort2<int_4,int_4>
225#pragma define_template TabSort2<uint_4,uint_4>
226#pragma define_template TabSort2<int_8,int_8>
227#pragma define_template TabSort2<uint_8,uint_8>
228#pragma define_template TabSort2<r_4,r_4>
229#pragma define_template TabSort2<r_8,r_8>
230#pragma define_template TabSort2<r_4,int_4>
231#pragma define_template TabSort2<r_8,int_4>
[3678]232
233#pragma define_template TabSortInd<int_2>
234#pragma define_template TabSortInd<uint_2>
235#pragma define_template TabSortInd<int_4>
236#pragma define_template TabSortInd<uint_4>
237#pragma define_template TabSortInd<int_8>
238#pragma define_template TabSortInd<uint_8>
239#pragma define_template TabSortInd<r_4>
240#pragma define_template TabSortInd<r_8>
241#endif
242
243#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
244template void TabSort(uint_4, int_2*);
245template void TabSort(uint_4, uint_2*);
246template void TabSort(uint_4, int_4*);
247template void TabSort(uint_4, uint_4*);
248template void TabSort(uint_4, int_8*);
249template void TabSort(uint_4, uint_8*);
250template void TabSort(uint_4, r_4*);
251template void TabSort(uint_4, r_8*);
252
253template void TabSort2(uint_4, int_2*, int_2*);
254template void TabSort2(uint_4, uint_2*, uint_2*);
255template void TabSort2(uint_4, int_4*, int_4*);
256template void TabSort2(uint_4, uint_4*, uint_4*);
257template void TabSort2(uint_4, int_8*, int_8*);
258template void TabSort2(uint_4, uint_8*, uint_8*);
259template void TabSort2(uint_4, r_4*, r_4*);
260template void TabSort2(uint_4, r_8*, r_8*);
[3685]261template void TabSort2(uint_4, r_4*, int_4*);
262template void TabSort2(uint_4, r_8*, int_4*);
[3678]263
264template void TabSortInd(uint_4, int_2*, uint_4*);
265template void TabSortInd(uint_4, uint_2*, uint_4*);
266template void TabSortInd(uint_4, int_4*, uint_4*);
267template void TabSortInd(uint_4, uint_4*, uint_4*);
268template void TabSortInd(uint_4, int_8*, uint_4*);
269template void TabSortInd(uint_4, uint_8*, uint_4*);
270template void TabSortInd(uint_4, r_4*, uint_4*);
271template void TabSortInd(uint_4, r_8*, uint_4*);
272#endif
273
274} // FIN namespace SOPHYA
Note: See TracBrowser for help on using the repository browser.