1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
*DECK ALINDX
SUBROUTINE ALINDX(n,arr,indx)
*
*-----------------------------------------------------------------------
*
*Purpose:
* indexes an array arr(1:n) such that abs(arr(indx(j))) is in descending
* order for j=1:n.
*
*Copyright:
* Copyright (C) 2020 Ecole Polytechnique de MontCOMPLEX(KIND=8)
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): A. Hebert
*
*Parameters: input
* n size of the array a.
* a array a.
*
*Parameters: output
* indx indexing vector.
*
*-----------------------------------------------------------------------
*
*----
* subroutine arguments
*----
INTEGER n,indx(n)
COMPLEX(KIND=8) arr(n)
*----
* local variables
*----
INTEGER nstack,i,indxt,ir,itemp,j,jstack,k,l
COMPLEX(KIND=8) a
INTEGER, parameter :: M=7
INTEGER, allocatable, dimension(:) :: istack
*
nstack=2*ceiling(log(real(n))/log(2.0))
allocate(istack(nstack))
do j=1,n
indx(j)=j
enddo
jstack=0
l=1
ir=n
1 if(ir-l.lt.M) then
do j=l+1,ir
indxt=indx(j)
a=arr(indxt)
do i=j-1,l,-1
if(absc(arr(indx(i))).ge.absc(a)) goto 2
indx(i+1)=indx(i)
enddo
i=l-1
2 indx(i+1)=indxt
enddo
if(jstack.eq.0)go to 5
ir=istack(jstack)
l=istack(jstack-1)
jstack=jstack-2
else
k=(l+ir)/2
itemp=indx(k)
indx(k)=indx(l+1)
indx(l+1)=itemp
if(absc(arr(indx(l))).lt.absc(arr(indx(ir))))then
itemp=indx(l)
indx(l)=indx(ir)
indx(ir)=itemp
endif
if(absc(arr(indx(l+1))).lt.absc(arr(indx(ir)))) then
itemp=indx(l+1)
indx(l+1)=indx(ir)
indx(ir)=itemp
endif
if(absc(arr(indx(l))).lt.absc(arr(indx(l+1)))) then
itemp=indx(l)
indx(l)=indx(l+1)
indx(l+1)=itemp
endif
i=l+1
j=ir
indxt=indx(l+1)
a=arr(indxt)
3 i=i+1
if(absc(arr(indx(i))).gt.absc(a)) goto 3
4 j=j-1
if(absc(arr(indx(j))).lt.absc(a)) goto 4
if(j.ge.i) then
itemp=indx(i)
indx(i)=indx(j)
indx(j)=itemp
goto 3
endif
indx(l+1)=indx(j)
indx(j)=indxt
jstack=jstack+2
if(jstack.gt.nstack) call XABORT('ALINDX: nstack too small.')
if(ir-i+1.ge.j-l) then
istack(jstack)=ir
istack(jstack-1)=i
ir=j-1
else
istack(jstack)=j-1
istack(jstack-1)=l
l=i
endif
endif
goto 1
5 deallocate(istack)
return
*
contains
function absc(val) result(aval)
! definition of combined value for complex numbers
COMPLEX(kind=8), intent(in) :: val
REAL(kind=8) :: aval
aval=100.0d0*real(val)+aimag(val)
end function absc
END
|
