Changeset 3685 in Sophya for trunk/SophyaLib


Ignore:
Timestamp:
Nov 27, 2009, 2:54:18 PM (16 years ago)
Author:
cmv
Message:
  • changement du code C-Fortran like addressing par du code C++ purement C-like addressing (NR C++ version)

cmv, 27/11/2009

Location:
trunk/SophyaLib/NTools
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/NTools/nbtrixx.cc

    r3678 r3685  
    77namespace SOPHYA {
    88
    9 // attention, Numerical recipes prend des tableaux de 1 a n on remet
    10 // de 0 a n-1 en decramentant le pointeur du tableau d'entree
    11  
    12 //-------------------------------------------------------------
    13 #define SWAP_NBTRI(a,b) temp=(a);(a)=(b);(b)=temp;
    14 #define M 7
    15 #define NSTACK 50
     9typedef int_4 Int; // 32 bit integer
     10
     11//-------------------------------------------------------------
    1612template <class T>
    17 void TabSort(uint_4 n, T* arr_c)
     13void TabSort(uint_4 n, T* arr)
    1814{
    19         uint_4 i,ir=n,j,k,l=1,*istack;
    20         int jstack=0;
    21         T a,temp;
    22         T *arr = arr_c - 1; // MODIF
    23 
    24         //istack=lvector(1,NSTACK);
    25         istack = new uint_4[NSTACK+1]; // MODIF
     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;
    2620        for (;;) {
    2721                if (ir-l < M) {
    … …  
    3428                                arr[i+1]=a;
    3529                        }
    36                         if (jstack == 0) break;
     30                        if (jstack < 0) break;
    3731                        ir=istack[jstack--];
    3832                        l=istack[jstack--];
    3933                } else {
    4034                        k=(l+ir) >> 1;
    41                         SWAP_NBTRI(arr[k],arr[l+1])
     35                        _NR_SWAP(arr[k],arr[l+1]);
    4236                        if (arr[l] > arr[ir]) {
    43                                 SWAP_NBTRI(arr[l],arr[ir])
     37                                _NR_SWAP(arr[l],arr[ir]);
    4438                        }
    4539                        if (arr[l+1] > arr[ir]) {
    46                                 SWAP_NBTRI(arr[l+1],arr[ir])
     40                                _NR_SWAP(arr[l+1],arr[ir]);
    4741                        }
    4842                        if (arr[l] > arr[l+1]) {
    49                                 SWAP_NBTRI(arr[l],arr[l+1])
     43                                _NR_SWAP(arr[l],arr[l+1]);
    5044                        }
    5145                        i=l+1;
    … …  
    5650                                do j--; while (arr[j] > a);
    5751                                if (j < i) break;
    58                                 SWAP_NBTRI(arr[i],arr[j]);
     52                                _NR_SWAP(arr[i],arr[j]);
    5953                        }
    6054                        arr[l+1]=arr[j];
    6155                        arr[j]=a;
    6256                        jstack += 2;
    63                         if (jstack > NSTACK) throw("NSTACK too small in sort.");
     57                        if (jstack >= NSTACK) throw("NSTACK too small in sort.");
    6458                        if (ir-i+1 >= j-l) {
    6559                                istack[jstack]=ir;
    … …  
    7367                }
    7468        }
    75         //free_lvector(istack,1,NSTACK);
    76         delete [] istack;  // MODIF
     69        delete [] istack;
    7770}
    78 #undef M
    79 #undef NSTACK
    80 #undef SWAP_NBTRI
    81 
    82 //-------------------------------------------------------------
    83 #define SWAP_NBTRI(a,b) temp=(a);(a)=(b);(b)=temp;
    84 #define M 7
    85 #define NSTACK 50
    86 template <class T>
    87 void TabSort2(uint_4 n, T *arr_c, T *brr_c)
     71
     72//-------------------------------------------------------------
     73template <class T, class U>
     74void TabSort2(uint_4 n, T *arr, U *brr)
    8875{
    89         uint_4 i,ir=n,j,k,l=1,*istack;
    90         int jstack=0;
    91         T a,b,temp;
    92         T *arr = arr_c - 1; // MODIF
    93         T *brr = brr_c - 1; // MODIF
    94 
    95         //istack=lvector(1,NSTACK);
    96         istack = new uint_4[NSTACK+1]; // MODIF
     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;
    9782        for (;;) {
    9883                if (ir-l < M) {
    … …  
    10893                                brr[i+1]=b;
    10994                        }
    110                         if (!jstack) {
    111                                 //free_lvector(istack,1,NSTACK);
    112                                 delete [] istack;  // MODIF
    113                                 return;
    114                         }
    115                         ir=istack[jstack];
    116                         l=istack[jstack-1];
    117                         jstack -= 2;
     95                        if (jstack < 0) break;
     96                        ir=istack[jstack--];
     97                        l=istack[jstack--];
    11898                } else {
    11999                        k=(l+ir) >> 1;
    120                         SWAP_NBTRI(arr[k],arr[l+1])
    121                         SWAP_NBTRI(brr[k],brr[l+1])
     100                        _NR_SWAP(arr[k],arr[l+1]);
     101                        _NR_SWAP(brr[k],brr[l+1]);
    122102                        if (arr[l] > arr[ir]) {
    123                                 SWAP_NBTRI(arr[l],arr[ir])
    124                                 SWAP_NBTRI(brr[l],brr[ir])
     103                                _NR_SWAP(arr[l],arr[ir]);
     104                                _NR_SWAP(brr[l],brr[ir]);
    125105                        }
    126106                        if (arr[l+1] > arr[ir]) {
    127                                 SWAP_NBTRI(arr[l+1],arr[ir])
    128                                 SWAP_NBTRI(brr[l+1],brr[ir])
     107                                _NR_SWAP(arr[l+1],arr[ir]);
     108                                _NR_SWAP(brr[l+1],brr[ir]);
    129109                        }
    130110                        if (arr[l] > arr[l+1]) {
    131                                 SWAP_NBTRI(arr[l],arr[l+1])
    132                                 SWAP_NBTRI(brr[l],brr[l+1])
     111                                _NR_SWAP(arr[l],arr[l+1]);
     112                                _NR_SWAP(brr[l],brr[l+1]);
    133113                        }
    134114                        i=l+1;
    … …  
    140120                                do j--; while (arr[j] > a);
    141121                                if (j < i) break;
    142                                 SWAP_NBTRI(arr[i],arr[j])
    143                                 SWAP_NBTRI(brr[i],brr[j])
     122                                _NR_SWAP(arr[i],arr[j]);
     123                                _NR_SWAP(brr[i],brr[j]);
    144124                        }
    145125                        arr[l+1]=arr[j];
    … …  
    148128                        brr[j]=b;
    149129                        jstack += 2;
    150                         if (jstack > NSTACK) throw("NSTACK too small in sort2.");
     130                        if (jstack >= NSTACK) throw("NSTACK too small in sort2.");
    151131                        if (ir-i+1 >= j-l) {
    152132                                istack[jstack]=ir;
    … …  
    160140                }
    161141        }
     142        delete [] istack;
    162143}
    163 #undef M
    164 #undef NSTACK
    165 #undef SWAP_NBTRI
    166 
    167 //-------------------------------------------------------------
    168 #define SWAP_NBTRI(a,b) itemp=(a);(a)=(b);(b)=itemp;
    169 #define M 7
    170 #define NSTACK 50
     144
     145
     146//-------------------------------------------------------------
    171147template <class T>
    172 void TabSortInd(uint_4 n, T *arr_c, uint_4 *indx_c)
     148void TabSortInd(uint_4 n, T *arr, uint_4 *indx)
    173149{
    174         uint_4 i,indxt,ir=n,itemp,j,k,l=1;
    175         int jstack=0,*istack;
     150        const Int M=7,NSTACK=64;
     151        Int i,indxt,ir,j,k,jstack=-1,l=0;
    176152        T a;
    177         T *arr = arr_c - 1; // MODIF
    178         uint_4 *indx = indx_c - 1; // MODIF
    179 
    180         //istack=ivector(1,NSTACK);
    181         istack = new int[NSTACK+1]; // MODIF
    182         for (j=1;j<=n;j++) indx[j]=j;
     153        Int* istack = new Int[NSTACK];
     154        ir=n-1;
     155        for (j=0;j<(Int)n;j++) indx[j]=j;
    183156        for (;;) {
    184157                if (ir-l < M) {
    … …  
    192165                                indx[i+1]=indxt;
    193166                        }
    194                         if (jstack == 0) break;
     167                        if (jstack < 0) break;
    195168                        ir=istack[jstack--];
    196169                        l=istack[jstack--];
    197170                } else {
    198171                        k=(l+ir) >> 1;
    199                         SWAP_NBTRI(indx[k],indx[l+1]);
     172                        _NR_SWAP(indx[k],indx[l+1]);
    200173                        if (arr[indx[l]] > arr[indx[ir]]) {
    201                                 SWAP_NBTRI(indx[l],indx[ir])
     174                                _NR_SWAP(indx[l],indx[ir]);
    202175                        }
    203176                        if (arr[indx[l+1]] > arr[indx[ir]]) {
    204                                 SWAP_NBTRI(indx[l+1],indx[ir])
     177                                _NR_SWAP(indx[l+1],indx[ir]);
    205178                        }
    206179                        if (arr[indx[l]] > arr[indx[l+1]]) {
    207                                 SWAP_NBTRI(indx[l],indx[l+1])
     180                                _NR_SWAP(indx[l],indx[l+1]);
    208181                        }
    209182                        i=l+1;
    … …  
    215188                                do j--; while (arr[indx[j]] > a);
    216189                                if (j < i) break;
    217                                 SWAP_NBTRI(indx[i],indx[j])
     190                                _NR_SWAP(indx[i],indx[j]);
    218191                        }
    219192                        indx[l+1]=indx[j];
    220193                        indx[j]=indxt;
    221194                        jstack += 2;
    222                         if (jstack > NSTACK) throw("NSTACK too small in indexx.");
     195                        if (jstack >= NSTACK) throw("NSTACK too small in index.");
    223196                        if (ir-i+1 >= j-l) {
    224197                                istack[jstack]=ir;
    … …  
    232205                }
    233206        }
    234         //free_ivector(istack,1,NSTACK);
    235         delete [] istack;  // MODIF
    236         /* pour avoir un retour d'indice entre 0 et n-1 */
    237         for(j=0;j<n;j++) indx_c[j]--;  // MODIF
     207        delete [] istack;
    238208}
    239 #undef M
    240 #undef NSTACK
    241 #undef SWAP_NBTRI
     209
    242210
    243211//-------------------------------------------------------------
    … …  
    252220#pragma define_template TabSort<r_8>
    253221
    254 #pragma define_template TabSort2<int_2>
    255 #pragma define_template TabSort2<uint_2>
    256 #pragma define_template TabSort2<int_4>
    257 #pragma define_template TabSort2<uint_4>
    258 #pragma define_template TabSort2<int_8>
    259 #pragma define_template TabSort2<uint_8>
    260 #pragma define_template TabSort2<r_4>
    261 #pragma define_template TabSort2<r_8>
     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>
    262232
    263233#pragma define_template TabSortInd<int_2>
    … …  
    289259template void TabSort2(uint_4, r_4*, r_4*);
    290260template void TabSort2(uint_4, r_8*, r_8*);
     261template void TabSort2(uint_4, r_4*, int_4*);
     262template void TabSort2(uint_4, r_8*, int_4*);
    291263
    292264template void TabSortInd(uint_4, int_2*, uint_4*);

  • trunk/SophyaLib/NTools/nbtrixx.h

    r3678 r3685  
    66namespace SOPHYA {
    77
    8 template <class T>
    9 void TabSort(uint_4 n, T *arr_c);
     8template<class T>
     9inline void _NR_SWAP(T &a, T &b) {T dum=a; a=b; b=dum;}
    1010
    1111template <class T>
    12 void TabSort2(uint_4 n, T *arr_c, T *brr_c);
     12void TabSort(uint_4 n, T *arr);
     13
     14  template <class T, class U>
     15void TabSort2(uint_4 n, T *arr, U *brr);
    1316
    1417template <class T>
    15 void TabSortInd(uint_4 n, T *arr_c, uint_4 *indx_c);
     18void TabSortInd(uint_4 n, T *arr, uint_4 *indx);
    1619
    1720}
Note: See TracChangeset for help on using the changeset viewer.