Context Navigation

mrgrnk.F90 @ 5159

Last change on this file since 5159 was 5159, checked in by abarral, 6 months ago
Put dimensions.h and paramet.h into modules
File size: 18.5 KB

Rev	Line
[3358]	1	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	2	! Copyright (c) 2015, Regents of the University of Colorado
	3	! All rights reserved.
[5159]	4
[3358]	5	! Redistribution and use in source and binary forms, with or without modification, are
	6	! permitted provided that the following conditions are met:
[5159]	7
[3358]	8	! 1. Redistributions of source code must retain the above copyright notice, this list of
	9	! conditions and the following disclaimer.
[5159]	10
[3358]	11	! 2. Redistributions in binary form must reproduce the above copyright notice, this list
	12	! of conditions and the following disclaimer in the documentation and/or other
	13	! materials provided with the distribution.
[5159]	14
[3358]	15	! 3. Neither the name of the copyright holder nor the names of its contributors may be
	16	! used to endorse or promote products derived from this software without specific prior
	17	! written permission.
[5159]	18
[3358]	19	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
	20	! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
	21	! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
	22	! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
	23	! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
	24	! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
	25	! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
	26	! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
	27	! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
[5159]	28
[3358]	29	! History:
	30	! May 2015: Dustin Swales - Modified for COSPv2.0
	31	!
	32	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	33	Module m_mrgrnk
	34	USE COSP_KINDS, ONLY: wp
	35	Integer, Parameter :: kdp = selected_real_kind(15)
	36	public :: mrgrnk
	37	private :: kdp
	38	private :: R_mrgrnk, I_mrgrnk, D_mrgrnk
	39
	40	interface mrgrnk
	41	! module procedure D_mrgrnk, R_mrgrnk, I_mrgrnk
	42	module procedure R_mrgrnk, I_mrgrnk
	43
	44	end interface
	45	contains
	46
	47	Subroutine D_mrgrnk (XDONT, IRNGT)
	48	! __________________________________________________________
	49	! MRGRNK = Merge-sort ranking of an array
	50	! For performance reasons, the first 2 passes are taken
	51	! out of the standard loop, and use dedicated coding.
	52	! __________________________________________________________
	53	! __________________________________________________________
	54	Real (wp), Dimension (:), Intent (In) :: XDONT
	55	Integer, Dimension (:), Intent (Out) :: IRNGT
	56	! __________________________________________________________
	57	Real (wp) :: XVALA, XVALB
[5159]	58
[3358]	59	Integer, Dimension (SIZE(IRNGT)) :: JWRKT
	60	Integer :: LMTNA, LMTNC, IRNG1, IRNG2
	61	Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
[5159]	62
[3358]	63	NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
	64	Select Case (NVAL)
	65	Case (:0)
	66	Return
	67	Case (1)
	68	IRNGT (1) = 1
	69	Return
[5095]	70	Case Default
	71	Continue
[3358]	72	End Select
[5159]	73
[3358]	74	! Fill-in the index array, creating ordered couples
[5159]	75
[3358]	76	Do IIND = 2, NVAL, 2
	77	If (XDONT(IIND-1) <= XDONT(IIND)) Then
	78	IRNGT (IIND-1) = IIND - 1
	79	IRNGT (IIND) = IIND
	80	Else
	81	IRNGT (IIND-1) = IIND
	82	IRNGT (IIND) = IIND - 1
	83	End If
	84	End Do
	85	If (Modulo(NVAL, 2) /= 0) Then
	86	IRNGT (NVAL) = NVAL
	87	End If
[5159]	88
[3358]	89	! We will now have ordered subsets A - B - A - B - ...
	90	! and merge A and B couples into C - C - ...
[5159]	91
[3358]	92	LMTNA = 2
	93	LMTNC = 4
[5159]	94
[3358]	95	! First iteration. The length of the ordered subsets goes from 2 to 4
[5159]	96
[3358]	97	Do
	98	If (NVAL <= 2) Exit
[5159]	99
[3358]	100	! Loop on merges of A and B into C
[5159]	101
[3358]	102	Do IWRKD = 0, NVAL - 1, 4
	103	If ((IWRKD+4) > NVAL) Then
	104	If ((IWRKD+2) >= NVAL) Exit
[5159]	105
[3358]	106	! 1 2 3
[5159]	107
[3358]	108	If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
[5159]	109
[3358]	110	! 1 3 2
[5159]	111
[3358]	112	If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
	113	IRNG2 = IRNGT (IWRKD+2)
	114	IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
	115	IRNGT (IWRKD+3) = IRNG2
[5159]	116
[3358]	117	! 3 1 2
[5159]	118
[3358]	119	Else
	120	IRNG1 = IRNGT (IWRKD+1)
	121	IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
	122	IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
	123	IRNGT (IWRKD+2) = IRNG1
	124	End If
	125	Exit
	126	End If
[5159]	127
[3358]	128	! 1 2 3 4
[5159]	129
[3358]	130	If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
[5159]	131
[3358]	132	! 1 3 x x
[5159]	133
[3358]	134	If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
	135	IRNG2 = IRNGT (IWRKD+2)
	136	IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
	137	If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
	138	! 1 3 2 4
	139	IRNGT (IWRKD+3) = IRNG2
	140	Else
	141	! 1 3 4 2
	142	IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
	143	IRNGT (IWRKD+4) = IRNG2
	144	End If
[5159]	145
[3358]	146	! 3 x x x
[5159]	147
[3358]	148	Else
	149	IRNG1 = IRNGT (IWRKD+1)
	150	IRNG2 = IRNGT (IWRKD+2)
	151	IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
	152	If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
	153	IRNGT (IWRKD+2) = IRNG1
	154	If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
	155	! 3 1 2 4
	156	IRNGT (IWRKD+3) = IRNG2
	157	Else
	158	! 3 1 4 2
	159	IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
	160	IRNGT (IWRKD+4) = IRNG2
	161	End If
	162	Else
	163	! 3 4 1 2
	164	IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
	165	IRNGT (IWRKD+3) = IRNG1
	166	IRNGT (IWRKD+4) = IRNG2
	167	End If
	168	End If
	169	End Do
[5159]	170
[3358]	171	! The Cs become As and Bs
[5159]	172
[3358]	173	LMTNA = 4
	174	Exit
	175	End Do
[5159]	176
[3358]	177	! Iteration loop. Each time, the length of the ordered subsets
	178	! is doubled.
[5159]	179
[3358]	180	Do
	181	If (LMTNA >= NVAL) Exit
	182	IWRKF = 0
	183	LMTNC = 2 * LMTNC
[5159]	184
[3358]	185	! Loop on merges of A and B into C
[5159]	186
[3358]	187	Do
	188	IWRK = IWRKF
	189	IWRKD = IWRKF + 1
	190	JINDA = IWRKF + LMTNA
	191	IWRKF = IWRKF + LMTNC
	192	If (IWRKF >= NVAL) Then
	193	If (JINDA >= NVAL) Exit
	194	IWRKF = NVAL
	195	End If
	196	IINDA = 1
	197	IINDB = JINDA + 1
[5159]	198
[3358]	199	! Shortcut for the case when the max of A is smaller
	200	! than the min of B. This line may be activated when the
	201	! initial set is already close to sorted.
[5159]	202
[3358]	203	! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
[5159]	204
[3358]	205	! One steps in the C subset, that we build in the final rank array
[5159]	206
[3358]	207	! Make a copy of the rank array for the merge iteration
[5159]	208
[3358]	209	JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
[5159]	210
[3358]	211	XVALA = XDONT (JWRKT(IINDA))
	212	XVALB = XDONT (IRNGT(IINDB))
[5159]	213
[3358]	214	Do
	215	IWRK = IWRK + 1
[5159]	216
[3358]	217	! We still have unprocessed values in both A and B
[5159]	218
[3358]	219	If (XVALA > XVALB) Then
	220	IRNGT (IWRK) = IRNGT (IINDB)
	221	IINDB = IINDB + 1
	222	If (IINDB > IWRKF) Then
	223	! Only A still with unprocessed values
	224	IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
	225	Exit
	226	End If
	227	XVALB = XDONT (IRNGT(IINDB))
	228	Else
	229	IRNGT (IWRK) = JWRKT (IINDA)
	230	IINDA = IINDA + 1
	231	If (IINDA > LMTNA) Exit! Only B still with unprocessed values
	232	XVALA = XDONT (JWRKT(IINDA))
	233	End If
[5159]	234
[3358]	235	End Do
	236	End Do
[5159]	237
[3358]	238	! The Cs become As and Bs
[5159]	239
[3358]	240	LMTNA = 2 * LMTNA
	241	End Do
[5159]	242
[3358]	243	Return
[5159]	244
[3358]	245	End Subroutine D_mrgrnk
	246
	247	Subroutine R_mrgrnk (XDONT, IRNGT)
	248	! __________________________________________________________
	249	! MRGRNK = Merge-sort ranking of an array
	250	! For performance reasons, the first 2 passes are taken
	251	! out of the standard loop, and use dedicated coding.
	252	! __________________________________________________________
	253	! _________________________________________________________
	254	Real(wp), Dimension (:), Intent (In) :: XDONT
	255	Integer, Dimension (:), Intent (Out) :: IRNGT
	256	! __________________________________________________________
	257	Real(wp) :: XVALA, XVALB
[5159]	258
[3358]	259	Integer, Dimension (SIZE(IRNGT)) :: JWRKT
	260	Integer :: LMTNA, LMTNC, IRNG1, IRNG2
	261	Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
[5159]	262
[3358]	263	NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
	264	Select Case (NVAL)
	265	Case (:0)
	266	Return
	267	Case (1)
	268	IRNGT (1) = 1
	269	Return
[5095]	270	Case Default
	271	Continue
[3358]	272	End Select
[5159]	273
[3358]	274	! Fill-in the index array, creating ordered couples
[5159]	275
[3358]	276	Do IIND = 2, NVAL, 2
	277	If (XDONT(IIND-1) <= XDONT(IIND)) Then
	278	IRNGT (IIND-1) = IIND - 1
	279	IRNGT (IIND) = IIND
	280	Else
	281	IRNGT (IIND-1) = IIND
	282	IRNGT (IIND) = IIND - 1
	283	End If
	284	End Do
	285	If (Modulo(NVAL, 2) /= 0) Then
	286	IRNGT (NVAL) = NVAL
	287	End If
[5159]	288
[3358]	289	! We will now have ordered subsets A - B - A - B - ...
	290	! and merge A and B couples into C - C - ...
[5159]	291
[3358]	292	LMTNA = 2
	293	LMTNC = 4
[5159]	294
[3358]	295	! First iteration. The length of the ordered subsets goes from 2 to 4
[5159]	296
[3358]	297	Do
	298	If (NVAL <= 2) Exit
[5159]	299
[3358]	300	! Loop on merges of A and B into C
[5159]	301
[3358]	302	Do IWRKD = 0, NVAL - 1, 4
	303	If ((IWRKD+4) > NVAL) Then
	304	If ((IWRKD+2) >= NVAL) Exit
[5159]	305
[3358]	306	! 1 2 3
[5159]	307
[3358]	308	If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
[5159]	309
[3358]	310	! 1 3 2
[5159]	311
[3358]	312	If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
	313	IRNG2 = IRNGT (IWRKD+2)
	314	IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
	315	IRNGT (IWRKD+3) = IRNG2
[5159]	316
[3358]	317	! 3 1 2
[5159]	318
[3358]	319	Else
	320	IRNG1 = IRNGT (IWRKD+1)
	321	IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
	322	IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
	323	IRNGT (IWRKD+2) = IRNG1
	324	End If
	325	Exit
	326	End If
[5159]	327
[3358]	328	! 1 2 3 4
[5159]	329
[3358]	330	If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
[5159]	331
[3358]	332	! 1 3 x x
[5159]	333
[3358]	334	If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
	335	IRNG2 = IRNGT (IWRKD+2)
	336	IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
	337	If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
	338	! 1 3 2 4
	339	IRNGT (IWRKD+3) = IRNG2
	340	Else
	341	! 1 3 4 2
	342	IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
	343	IRNGT (IWRKD+4) = IRNG2
	344	End If
[5159]	345
[3358]	346	! 3 x x x
[5159]	347
[3358]	348	Else
	349	IRNG1 = IRNGT (IWRKD+1)
	350	IRNG2 = IRNGT (IWRKD+2)
	351	IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
	352	If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
	353	IRNGT (IWRKD+2) = IRNG1
	354	If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
	355	! 3 1 2 4
	356	IRNGT (IWRKD+3) = IRNG2
	357	Else
	358	! 3 1 4 2
	359	IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
	360	IRNGT (IWRKD+4) = IRNG2
	361	End If
	362	Else
	363	! 3 4 1 2
	364	IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
	365	IRNGT (IWRKD+3) = IRNG1
	366	IRNGT (IWRKD+4) = IRNG2
	367	End If
	368	End If
	369	End Do
[5159]	370
[3358]	371	! The Cs become As and Bs
[5159]	372
[3358]	373	LMTNA = 4
	374	Exit
	375	End Do
[5159]	376
[3358]	377	! Iteration loop. Each time, the length of the ordered subsets
	378	! is doubled.
[5159]	379
[3358]	380	Do
	381	If (LMTNA >= NVAL) Exit
	382	IWRKF = 0
	383	LMTNC = 2 * LMTNC
[5159]	384
[3358]	385	! Loop on merges of A and B into C
[5159]	386
[3358]	387	Do
	388	IWRK = IWRKF
	389	IWRKD = IWRKF + 1
	390	JINDA = IWRKF + LMTNA
	391	IWRKF = IWRKF + LMTNC
	392	If (IWRKF >= NVAL) Then
	393	If (JINDA >= NVAL) Exit
	394	IWRKF = NVAL
	395	End If
	396	IINDA = 1
	397	IINDB = JINDA + 1
[5159]	398
[3358]	399	! Shortcut for the case when the max of A is smaller
	400	! than the min of B. This line may be activated when the
	401	! initial set is already close to sorted.
[5159]	402
[3358]	403	! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
[5159]	404
[3358]	405	! One steps in the C subset, that we build in the final rank array
[5159]	406
[3358]	407	! Make a copy of the rank array for the merge iteration
[5159]	408
[3358]	409	JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
[5159]	410
[3358]	411	XVALA = XDONT (JWRKT(IINDA))
	412	XVALB = XDONT (IRNGT(IINDB))
[5159]	413
[3358]	414	Do
	415	IWRK = IWRK + 1
[5159]	416
[3358]	417	! We still have unprocessed values in both A and B
[5159]	418
[3358]	419	If (XVALA > XVALB) Then
	420	IRNGT (IWRK) = IRNGT (IINDB)
	421	IINDB = IINDB + 1
	422	If (IINDB > IWRKF) Then
	423	! Only A still with unprocessed values
	424	IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
	425	Exit
	426	End If
	427	XVALB = XDONT (IRNGT(IINDB))
	428	Else
	429	IRNGT (IWRK) = JWRKT (IINDA)
	430	IINDA = IINDA + 1
	431	If (IINDA > LMTNA) Exit! Only B still with unprocessed values
	432	XVALA = XDONT (JWRKT(IINDA))
	433	End If
[5159]	434
[3358]	435	End Do
	436	End Do
[5159]	437
[3358]	438	! The Cs become As and Bs
[5159]	439
[3358]	440	LMTNA = 2 * LMTNA
	441	End Do
[5159]	442
[3358]	443	Return
[5159]	444
[3358]	445	End Subroutine R_mrgrnk
	446	Subroutine I_mrgrnk (XDONT, IRNGT)
	447	! __________________________________________________________
	448	! MRGRNK = Merge-sort ranking of an array
	449	! For performance reasons, the first 2 passes are taken
	450	! out of the standard loop, and use dedicated coding.
	451	! __________________________________________________________
	452	! __________________________________________________________
	453	Integer, Dimension (:), Intent (In) :: XDONT
	454	Integer, Dimension (:), Intent (Out) :: IRNGT
	455	! __________________________________________________________
	456	Integer :: XVALA, XVALB
[5159]	457
[3358]	458	Integer, Dimension (SIZE(IRNGT)) :: JWRKT
	459	Integer :: LMTNA, LMTNC, IRNG1, IRNG2
	460	Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
[5159]	461
[3358]	462	NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
	463	Select Case (NVAL)
	464	Case (:0)
	465	Return
	466	Case (1)
	467	IRNGT (1) = 1
	468	Return
[5095]	469	Case Default
	470	Continue
[3358]	471	End Select
[5159]	472
[3358]	473	! Fill-in the index array, creating ordered couples
[5159]	474
[3358]	475	Do IIND = 2, NVAL, 2
	476	If (XDONT(IIND-1) <= XDONT(IIND)) Then
	477	IRNGT (IIND-1) = IIND - 1
	478	IRNGT (IIND) = IIND
	479	Else
	480	IRNGT (IIND-1) = IIND
	481	IRNGT (IIND) = IIND - 1
	482	End If
	483	End Do
	484	If (Modulo(NVAL, 2) /= 0) Then
	485	IRNGT (NVAL) = NVAL
	486	End If
[5159]	487
[3358]	488	! We will now have ordered subsets A - B - A - B - ...
	489	! and merge A and B couples into C - C - ...
[5159]	490
[3358]	491	LMTNA = 2
	492	LMTNC = 4
[5159]	493
[3358]	494	! First iteration. The length of the ordered subsets goes from 2 to 4
[5159]	495
[3358]	496	Do
	497	If (NVAL <= 2) Exit
[5159]	498
[3358]	499	! Loop on merges of A and B into C
[5159]	500
[3358]	501	Do IWRKD = 0, NVAL - 1, 4
	502	If ((IWRKD+4) > NVAL) Then
	503	If ((IWRKD+2) >= NVAL) Exit
[5159]	504
[3358]	505	! 1 2 3
[5159]	506
[3358]	507	If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
[5159]	508
[3358]	509	! 1 3 2
[5159]	510
[3358]	511	If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
	512	IRNG2 = IRNGT (IWRKD+2)
	513	IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
	514	IRNGT (IWRKD+3) = IRNG2
[5159]	515
[3358]	516	! 3 1 2
[5159]	517
[3358]	518	Else
	519	IRNG1 = IRNGT (IWRKD+1)
	520	IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
	521	IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
	522	IRNGT (IWRKD+2) = IRNG1
	523	End If
	524	Exit
	525	End If
[5159]	526
[3358]	527	! 1 2 3 4
[5159]	528
[3358]	529	If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
[5159]	530
[3358]	531	! 1 3 x x
[5159]	532
[3358]	533	If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
	534	IRNG2 = IRNGT (IWRKD+2)
	535	IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
	536	If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
	537	! 1 3 2 4
	538	IRNGT (IWRKD+3) = IRNG2
	539	Else
	540	! 1 3 4 2
	541	IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
	542	IRNGT (IWRKD+4) = IRNG2
	543	End If
[5159]	544
[3358]	545	! 3 x x x
[5159]	546
[3358]	547	Else
	548	IRNG1 = IRNGT (IWRKD+1)
	549	IRNG2 = IRNGT (IWRKD+2)
	550	IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
	551	If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
	552	IRNGT (IWRKD+2) = IRNG1
	553	If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
	554	! 3 1 2 4
	555	IRNGT (IWRKD+3) = IRNG2
	556	Else
	557	! 3 1 4 2
	558	IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
	559	IRNGT (IWRKD+4) = IRNG2
	560	End If
	561	Else
	562	! 3 4 1 2
	563	IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
	564	IRNGT (IWRKD+3) = IRNG1
	565	IRNGT (IWRKD+4) = IRNG2
	566	End If
	567	End If
	568	End Do
[5159]	569
[3358]	570	! The Cs become As and Bs
[5159]	571
[3358]	572	LMTNA = 4
	573	Exit
	574	End Do
[5159]	575
[3358]	576	! Iteration loop. Each time, the length of the ordered subsets
	577	! is doubled.
[5159]	578
[3358]	579	Do
	580	If (LMTNA >= NVAL) Exit
	581	IWRKF = 0
	582	LMTNC = 2 * LMTNC
[5159]	583
[3358]	584	! Loop on merges of A and B into C
[5159]	585
[3358]	586	Do
	587	IWRK = IWRKF
	588	IWRKD = IWRKF + 1
	589	JINDA = IWRKF + LMTNA
	590	IWRKF = IWRKF + LMTNC
	591	If (IWRKF >= NVAL) Then
	592	If (JINDA >= NVAL) Exit
	593	IWRKF = NVAL
	594	End If
	595	IINDA = 1
	596	IINDB = JINDA + 1
[5159]	597
[3358]	598	! Shortcut for the case when the max of A is smaller
	599	! than the min of B. This line may be activated when the
	600	! initial set is already close to sorted.
[5159]	601
[3358]	602	! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
[5159]	603
[3358]	604	! One steps in the C subset, that we build in the final rank array
[5159]	605
[3358]	606	! Make a copy of the rank array for the merge iteration
[5159]	607
[3358]	608	JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
[5159]	609
[3358]	610	XVALA = XDONT (JWRKT(IINDA))
	611	XVALB = XDONT (IRNGT(IINDB))
[5159]	612
[3358]	613	Do
	614	IWRK = IWRK + 1
[5159]	615
[3358]	616	! We still have unprocessed values in both A and B
[5159]	617
[3358]	618	If (XVALA > XVALB) Then
	619	IRNGT (IWRK) = IRNGT (IINDB)
	620	IINDB = IINDB + 1
	621	If (IINDB > IWRKF) Then
	622	! Only A still with unprocessed values
	623	IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
	624	Exit
	625	End If
	626	XVALB = XDONT (IRNGT(IINDB))
	627	Else
	628	IRNGT (IWRK) = JWRKT (IINDA)
	629	IINDA = IINDA + 1
	630	If (IINDA > LMTNA) Exit! Only B still with unprocessed values
	631	XVALA = XDONT (JWRKT(IINDA))
	632	End If
[5159]	633
[3358]	634	End Do
	635	End Do
[5159]	636
[3358]	637	! The Cs become As and Bs
[5159]	638
[3358]	639	LMTNA = 2 * LMTNA
	640	End Do
[5159]	641
[3358]	642	Return
[5159]	643
[3358]	644	End Subroutine I_mrgrnk
	645	end module m_mrgrnk

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format