source: LMDZ4/trunk/libf/cosp/mrgrnk.F90 @ 2615

Last change on this file since 2615 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
File size: 11.2 KB
Line 
1Module m_mrgrnk
2Integer, Parameter :: kdp = selected_real_kind(15)
3public :: mrgrnk
4private :: kdp
5private :: I_mrgrnk, D_mrgrnk
6interface mrgrnk
7  module procedure D_mrgrnk, I_mrgrnk
8end interface mrgrnk
9contains
10
11Subroutine D_mrgrnk (XDONT, IRNGT)
12! __________________________________________________________
13!   MRGRNK = Merge-sort ranking of an array
14!   For performance reasons, the first 2 passes are taken
15!   out of the standard loop, and use dedicated coding.
16! __________________________________________________________
17! __________________________________________________________
18      Real (kind=kdp), Dimension (:), Intent (In) :: XDONT
19      Integer, Dimension (:), Intent (Out) :: IRNGT
20! __________________________________________________________
21      Real (kind=kdp) :: XVALA, XVALB
22!
23      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
24      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
25      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
26!
27      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
28      Select Case (NVAL)
29      Case (:0)
30         Return
31      Case (1)
32         IRNGT (1) = 1
33         Return
34      Case Default
35         Continue
36      End Select
37!
38!  Fill-in the index array, creating ordered couples
39!
40      Do IIND = 2, NVAL, 2
41         If (XDONT(IIND-1) <= XDONT(IIND)) Then
42            IRNGT (IIND-1) = IIND - 1
43            IRNGT (IIND) = IIND
44         Else
45            IRNGT (IIND-1) = IIND
46            IRNGT (IIND) = IIND - 1
47         End If
48      End Do
49      If (Modulo(NVAL, 2) /= 0) Then
50         IRNGT (NVAL) = NVAL
51      End If
52!
53!  We will now have ordered subsets A - B - A - B - ...
54!  and merge A and B couples into     C   -   C   - ...
55!
56      LMTNA = 2
57      LMTNC = 4
58!
59!  First iteration. The length of the ordered subsets goes from 2 to 4
60!
61      Do
62         If (NVAL <= 2) Exit
63!
64!   Loop on merges of A and B into C
65!
66         Do IWRKD = 0, NVAL - 1, 4
67            If ((IWRKD+4) > NVAL) Then
68               If ((IWRKD+2) >= NVAL) Exit
69!
70!   1 2 3
71!
72               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
73!
74!   1 3 2
75!
76               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
77                  IRNG2 = IRNGT (IWRKD+2)
78                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
79                  IRNGT (IWRKD+3) = IRNG2
80!
81!   3 1 2
82!
83               Else
84                  IRNG1 = IRNGT (IWRKD+1)
85                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
86                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
87                  IRNGT (IWRKD+2) = IRNG1
88               End If
89               Exit
90            End If
91!
92!   1 2 3 4
93!
94            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
95!
96!   1 3 x x
97!
98            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
99               IRNG2 = IRNGT (IWRKD+2)
100               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
101               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
102!   1 3 2 4
103                  IRNGT (IWRKD+3) = IRNG2
104               Else
105!   1 3 4 2
106                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
107                  IRNGT (IWRKD+4) = IRNG2
108               End If
109!
110!   3 x x x
111!
112            Else
113               IRNG1 = IRNGT (IWRKD+1)
114               IRNG2 = IRNGT (IWRKD+2)
115               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
116               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
117                  IRNGT (IWRKD+2) = IRNG1
118                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
119!   3 1 2 4
120                     IRNGT (IWRKD+3) = IRNG2
121                  Else
122!   3 1 4 2
123                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
124                     IRNGT (IWRKD+4) = IRNG2
125                  End If
126               Else
127!   3 4 1 2
128                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
129                  IRNGT (IWRKD+3) = IRNG1
130                  IRNGT (IWRKD+4) = IRNG2
131               End If
132            End If
133         End Do
134!
135!  The Cs become As and Bs
136!
137         LMTNA = 4
138         Exit
139      End Do
140!
141!  Iteration loop. Each time, the length of the ordered subsets
142!  is doubled.
143!
144      Do
145         If (LMTNA >= NVAL) Exit
146         IWRKF = 0
147         LMTNC = 2 * LMTNC
148!
149!   Loop on merges of A and B into C
150!
151         Do
152            IWRK = IWRKF
153            IWRKD = IWRKF + 1
154            JINDA = IWRKF + LMTNA
155            IWRKF = IWRKF + LMTNC
156            If (IWRKF >= NVAL) Then
157               If (JINDA >= NVAL) Exit
158               IWRKF = NVAL
159            End If
160            IINDA = 1
161            IINDB = JINDA + 1
162!
163!   Shortcut for the case when the max of A is smaller
164!   than the min of B. This line may be activated when the
165!   initial set is already close to sorted.
166!
167!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
168!
169!  One steps in the C subset, that we build in the final rank array
170!
171!  Make a copy of the rank array for the merge iteration
172!
173            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
174!
175            XVALA = XDONT (JWRKT(IINDA))
176            XVALB = XDONT (IRNGT(IINDB))
177!
178            Do
179               IWRK = IWRK + 1
180!
181!  We still have unprocessed values in both A and B
182!
183               If (XVALA > XVALB) Then
184                  IRNGT (IWRK) = IRNGT (IINDB)
185                  IINDB = IINDB + 1
186                  If (IINDB > IWRKF) Then
187!  Only A still with unprocessed values
188                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
189                     Exit
190                  End If
191                  XVALB = XDONT (IRNGT(IINDB))
192               Else
193                  IRNGT (IWRK) = JWRKT (IINDA)
194                  IINDA = IINDA + 1
195                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
196                  XVALA = XDONT (JWRKT(IINDA))
197               End If
198!
199            End Do
200         End Do
201!
202!  The Cs become As and Bs
203!
204         LMTNA = 2 * LMTNA
205      End Do
206!
207      Return
208!
209End Subroutine D_mrgrnk
210
211Subroutine I_mrgrnk (XDONT, IRNGT)
212! __________________________________________________________
213!   MRGRNK = Merge-sort ranking of an array
214!   For performance reasons, the first 2 passes are taken
215!   out of the standard loop, and use dedicated coding.
216! __________________________________________________________
217! __________________________________________________________
218      Integer, Dimension (:), Intent (In)  :: XDONT
219      Integer, Dimension (:), Intent (Out) :: IRNGT
220! __________________________________________________________
221      Integer :: XVALA, XVALB
222!
223      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
224      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
225      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
226!
227      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
228      Select Case (NVAL)
229      Case (:0)
230         Return
231      Case (1)
232         IRNGT (1) = 1
233         Return
234      Case Default
235         Continue
236      End Select
237!
238!  Fill-in the index array, creating ordered couples
239!
240      Do IIND = 2, NVAL, 2
241         If (XDONT(IIND-1) <= XDONT(IIND)) Then
242            IRNGT (IIND-1) = IIND - 1
243            IRNGT (IIND) = IIND
244         Else
245            IRNGT (IIND-1) = IIND
246            IRNGT (IIND) = IIND - 1
247         End If
248      End Do
249      If (Modulo(NVAL, 2) /= 0) Then
250         IRNGT (NVAL) = NVAL
251      End If
252!
253!  We will now have ordered subsets A - B - A - B - ...
254!  and merge A and B couples into     C   -   C   - ...
255!
256      LMTNA = 2
257      LMTNC = 4
258!
259!  First iteration. The length of the ordered subsets goes from 2 to 4
260!
261      Do
262         If (NVAL <= 2) Exit
263!
264!   Loop on merges of A and B into C
265!
266         Do IWRKD = 0, NVAL - 1, 4
267            If ((IWRKD+4) > NVAL) Then
268               If ((IWRKD+2) >= NVAL) Exit
269!
270!   1 2 3
271!
272               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
273!
274!   1 3 2
275!
276               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
277                  IRNG2 = IRNGT (IWRKD+2)
278                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
279                  IRNGT (IWRKD+3) = IRNG2
280!
281!   3 1 2
282!
283               Else
284                  IRNG1 = IRNGT (IWRKD+1)
285                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
286                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
287                  IRNGT (IWRKD+2) = IRNG1
288               End If
289               Exit
290            End If
291!
292!   1 2 3 4
293!
294            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
295!
296!   1 3 x x
297!
298            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
299               IRNG2 = IRNGT (IWRKD+2)
300               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
301               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
302!   1 3 2 4
303                  IRNGT (IWRKD+3) = IRNG2
304               Else
305!   1 3 4 2
306                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
307                  IRNGT (IWRKD+4) = IRNG2
308               End If
309!
310!   3 x x x
311!
312            Else
313               IRNG1 = IRNGT (IWRKD+1)
314               IRNG2 = IRNGT (IWRKD+2)
315               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
316               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
317                  IRNGT (IWRKD+2) = IRNG1
318                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
319!   3 1 2 4
320                     IRNGT (IWRKD+3) = IRNG2
321                  Else
322!   3 1 4 2
323                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
324                     IRNGT (IWRKD+4) = IRNG2
325                  End If
326               Else
327!   3 4 1 2
328                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
329                  IRNGT (IWRKD+3) = IRNG1
330                  IRNGT (IWRKD+4) = IRNG2
331               End If
332            End If
333         End Do
334!
335!  The Cs become As and Bs
336!
337         LMTNA = 4
338         Exit
339      End Do
340!
341!  Iteration loop. Each time, the length of the ordered subsets
342!  is doubled.
343!
344      Do
345         If (LMTNA >= NVAL) Exit
346         IWRKF = 0
347         LMTNC = 2 * LMTNC
348!
349!   Loop on merges of A and B into C
350!
351         Do
352            IWRK = IWRKF
353            IWRKD = IWRKF + 1
354            JINDA = IWRKF + LMTNA
355            IWRKF = IWRKF + LMTNC
356            If (IWRKF >= NVAL) Then
357               If (JINDA >= NVAL) Exit
358               IWRKF = NVAL
359            End If
360            IINDA = 1
361            IINDB = JINDA + 1
362!
363!   Shortcut for the case when the max of A is smaller
364!   than the min of B. This line may be activated when the
365!   initial set is already close to sorted.
366!
367!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
368!
369!  One steps in the C subset, that we build in the final rank array
370!
371!  Make a copy of the rank array for the merge iteration
372!
373            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
374!
375            XVALA = XDONT (JWRKT(IINDA))
376            XVALB = XDONT (IRNGT(IINDB))
377!
378            Do
379               IWRK = IWRK + 1
380!
381!  We still have unprocessed values in both A and B
382!
383               If (XVALA > XVALB) Then
384                  IRNGT (IWRK) = IRNGT (IINDB)
385                  IINDB = IINDB + 1
386                  If (IINDB > IWRKF) Then
387!  Only A still with unprocessed values
388                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
389                     Exit
390                  End If
391                  XVALB = XDONT (IRNGT(IINDB))
392               Else
393                  IRNGT (IWRK) = JWRKT (IINDA)
394                  IINDA = IINDA + 1
395                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
396                  XVALA = XDONT (JWRKT(IINDA))
397               End If
398!
399            End Do
400         End Do
401!
402!  The Cs become As and Bs
403!
404         LMTNA = 2 * LMTNA
405      End Do
406!
407      Return
408!
409End Subroutine I_mrgrnk
410end module m_mrgrnk
Note: See TracBrowser for help on using the repository browser.