This document is "essentially" an annotated version of a simple sort code. It produces a vector norder such that an intial vector Energy is ordered. That is
Energy(norder(i))<=Energy(norderIi+1)) for i = 1 to Ntotal –1.
The routines are SORT.FOR and sort.c are self contained. They are tested by tsort.c and TSORT.FOR which depend on random numbers. In addtion they are derived from locate routines. Short descriptions of these are in ..\..\definitions\Random\Random.doc and Locate.htm.
1. #Description of the routine sort
2. #Description of the
calling routine
3. #Annotated test code
TSORT.FOR
4. #C code
5. #File Sort
The subroutine sort(norder,energy,ntotal)
produces a vector norder(n) such that energy(norder(j)).<=energy(norder(j+1)) for 1<=j<=ntotal-1. The routine
begins with a modification of the routine Locate, a short, fast look up routine for unevenly
spaced data given in Press[1] (Fortran 1986) p 90.
In C notation, but energy having Fortran spacing from 1 to n, locate(double x,double energy[],int j) returns j such that energy[j]<=x<=energy[j+1]. The value j=0 ==> x<energy[1]. The value
j=n ==> energy[n] < x.
The
routine TLOC has the vector norder[n] as an
additional input so that energy(norder(j)).<=x<=energy(norder(j+1))
with the same meanings for 0 and n as before.
The routine sort uses these in succession for each value of energy.
SUBROUTINE SORT(NORDER,ENERGY,NTOTAL)
INTEGER*2 NORDER(NTOTAL) The
size of this effects the speed
DIMENSION ENERGY(NTOTAL) It may
be desirable to make this real*8 for many applications.
IF(NTOTAL.GT.32768)THEN
PRINT*,' FROM SORT NTOTAL =',NTOTAL This is to warn about an unexpected event
PRINT*,' IS TOO LARGE FOR INTEGER*2' The single if
requires negligible time
STOP
ENDIF
NORDER(1)=1
DO I=2,NTOTAL
CALL TLOC(NORDER,ENERGY,ENERGY(I),I-1,J)
DO K=I-1,J+1,-1 This is the reason for using integer*2, if all elements
NORDER(K+1)=NORDER(K) of norder fit in
a very small space the NTOTAL2 updating
ENDDO is much faster.
NORDER(J+1)=I
ENDDO
RETURN
END
SUBROUTINE TLOC(NORDER,ENERGY,X,NMAX,J)
C
CODE RETURNS J SUCH THAT ENERGY(NORDER(J))<=X<=ENERGY(NORDER(J+1)
C
J=0 ==> X < ENERGY(NORDER(1)
C
J=NMAX ==> X > ENERGY(NORDER(NMAX)) for
details see Locate
INTEGER*2 NORDER(NMAX)
DIMENSION ENERGY(NMAX)
JL=0
JU=NMAX+1
10 IF(JU-JL.LE.1)GOTO
20
JM=(JU+JL)/2
IF(X.GT.ENERGY(NORDER(JM)))THEN
JL=JM
ELSE
JU=JM
ENDIF
GOTO 10
20 J=JL
RETURN
END
The advantage of using the
auxiliary vector norder, is in the inner loop where it is only norder
that gets moved, rather than energy and other vectors in the same original
order as energy. Thus if peak areas and
the errors in them are orignally such that energy(I), area(I), and error(I) refer to the same peak,
these can be output in energy order by
Call sort(norder,energy,ntotal)
Do I=1,ntotal
Print*,energy(norder(I)),area(norder(I)),error(norder(I))
Enddo
If one wishes to return only
energy, area and error to a main code with these ordered by energy, there is
still a saving to the norder step. This comes about because the inner loop in
sort which is 2d involves many fewer moves and is much smaller in size and thus
requires fewer memory calls. The final
step is most easily made with a single temp array so that in the calling
routine there is originally
Dimension
E(ntotal),area(ntotal),error(ntotal),itype(ntotal)
These are presumed ordered in
an arbitrary way, but E(3), area(3), error(3), and itype(3) still all refer to the same peak of interest. Note that there must also be in this routine
Integer*2 norder(ntotal)
And
dimension temp(ntotal)
Then the code finds values
for norder by
call sort(norder,energy,ntotal)
The reordering is in the two
steps
Do I=1,ntotal
Temp(I)=e(norder(I))
Enddo
Do I=1,ntotal
E(I)=temp(I)
Enddo
Two steps are used owing to
the fact that the E's used in temp are at numerous different I's making it extremely difficult, if possible at all, to
reorder as they are moved --- see Press about the bubble sort. The vector temp can then be re-used for area,
and error. If space is tight the vector ntemp can be defined as
Dimension ntemp(ntotal)
Equivalence (ntemp(1),temp(1))
Now the loop structure above
with temp è ntemp and E è itype works for the last
vector also. Each of these steps is
linear in ntotal and thus fast compared to the sort.
parameter (ntotal=35) This
is the number of values that fit easily on my screen.
integer*2 norder(ntotal)
dimension e(ntotal),temp(ntotal)
COMMON/RAN/IX,IY,ITAB(128)
call rseed(3157891,4821359) See ../../Fittery/ArtDat/Random.doc for
documentation
5 print*,ix,iy This prints the starting values of the seed for restarting at the
call rseed(ix,iy) error if things go wrong.
do i=1,ntotal
e(i)=rndmf(1.)
This makes a random set of e's for testing.
enddo
call SORT(NORDER,E,NTOTAL)
do i=1,ntotal
temp(i)=e(norder(i))
enddo
do i=1,ntotal
e(i)=temp(i)
enddo
do i=1,ntotal
print*,i,e(i)
enddo
read(*,*)itest
if(itest.eq.1)goto
5 This
repeats the test until the author gets tired.
end
c$include sort sort.for
c$include random random.for
1386206595 4821359
1 0.0804947
2 0.0815424
3 0.0858933
4 0.0916346
5 0.0920604
6 0.1017884
7 0.1122877
8 0.1290050
9 0.1450515
10 0.1472510
11 0.1550671
12 0.2034720
13 0.2243175
14 0.2760789
15 0.2859344
16 0.3447586
17 0.3797226
18 0.4275883
19 0.4374947
20 0.4976539
21 0.5023293
22 0.5431748
23 0.5517566
24 0.5822303
25 0.6324221
26 0.6494815
27 0.7218786
28 0.7326567
29 0.7695072
30 0.7968472
31 0.8648781
32 0.9053845
33 0.9581237
34 0.9732007
35 0.9733715
#include <string.h>
#include <stdio.h>
#define NTOTAL 23 ß
Fills the screen
int iloc(int norder[],double energy[],double x,int nmax);
int rseed(int ia,int ib);
double rndmf();
void sort(int norder[],double energy[],int ntotal);
void sort(int norder[],double energy[],int ntotal)
{int i,j,k;
norder[0]=0;
for (i=1;i<ntotal;++i)
{j=iloc(norder,energy,energy[i],i);
if(i-1 >= j+1)
{for(k=i-1;k>=j+1;--k){norder[k+1]=norder[k];}}
norder[j+1]=i;}
return;}
int iloc(int norder[],double energy[],double x,int nmax)
/* code returns iloc such that energy(norder(iloc))<=x<=energy(norder(iloc+1)
j=-1 ==> x < energy(norder(0)
j=nmax ==> x > energy(norder(nmax))*/
{int jl,ju,jm;
jl=-1;
ju=nmax;
while (ju-jl>1)
{jm=(ju+jl) >> 1; /* the /2 has become a right shift*/
if(x > energy[norder[jm]])jl=jm;
else ju=jm;
}
return jl;}
void main(void)
{int norder[NTOTAL];
int ix=3157891,iy=4821359;
int i;
double e[NTOTAL],temp[NTOTAL]; ß will complain about temp – it used elsewhere
char ans[80];
rseed(ix,iy);
begcheck:
for(i=0;i<NTOTAL;++i)e[i]=rndmf();
sort(norder,e,NTOTAL);
for(i=0;i<NTOTAL;++i)
{printf(" %d %lg \n",i,e[norder[i]]);}
scanf("%s",ans);
if(strncmp("stop",ans,4)!=0)goto begcheck;
return;}
The data FT is stored as binary data in a file named
FT. This file can be opened, all the
data read and the 1037’th value placed in the variable
F by the statements.
DIMENSION
FT(32756)
OPEN( 1,
FILE=’FT’,ACCESS='SEQUENTIAL',FORM='UNFORMATTED',RECORDTYPE='FIXED')
READ(1)FT
F=FT(1037)
This same value can be found
by opening the file for ‘direct’ access.
Pascal direct access is described in PascalDA.htm.
OPEN(1,FILE=’FT’,ACCESS='DIRECT',STATUS='UNKNOWN',FORM='UNFORMATTED',RECL=8)
READ(1,REC=I)F
If only a few values of F are desired, the second method
avoids tying up a large amount of the codes memory space for dead storage and
ultimately reads only a few values from the disk. If many values are needed, memory allocation
and de-allocation appear necessary owing to the slow speed of disk reading relative
to memory look-ups. In the real world,
the computer makes a buffer for the reads and the user over-dimensions the
files, so the competition is actually between the computer’s ability to find
which sections of the file should be set into memory and the computer’s ability
to page out the code memory that is not used.
A big advantage in debugging comes from the fact that ..\Watfor.doc .htm can be used with the smaller direct access code to make
sure that the code has no undefined variables.
The code DAsort.wpj
begins with a single direct access file that consists of the energies and
intensities of 4 nuclide combinations. Raw.UNF This file contains 606 entries and has a length of 14,544
bytes. These numbers are important. Direct access files are pure binary. The end of file, EOF, marker tends to be
ignored and one can end up sorting the code or even the operating system.
1. K-40.IIS
2. ~Th-232.IIS
3. ~Ra-226.IIS
4. Eu-152.IIS
tdasort.for - Lines transferred = 30 – code in TEST1.zip
RUN
ENTER THE NAME OF THE INPUT FILE
RAW.UNF
I
EP, AJ, I1,I2
1
0.9498310000000000E+03
0.5758477568811740E-01 1 61
2
0.1460830000000000E+04
0.1100000000000000E+02 1 61
3
0.1539770000000000E+03
0.7220000000000000E+00 1 62
4
0.1739640000000000E+03
0.3500000000000000E-01 1 62
5
0.1845400000000000E+03
0.7000000000000000E-01 1 62
6
0.1913530000000000E+03
0.1230000000000000E+00 1 62
The first six lines of the
data are shown above. The first two are
from K-40. the next four are from the Th-232 decay
chain. The first column is the energy of
the gamma ray, while the second is the relative intensity (the number of gammas
of this energy per 100 nuclear decays).
The third and fourth columns are the attenuation type and the constant
number used in ..\..\Fittery\robfit\TemplFit\Welcome.htm, The code described
here will produce either a file ordered by energy or by relative intensity.
The final code along with RAW.UNF is in DAsort.zip. This is set up to make a stand alone sort
routine. The system calls are set for
Watcom. Wsystem
converts these to Watfor and PSSYS converts these to
Power Station Fortran.
[1] W.H.
Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical
Recipes,