Merge lp:~maddevelopers/mg5amcnlo/HwU_allocated into lp:~mg5core2/mg5amcnlo/2.5.5
- HwU_allocated
- Merge into 2.5.5
Proposed by
Rikkert Frederix
Status: | Merged |
---|---|
Merge reported by: | Olivier Mattelaer |
Merged at revision: | not available |
Proposed branch: | lp:~maddevelopers/mg5amcnlo/HwU_allocated |
Merge into: | lp:~mg5core2/mg5amcnlo/2.5.5 |
Diff against target: |
449 lines (+202/-83) 3 files modified
Template/NLO/FixedOrderAnalysis/HwU.f (+200/-59) Template/NLO/FixedOrderAnalysis/HwU.inc (+0/-22) Template/NLO/FixedOrderAnalysis/analysis_HwU_template.f (+2/-2) |
To merge this branch: | bzr merge lp:~maddevelopers/mg5amcnlo/HwU_allocated |
Related bugs: |
Reviewer | Review Type | Date Requested | Status |
---|---|---|---|
Valentin Hirschi (community) | Approve | ||
Review via email: mp+323759@code.launchpad.net |
Commit message
Description of the change
Converted a large HwU common block to a module. This allows for dynamic (de-)allocation of the memory required by all variables in the common block, reducing the footprint and size of the executable in case large analyses are used. Moreover, this leads to more flexible plotting options (more histograms can be 'booked' with as many bins as one requires).
Changes are very much contained to (almost) a single file, hence, review should be straight-forward.
To post a comment you must log in.
- 295. By Rikkert Frederix
-
small fix for linux
Preview Diff
[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1 | === modified file 'Template/NLO/FixedOrderAnalysis/HwU.f' |
2 | --- Template/NLO/FixedOrderAnalysis/HwU.f 2017-02-03 13:37:36 +0000 |
3 | +++ Template/NLO/FixedOrderAnalysis/HwU.f 2017-05-09 14:02:00 +0000 |
4 | @@ -1,7 +1,7 @@ |
5 | CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
6 | C C |
7 | C HwU: Histograms with Uncertainties C |
8 | -C By Rikkert Frederix, Dec. 2014 C |
9 | +C By Rikkert Frederix, 12-2014--05-2017 C |
10 | C C |
11 | C Book, fill and write out histograms. Particularly suited for NLO C |
12 | C computations with correlations between points (ie. event and C |
13 | @@ -9,26 +9,41 @@ |
14 | C and PDF uncertainties through reweighting). C |
15 | C C |
16 | CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
17 | + |
18 | +c The module contains effectively the common block with allocatable |
19 | +c variables (something not possible in old fortran version) |
20 | + module HwU_variables |
21 | + implicit none |
22 | + integer :: max_plots,max_points,max_bins,nwgts,np |
23 | + integer :: error_estimation=3 |
24 | + logical, allocatable :: booked(:) |
25 | + integer, allocatable :: nbin(:),histi(:,:),p_bin(:),p_label(:) |
26 | + character(len=50), allocatable :: title(:) |
27 | + character(len=50), allocatable :: wgts_info(:) |
28 | + double precision, allocatable :: histy(:,:,:),histy_acc(:,:,:) |
29 | + $ ,histy2(:,:),histy_err(:,:),histxl(:,:),histxm(:,:) |
30 | + $ ,step(:),p_wgts(:,:) |
31 | + save |
32 | + end module HwU_variables |
33 | + |
34 | |
35 | c To be called once at the start of each run. Initialises the packages |
36 | c and sets the number of weights that need to be included for each point. |
37 | subroutine HwU_inithist(nweights,wgt_info) |
38 | + use HwU_variables |
39 | implicit none |
40 | - include "HwU.inc" |
41 | integer i,nweights |
42 | character*(*) wgt_info(*) |
43 | - do i=1,max_plots |
44 | - booked(i)=.false. |
45 | - enddo |
46 | + call HwU_deallocate_all |
47 | + max_plots=0 |
48 | + max_points=0 |
49 | + max_bins=0 |
50 | + np=0 |
51 | c Number of weights associated to each point. Note that the first |
52 | c weight should always be the 'central value' and it should not be |
53 | c zero if any of the other weights are non-zero. |
54 | nwgts=nweights |
55 | - if (nwgts.gt.max_wgts) then |
56 | - write (*,*) 'ERROR: increase max_wgts in HwU histogramming' |
57 | - $ ,max_wgts,nwgts |
58 | - stop 1 |
59 | - endif |
60 | + allocate(wgts_info(nwgts)) |
61 | do i=1,nwgts |
62 | wgts_info(i)=wgt_info(i) |
63 | enddo |
64 | @@ -48,10 +63,9 @@ |
65 | c square root of the sum of the squares. Perform a weighted average |
66 | c iteration-by-iteration |
67 | c input 3: Same as input 2, but weighted average is same as from MINT |
68 | + use HwU_variables |
69 | implicit none |
70 | integer input |
71 | - integer error_estimation |
72 | - common /HwU_common2/ error_estimation |
73 | if (input.ge.0 .and. input.le.3) then |
74 | error_estimation=input |
75 | else |
76 | @@ -60,39 +74,20 @@ |
77 | endif |
78 | return |
79 | end |
80 | - |
81 | - block data HwU |
82 | -c set the default for the error estimation method. To reduce the size of |
83 | -c the executable, put the error_estimation variable in a separate common |
84 | -c block. If we would have included the 'HwU.inc' file here, that |
85 | -c complete common block seems to be included in the size executable |
86 | -c (approx. 110 MB). |
87 | - integer error_estimation |
88 | - common /HwU_common2/ error_estimation |
89 | - data error_estimation /3/ |
90 | - end |
91 | - |
92 | + |
93 | c Book the histograms at the start of the run. Give a 'label' (an |
94 | c integer) that identifies the plot when filling it and a title |
95 | c ('title_l') for each plot. Also the number of bins ('nbin_l') and the |
96 | c plot range (from 'xmin' to 'xmax') should be given. |
97 | subroutine HwU_book(label,title_l,nbin_l,xmin,xmax) |
98 | + use HwU_variables |
99 | implicit none |
100 | - include "HwU.inc" |
101 | integer label,nbin_l,i,j |
102 | character*(*) title_l |
103 | double precision xmin,xmax |
104 | -c Check that label and number of bins are reasonable |
105 | - if (label.gt.max_plots) then |
106 | - write (*,*) 'ERROR: increase max_plots in HwU histogramming' |
107 | - $ ,max_plots, label |
108 | - stop 1 |
109 | - endif |
110 | - if (nbin_l.gt.max_bins) then |
111 | - write (*,*) 'ERROR: increase max_bins in HwU histogramming' |
112 | - $ ,max_bins,nbin_l |
113 | - stop 1 |
114 | - endif |
115 | +c Allocate space for new histograms if needed |
116 | + call HwU_allocate_histo(label,nbin_l) |
117 | +c Setup the histogram |
118 | booked(label)=.true. |
119 | title(label)=title_l |
120 | nbin(label)=nbin_l |
121 | @@ -111,7 +106,6 @@ |
122 | histy2(label,i)=0d0 |
123 | histy_err(label,i)=0d0 |
124 | enddo |
125 | - np=0 |
126 | return |
127 | end |
128 | |
129 | @@ -121,8 +115,8 @@ |
130 | c the 'HwU_inithist' subroutine. That means that each point should have |
131 | c the same number of weights. |
132 | subroutine HwU_fill(label,x,wgts) |
133 | + use HwU_variables |
134 | implicit none |
135 | - include "HwU.inc" |
136 | integer label,i,j,bin |
137 | double precision x, wgts(*) |
138 | c If central weight is zero do not add this point. |
139 | @@ -147,11 +141,7 @@ |
140 | enddo |
141 | c If a new bin, add it to the list of points |
142 | np=np+1 |
143 | - if (np.gt.max_points) then |
144 | - write (*,*) 'ERROR: increase max_points in HwU histogramming' |
145 | - $ ,max_points |
146 | - stop 1 |
147 | - endif |
148 | + call HwU_allocate_p |
149 | p_label(np)=label |
150 | p_bin(np)=bin |
151 | do j=1,nwgts |
152 | @@ -159,7 +149,7 @@ |
153 | enddo |
154 | return |
155 | end |
156 | - |
157 | + |
158 | c Call after all correlated contributions for a give phase-space |
159 | c point. I.e., every time you get a new set of random numbers from |
160 | c MINT/VEGAS. It adds the current list of points to the histograms. Add |
161 | @@ -168,8 +158,8 @@ |
162 | c this way, correlations between events and counter-events can be |
163 | c correctly taken into account. |
164 | subroutine HwU_add_points |
165 | + use HwU_variables |
166 | implicit none |
167 | - include "HwU.inc" |
168 | integer i,j |
169 | do i=1,np |
170 | do j=1,nwgts |
171 | @@ -192,12 +182,11 @@ |
172 | c the current iteration so that they can be filled with the next |
173 | c iteration. |
174 | subroutine HwU_accum_iter(inclde,nPSpoints,values) |
175 | + use HwU_variables |
176 | implicit none |
177 | - include "HwU.inc" |
178 | logical inclde |
179 | integer nPSpoints,label,i,j |
180 | - double precision nPSinv,etot,vtot(max_wgts),niter,y_squared |
181 | - $ ,values(2) |
182 | + double precision nPSinv,etot,niter,y_squared,values(2) |
183 | data niter /0d0/ |
184 | nPSinv = 1d0/dble(nPSpoints) |
185 | if (inclde) niter = niter+1d0 |
186 | @@ -226,8 +215,8 @@ |
187 | c intermediate stages this function can be called (together with |
188 | c HwU_output) to write intermediate plots to disk. |
189 | subroutine finalize_histograms(nPSpoints) |
190 | + use HwU_variables |
191 | implicit none |
192 | - include "HwU.inc" |
193 | integer label,nPSpoints,i,j |
194 | double precision nPSinv,niter,dummy(2) |
195 | nPSinv=1d0/dble(nPSpoints) |
196 | @@ -257,13 +246,13 @@ |
197 | c histograms the central value should not be zero if any of the other |
198 | c weights are non-zero. |
199 | subroutine accumulate_results(label,nPSinv,niter,values) |
200 | + use HwU_variables |
201 | implicit none |
202 | - include "HwU.inc" |
203 | integer label,i,j |
204 | - double precision nPSinv,etot,vtot(max_wgts),niter,y_squared |
205 | + double precision nPSinv,etot,niter,y_squared |
206 | $ ,values(2),a1,a2 |
207 | - integer error_estimation |
208 | - common /HwU_common2/ error_estimation |
209 | + double precision,allocatable :: vtot(:) |
210 | + if (.not. allocated(vtot)) allocate(vtot(nwgts)) |
211 | if (error_estimation.eq.2) then |
212 | c Use the weighted average bin-by-bin. This is not really justified |
213 | c for fNLO computations, because for bins with low statistics, the |
214 | @@ -395,14 +384,15 @@ |
215 | c Write the histograms to disk at the end of the run, multiplying the |
216 | c output by 'xnorm' |
217 | subroutine HwU_output(unit,xnorm) |
218 | + use HwU_variables |
219 | implicit none |
220 | - include "HwU.inc" |
221 | integer unit,i,j,label |
222 | integer max_length |
223 | - parameter (max_length=(max_wgts+3)*17) |
224 | - character*(max_length) buffer |
225 | + character(len=:), allocatable :: buffer |
226 | character*4 str_nbin |
227 | double precision xnorm |
228 | + if (.not. allocated(buffer)) |
229 | + & allocate(character(len=(nwgts+3)*17) :: buffer) |
230 | c column info: x_min, x_max, y (central value), dy, {extra |
231 | c weights}. |
232 | write (unit,'(a$)') '##& xmin' |
233 | @@ -440,9 +430,157 @@ |
234 | write (unit,'(a)') '' |
235 | write (unit,'(a)') '' |
236 | enddo |
237 | - return |
238 | - end |
239 | - |
240 | + call HwU_deallocate_all |
241 | + deallocate(buffer) |
242 | + return |
243 | + end |
244 | + |
245 | +c Clean all the allocatable variables: |
246 | + subroutine HwU_deallocate_all |
247 | + use HwU_variables |
248 | + implicit none |
249 | + if (allocated(wgts_info)) deallocate(wgts_info) |
250 | + if (allocated(booked)) deallocate(booked) |
251 | + if (allocated(title)) deallocate(title) |
252 | + if (allocated(nbin)) deallocate(nbin) |
253 | + if (allocated(step)) deallocate(step) |
254 | + if (allocated(histxl)) deallocate(histxl) |
255 | + if (allocated(histxm)) deallocate(histxm) |
256 | + if (allocated(histy)) deallocate(histy) |
257 | + if (allocated(histy_acc)) deallocate(histy_acc) |
258 | + if (allocated(histi)) deallocate(histi) |
259 | + if (allocated(histy2)) deallocate(histy2) |
260 | + if (allocated(histy_err)) deallocate(histy_err) |
261 | + if (allocated(p_bin)) deallocate(p_bin) |
262 | + if (allocated(p_label)) deallocate(p_label) |
263 | + if (allocated(p_wgts)) deallocate(p_wgts) |
264 | + return |
265 | + end |
266 | + |
267 | + |
268 | + subroutine HwU_allocate_p |
269 | + use HwU_variables |
270 | + implicit none |
271 | + integer,allocatable :: itemp1(:) |
272 | + double precision, allocatable :: temp2(:,:) |
273 | + if (.not. allocated(p_bin)) then |
274 | + allocate(p_bin(max_plots)) |
275 | + allocate(p_label(max_plots)) |
276 | + allocate(p_wgts(nwgts,max_plots)) |
277 | + max_points=max_plots |
278 | + else |
279 | + if (np.gt.max_points) then |
280 | +c p_bin |
281 | + allocate(itemp1(np+max_plots)) |
282 | + itemp1(1:max_points)=p_bin |
283 | + call move_alloc(itemp1,p_bin) |
284 | + |
285 | +c p_label |
286 | + allocate(itemp1(np+max_plots)) |
287 | + itemp1(1:max_points)=p_label |
288 | + call move_alloc(itemp1,p_label) |
289 | +c p_wgts |
290 | + allocate(temp2(nwgts,np+max_plots)) |
291 | + temp2(1:nwgts,1:max_points)=p_wgts |
292 | + call move_alloc(temp2,p_wgts) |
293 | + max_points=np+max_plots |
294 | + endif |
295 | + endif |
296 | + return |
297 | + end |
298 | + |
299 | + subroutine HwU_allocate_histo(label,nbin_l) |
300 | + use HwU_variables |
301 | + implicit none |
302 | + logical,allocatable :: ltemp(:) |
303 | + integer,allocatable :: itemp1(:),itemp2(:,:) |
304 | + character(len=50),allocatable :: ctemp(:) |
305 | + double precision, allocatable :: temp1(:),temp2(:,:),temp3(:,:,:) |
306 | + integer label,i,nbin_l,label_max,nbin_max |
307 | + logical debug |
308 | + parameter (debug=.true.) |
309 | +c Check if variables are already allocated. If not, simply allocate a |
310 | +c single histogram |
311 | + if (.not. allocated(booked)) then |
312 | + allocate(booked(1)) |
313 | + booked(1)=.false. |
314 | + allocate(title(1)) |
315 | + allocate(nbin(1)) |
316 | + allocate(step(1)) |
317 | + allocate(histxl(1,nbin_l)) |
318 | + allocate(histxm(1,nbin_l)) |
319 | + allocate(histy(nwgts,1,nbin_l)) |
320 | + allocate(histy_acc(nwgts,1,nbin_l)) |
321 | + allocate(histi(1,nbin_l)) |
322 | + allocate(histy2(1,nbin_l)) |
323 | + allocate(histy_err(1,nbin_l)) |
324 | + max_plots=1 |
325 | + max_bins=nbin_l |
326 | + endif |
327 | +c If current label is greater than the plots already allocated, increase |
328 | +c the size of the allocated arrays. This is kind of slow, but shouldn't |
329 | +c really matter since it's only done at the start of a run. |
330 | + if (label.gt.max_plots .or. nbin_l.gt.max_bins) then |
331 | + label_max=max(label,max_plots) |
332 | + nbin_max=max(nbin_l,max_bins) |
333 | +c booked |
334 | + allocate(ltemp(label_max)) |
335 | + ltemp(1:max_plots)=booked |
336 | + call move_alloc(ltemp,booked) |
337 | + do i=max_plots+1,label_max |
338 | + booked(i)=.false. ! histos have not yet been setup |
339 | + enddo |
340 | +c title |
341 | + allocate(ctemp(label_max)) |
342 | + ctemp(1:max_plots)=title |
343 | + call move_alloc(ctemp,title) |
344 | +c nbin |
345 | + allocate(itemp1(label_max)) |
346 | + itemp1(1:max_plots)=nbin |
347 | + call move_alloc(itemp1,nbin) |
348 | +c step |
349 | + allocate(temp1(label_max)) |
350 | + temp1(1:max_plots)=step |
351 | + call move_alloc(temp1,step) |
352 | +c histxl |
353 | + allocate(temp2(label_max,nbin_max)) |
354 | + temp2(1:max_plots,1:max_bins)=histxl |
355 | + call move_alloc(temp2,histxl) |
356 | +c histxm |
357 | + allocate(temp2(label_max,nbin_max)) |
358 | + temp2(1:max_plots,1:max_bins)=histxm |
359 | + call move_alloc(temp2,histxm) |
360 | +c histy |
361 | + allocate(temp3(nwgts,label_max,nbin_max)) |
362 | + temp3(1:nwgts,1:max_plots,1:max_bins)=histy |
363 | + call move_alloc(temp3,histy) |
364 | +c histy_acc |
365 | + allocate(temp3(nwgts,label_max,nbin_max)) |
366 | + temp3(1:nwgts,1:max_plots,1:max_bins)=histy_acc |
367 | + call move_alloc(temp3,histy_acc) |
368 | +c histi |
369 | + allocate(itemp2(label_max,nbin_max)) |
370 | + itemp2(1:max_plots,1:max_bins)=histi |
371 | + call move_alloc(itemp2,histi) |
372 | +c histy2 |
373 | + allocate(temp2(label_max,nbin_max)) |
374 | + temp2(1:max_plots,1:max_bins)=histy2 |
375 | + call move_alloc(temp2,histy2) |
376 | +c histy_err |
377 | + allocate(temp2(label_max,nbin_max)) |
378 | + temp2(1:max_plots,1:max_bins)=histy_err |
379 | + call move_alloc(temp2,histy_err) |
380 | +c Update maximums |
381 | + max_plots=label_max |
382 | + max_bins=nbin_max |
383 | + elseif (booked(label)) then |
384 | + write (*,*) 'ERROR in HwU.f: histogram already booked',label |
385 | + stop |
386 | + endif |
387 | + return |
388 | + end |
389 | + |
390 | + |
391 | c dummy subroutine |
392 | subroutine accum(idummy) |
393 | integer idummy |
394 | @@ -451,3 +589,6 @@ |
395 | subroutine addfil(string) |
396 | character*(*) string |
397 | end |
398 | + |
399 | + |
400 | + |
401 | |
402 | === removed file 'Template/NLO/FixedOrderAnalysis/HwU.inc' |
403 | --- Template/NLO/FixedOrderAnalysis/HwU.inc 2016-02-18 14:05:45 +0000 |
404 | +++ Template/NLO/FixedOrderAnalysis/HwU.inc 1970-01-01 00:00:00 +0000 |
405 | @@ -1,22 +0,0 @@ |
406 | -* -*-fortran-*- |
407 | - |
408 | - integer max_plots,max_bins,max_wgts,max_points |
409 | - parameter (max_plots=200) |
410 | - parameter (max_bins=100) |
411 | - parameter (max_wgts=1024) |
412 | - parameter (max_points=max_plots*40) |
413 | - |
414 | - logical booked(max_plots) |
415 | - integer nbin(max_plots),nwgts,np,p_bin(max_points) |
416 | - & ,p_label(max_points),histi(max_plots,max_bins) |
417 | - character*50 title(max_plots) |
418 | - character*50 wgts_info(max_wgts) |
419 | - double precision histy(max_wgts,max_plots,max_bins) |
420 | - $ ,histy_acc(max_wgts,max_plots,max_bins),histy2(max_plots |
421 | - $ ,max_bins),histy_err(max_plots,max_bins),histxl(max_plots |
422 | - $ ,max_bins),histxm(max_plots,max_bins),step(max_plots) |
423 | - $ ,p_wgts(max_wgts,max_points) |
424 | - |
425 | - common/HwU_common/histy,histy_acc,histy2,histy_err,histxl,histxm |
426 | - & ,p_wgts,step,histi,nbin,p_bin,p_label,np,nwgts |
427 | - & ,booked,title,wgts_info |
428 | |
429 | === modified file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_template.f' |
430 | --- Template/NLO/FixedOrderAnalysis/analysis_HwU_template.f 2014-12-03 11:51:14 +0000 |
431 | +++ Template/NLO/FixedOrderAnalysis/analysis_HwU_template.f 2017-05-09 14:02:00 +0000 |
432 | @@ -17,12 +17,12 @@ |
433 | c o) The first argument is an integer that labels the histogram. In |
434 | c the analysis_end and analysis_fill subroutines this label is used |
435 | c to keep track of the histogram. The label should be a number |
436 | -c between 1 and max_plots=200 (can be increased in HwU.inc). |
437 | +c starting at 1 and be increased for each plot. |
438 | c o) The second argument is a string that will apear above the |
439 | c histogram. Do not use brackets "(" or ")" inside this string. |
440 | c o) The third, forth and fifth arguments are the number of bis, the |
441 | c lower edge of the first bin and the upper edge of the last |
442 | -c bin. There is a maximum of 100 bins per histogram. |
443 | +c bin. |
444 | c o) When including scale and/or PDF uncertainties, declare a |
445 | c histogram for each weight, and compute the uncertainties from the |
446 | c final set of histograms |
447 | |
448 | === removed symlink 'Template/NLO/MCatNLO/include/HwU.inc' |
449 | === target was u'../../FixedOrderAnalysis/HwU.inc' |
Hi Rik,
Sorry for the late review.
The code seems OK and definitely a welcome addition.
It seems that all arrays are correctly deallocated in HwU_deallocate_all, but I see that it is only called in the "output" subroutine, which is in principle not ideal.
The usual structure is to have a master "clean_up" function in MadFKS that is called at the end of any driver and takes care of calling the deallocation of all possible modules used.
Now in this case, it of course works since we know that the user will output once only and at the end of the run basically.
I think move_alloc is supported in gfortran 4.6 (our min. requirement), but it is always worth checking that it compiles with v4.6 when coding beyond f77.
Also, ideally, one could turn the whole HwU.f into a module (not just its variables) and use it as library, putting it in the Source directory, but it's not crucial (and would also mean adding 'use HwU' in all the existing FO analysis which is not optimal).
Also note that in fortran2003+, one can use "derived types" which are very useful to mimic C++ 'struct' in fortran.
This could be especially useful in cases like this one where one would really want to allocate and de-allocate arrays of histograms directly, which could themselves contain allocatable arrays for their various properties (such as bins).
I had done something similar for the the dimensions of my DiscreteSampler (where bins are themselves derived types):
type sampledDimension threshold probing_ points
integer :: grid_mode
real*8 :: small_contrib_
real*8 :: damping_power
integer :: min_bin_
real*8 :: norm
real*8 :: abs_norm
real*8 :: variance_norm
real*8 :: norm_sqr
integer :: n_tot_entries
character, dimension(:), allocatable :: dimension_name
type(bin) , dimension(:), allocatable :: bins
endtype sampledDimension
Anyway, none of this is essential here, so thanks for this improvement!