Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This file contains the implementations of the generic interfaces in [pm_sampling]@(ref pm_sampling).
19 : !>
20 : !> \test
21 : !> [test_pm_sampling](@ref test_pm_sampling)
22 : !>
23 : !> \finmain
24 : !>
25 : !> \author
26 : !> \AmirShahmoradi, September 1, 2012, 12:00 AM, National Institute for Fusion Studies, The University of Texas at Austin
27 :
28 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
29 :
30 : #define RETURN_IF_FAILED(LINE,FAILED,MSG) \
31 : if (FAILED) then; \
32 : err%occurred = .true._LK; \
33 : err%msg = PROCEDURE_NAME//getLine(LINE)//SKC_": "//NL1//MSG; \
34 : call spec%disp%stop%show(err%msg); \
35 : return; \
36 : end if;
37 : ! Define intel-specific shared property of files on Windows.
38 : #if INTEL_ENABLED && WINDOWS_ENABLED
39 : #define SHARED, shared
40 : #else
41 : #define SHARED
42 : #endif
43 : ! Define runtime error handling mode.
44 : #if TESTING_ENABLED && (CAF_ENABLED || MPI_ENABLED || OMP_ENABLED)
45 : #define FAILED_IMAGE(FAILED)isFailedImage(FAILED)
46 : #else
47 : #define FAILED_IMAGE(FAILED)FAILED
48 : #endif
49 : #if CAF_ENABLED || MPI_ENABLED
50 : #define SET_CAFMPI(X)X
51 : #define SET_SERIAL(X)
52 : #else
53 : #define SET_CAFMPI(X)
54 : #define SET_SERIAL(X)X
55 : #endif
56 : #if CAF_ENABLED
57 : #define SET_CAF(X)X
58 : #else
59 : #define SET_CAF(X)
60 : #endif
61 : #if MPI_ENABLED
62 : #define SET_MPI(X)X
63 : #else
64 : #define SET_MPI(X)
65 : #endif
66 : #if OMP_ENABLED
67 : #define SET_OMP(X)X
68 : #else
69 : #define SET_OMP(X)
70 : #endif
71 : #if DEBUG_ENABLED
72 : #define SET_DEBUG(X)X
73 : #else
74 : #define SET_DEBUG(X)
75 : #endif
76 : !#if CAF_ENABLED
77 : !#define COINDEX(I)[I]
78 : !#else
79 : !#define COINDEX(I)
80 : !#endif
81 : ! Set the traceback information.
82 : !#define TRACE(LINE) SKC_"@file::"//getStr(__FILE__)//getLine(LINE)//PROCEDURE_NAME
83 : ! ! Set the alert optional input arguments.
84 : !#define ALERT_OPTIONS unit = reportFileUnit, width = disp%width, prefix = method, tmsize = tmsize, bmsize = bmsize
85 :
86 : !%%%%%%%%%%%%%%%%%%
87 : #if runParaDRAM_ENABLED
88 : !%%%%%%%%%%%%%%%%%%
89 :
90 : abstract interface
91 : #if OMP_ENABLED && (MATLAB_ENABLED || PYTHON_ENABLED || R_ENABLED)
92 : recursive function getLogFuncPtr_proc(logFuncState, ndim, njob, avgTimePerFunCall, avgCommPerFunCall) result(mold) bind(C)
93 : import :: IK, RKC, RKD
94 : integer(IK), intent(in), value :: ndim, njob
95 : real(RKD), intent(inout) :: avgTimePerFunCall, avgCommPerFunCall
96 : real(RKC), intent(inout) :: logFuncState(ndim, njob)
97 : real(RKC) :: mold
98 : end function
99 : #else
100 : recursive function getLogFuncPtr_proc(state, ndim) result(logFunc) bind(C)
101 : import :: IK, RKC
102 : integer(IK), intent(in), value :: ndim
103 : real(RKC), intent(in) :: state(ndim)
104 : real(RKC) :: logFunc
105 : end function
106 : #endif
107 : end interface
108 : procedure(getLogFuncPtr_proc), pointer :: getLogFuncPtr
109 : !procedure(real(RKC)), pointer :: getLogFuncPtr
110 1 : type(paradram_type) :: sampler
111 1 : type(err_type) :: err
112 : integer(IK) :: lenin
113 : stat = 0_IK
114 :
115 : ! Reconstruct the input file.
116 :
117 1 : if (present(input)) then
118 : lenin = 1
119 97 : do
120 98 : if (input(lenin) == c_null_char) exit
121 97 : lenin = lenin + 1
122 : end do
123 1 : if (1 < lenin) sampler%inputFile = getCharSeq(input(1 : lenin - 1))
124 : end if
125 :
126 : !if (0 < lenin) then
127 : ! allocate(character(lenin) :: sampler%inputFile)
128 : ! do iell = 1, lenin
129 : ! sampler%inputFile(iell : iell) = input(iell)
130 : ! end do
131 : !end if
132 :
133 : ! Associate the input C procedure pointer to a Fortran procedure pointer.
134 :
135 1 : call c_f_procpointer(cptr = getLogFunc, fptr = getLogFuncPtr)
136 1 : err = getErrSampling(sampler, getLogFuncWrapper, ndim) ! run ParaDRAM.
137 1 : if (err%occurred) stat = 1_IK
138 1 : nullify(getLogFuncPtr)
139 :
140 : contains
141 :
142 : #if OMP_ENABLED && (MATLAB_ENABLED || PYTHON_ENABLED || R_ENABLED)
143 : recursive function getLogFuncWrapper(logFuncState, avgTimePerFunCall, avgCommPerFunCall) result(mold)
144 : real(RKC), intent(inout), contiguous :: logFuncState(0:,:)
145 : real(RKD), intent(inout) :: avgTimePerFunCall, avgCommPerFunCall
146 : real(RKC) :: mold
147 : mold = getLogFuncPtr(logFuncState, ndim, size(logFuncState, 2, IK), avgTimePerFunCall, avgCommPerFunCall)
148 : end function
149 : #else
150 91228 : recursive function getLogFuncWrapper(state) result(logFunc)
151 : real(RKC), intent(in), contiguous :: state(:)
152 : real(RKC) :: logFunc
153 91228 : logFunc = getLogFuncPtr(state, ndim)
154 : !write(*,"(1I5,5E30.20,:,', ')") precision(state), state!, logFunc
155 91228 : end function
156 : #endif
157 :
158 : !%%%%%%%%%%%%%%%%%%%%%
159 : #elif getErrSampling_ENABLED
160 : !%%%%%%%%%%%%%%%%%%%%%
161 :
162 : use pm_matrixClass, only: isMatClass, posdefmat
163 : use pm_sampling_kernel, only: getErrChainRead, getErrChainWrite, isFailedChainResize
164 : use pm_sampling_kernel, only: spec_type, stat_type, proposal_type
165 : use pm_sampling_kernel, only: chainFileColName
166 : use pm_sampling_kernel, only: getErrKernelRun
167 : use pm_sampling_kernel, only: NL1, NL2, RKC
168 : #if ParaDISE_ENABLED
169 : #define SET_DRAMDISE(X)X
170 : #define SET_ParaNest(X)
171 : character(*,SKC), parameter :: method = SKC_"ParaDISE"
172 : #elif ParaDRAM_ENABLED
173 : #define SET_DRAMDISE(X)X
174 : #define SET_ParaNest(X)
175 : character(*,SKC), parameter :: method = SKC_"ParaDRAM"
176 : #elif ParaNest_ENABLED
177 : #define SET_ParaNest(X)X
178 : #define SET_DRAMDISE(X)
179 : character(*,SKC), parameter :: method = SKC_"ParaNest"
180 : #else
181 : #error "Unrecognized interface."
182 : #endif
183 : character(*,SKC), parameter :: PROCEDURE_NAME = MODULE_NAME//SKC_"@getErrSampling()"
184 : !character(*,SKC), parameter :: NL1 = new_line(SKC_"a"), NL2 = NL1//NL1
185 : integer(IK) :: ndimp1, idim, iq
186 13 : type(proposal_type) :: proposal
187 65 : type(spec_type) :: spec
188 481 : type(stat_type) :: stat
189 13 : ndimp1 = ndim + 1_IK
190 : !!#if CAF_ENABLED
191 : !!! This must be coarray allocatable `[:]`, even though it will have only one element.
192 : !!! Otherwise, GNU Fortran 10 compiler results in runtime segmentation
193 : !!! fault in coarray mode, when the routine is called multiple times.
194 : !!! This used to be the case for a compound object. Does it also hold for intrinsic type objects?
195 : !!integer(IK), allocatable :: comv_unifRandState(:,:)[:]
196 : !!#else
197 : !!integer(IK), allocatable :: comv_unifRandState(:,:)
198 : !!#endif
199 : !
200 : !! Initialize the simulation auxiliary variables.
201 : !
202 : !!allocate(character(2**13 - 1, SK) :: errmsg) ! 8191: roughly 66Kb of memory for error message accumulation.
203 : !!reportFileUnit = output_unit ! Temporarily set the report file to stdout until the report file is set up.
204 :
205 : ! Setup the simulation specifications.
206 :
207 13 : spec = spec_type(modelr_type(0._RKC), method, ndim)
208 13 : err = spec%set(sampler)
209 13 : RETURN_IF_FAILED(__LINE__,FAILED_IMAGE(err%occurred),err%msg)
210 :
211 : !print *, """"//getStr([css_type(chainFileColName), css_type(spec%domainAxisName%val)], format = spec%chainFile%format%header)//achar(0, SKC)//""""
212 13 : if (spec%run%is%new .and. spec%image%is%leader) then
213 13 : if (spec%outputChainFileFormat%isBinary) then
214 0 : write(spec%chainFile%unit) getStr([css_type(chainFileColName), css_type(spec%domainAxisName%val)], format = spec%chainFile%format%header)//achar(0, SKC)
215 : else
216 144 : write(spec%chainFile%unit, spec%chainFile%format%header) (trim(chainFileColName(idim)), idim = 1, size(chainFileColName, 1, IK)), (trim(spec%domainAxisName%val(idim)), idim = 1, ndim)
217 : end if
218 : end if
219 :
220 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
221 :
222 13 : proposal = proposal_type(spec, err)
223 13 : RETURN_IF_FAILED(__LINE__,FAILED_IMAGE(err%occurred),err%msg)
224 :
225 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
226 :
227 13 : if (spec%image%is%leader) call spec%disp%text%wrap(NL1//spec%method%val//SKC_" simulation began on "//getDateTime(SK_"%c %z")//NL1)
228 13 : err = getErrKernelRun(getLogFunc, spec, stat, proposal)
229 13 : RETURN_IF_FAILED(__LINE__,FAILED_IMAGE(err%occurred),err%msg)
230 13 : if (spec%image%is%leader) call spec%disp%text%wrap(NL1//spec%method%val//SKC_" simulation ended on "//getDateTime(SK_"%c %z")//NL1)
231 13 : if (spec%image%is%first) then
232 13 : call spec%disp%skip(count = 2_IK, unit = output_unit)
233 : !call execute_command_line(" ")
234 13 : flush(output_unit)
235 : end if
236 13 : RETURN_IF_FAILED(__LINE__,FAILED_IMAGE(err%occurred),err%msg)
237 :
238 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
239 :
240 13 : blockLeaderPostProcessing: if (spec%image%is%leader) then
241 :
242 : ! Rules:
243 : ! - Statistics values that require more than 1 line of the output file
244 : ! must be in the form of a table, the first line always being the column names.
245 : ! If the subsequent lines have one more column than the header line,
246 : ! then the first column of the subsequent lines will be interpreted as row names.
247 : ! Respecting this rule is important for parsing the contents of the report file in dynamic programming languages.
248 :
249 : associate(format => spec%reportFile%format%generic)
250 :
251 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
252 :
253 13 : call spec%disp%show("stats.numFuncCall.accepted")
254 13 : call spec%disp%show(stat%numFunCallAccepted, format = format)
255 13 : call spec%disp%note%show("This is the total number of accepted function calls (unique samples).")
256 :
257 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
258 :
259 13 : call spec%disp%show("stats.numFuncCall.acceptedRejected")
260 13 : call spec%disp%show(stat%numFunCallAcceptedRejected, format = format)
261 13 : call spec%disp%note%show("This is the total number of accepted or rejected function calls.")
262 :
263 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
264 :
265 13 : call spec%disp%show("stats.numFuncCall.acceptedRejectedDelayed")
266 13 : call spec%disp%show(stat%numFunCallAcceptedRejectedDelayed, format = format)
267 13 : call spec%disp%note%show("This is the total number of accepted or rejected or delayed-rejection (if any requested) function calls.")
268 :
269 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
270 :
271 13 : call spec%disp%show("stats.numFuncCall.acceptedRejectedDelayedUnused")
272 13 : call spec%disp%show(stat%numFunCallAcceptedRejectedDelayedUnused, format = format)
273 13 : call spec%disp%note%show("This is the total number of accepted or rejected or unused function calls (by all processes, including delayed rejections, if any requested).")
274 :
275 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
276 :
277 13 : SET_DRAMDISE(call spec%disp%show("stats.chain.verbose.efficiency.meanAcceptanceRate"))
278 13 : SET_DRAMDISE(call spec%disp%show(stat%cfc%meanAcceptanceRate(stat%numFunCallAccepted), format = format))
279 13 : SET_DRAMDISE(call spec%disp%note%show(SKC_"This is the average MCMC acceptance rate of the "//spec%method%val//SKC_" sampler."))
280 :
281 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
282 :
283 13 : SET_DRAMDISE(call spec%disp%show("stats.chain.verbose.efficiency.acceptedOverAcceptedRejected"))
284 : SET_ParaNest(call spec%disp%show("stats.chain.uniques.efficiency.acceptedOverAcceptedRejected"))
285 13 : call spec%disp%show(real(stat%numFunCallAccepted, RKC) / real(stat%numFunCallAcceptedRejected, RKC), format = format) ! accepted2AcceptedRejected
286 : spec%msg = SKC_"This is the "//spec%method%val//SKC_" sampler efficiency given the accepted and rejected function calls, &
287 13 : &that is, the number of accepted function calls divided by the number of (accepted + rejected) function calls."
288 13 : call spec%disp%note%show(spec%msg)
289 :
290 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
291 :
292 13 : SET_DRAMDISE(call spec%disp%show("stats.chain.verbose.efficiency.acceptedOverAcceptedRejectedDelayed"))
293 : SET_ParaNest(call spec%disp%show("stats.chain.uniques.efficiency.acceptedOverAcceptedRejectedDelayed"))
294 13 : call spec%disp%show(real(stat%numFunCallAccepted, RKC) / real(stat%numFunCallAcceptedRejectedDelayed, RKC), format = format)
295 : spec%msg = SKC_"This is the "//spec%method%val//SKC_" sampler efficiency given the accepted, rejected, and delayed-rejection (if any requested) &
296 13 : &function calls, that is, the number of accepted function calls divided by the number of (accepted + rejected + delayed-rejection) function calls."
297 13 : call spec%disp%note%show(spec%msg)
298 :
299 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
300 :
301 13 : SET_DRAMDISE(call spec%disp%show("stats.chain.verbose.efficiency.acceptedOverAcceptedRejectedDelayedUnused"))
302 : SET_ParaNest(call spec%disp%show("stats.chain.uniques.efficiency.acceptedOverAcceptedRejectedDelayedUnused"))
303 13 : call spec%disp%show(real(stat%numFunCallAccepted, RKC) / real(stat%numFunCallAcceptedRejectedDelayedUnused, RKC), format = format)
304 : spec%msg = SKC_"This is the "//spec%method%val//SKC_" sampler efficiency given the accepted, rejected, delayed-rejection (if any requested), and unused function calls &
305 13 : &(in parallel simulations), that is, the number of accepted function calls divided by the number of (accepted + rejected + delayed-rejection + unused) function calls."
306 13 : call spec%disp%note%show(spec%msg)
307 :
308 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
309 :
310 13 : call spec%disp%show("stats.time.total")
311 13 : call spec%disp%show(stat%timer%clock, format = format)
312 13 : call spec%disp%note%show("This is the total runtime in seconds.")
313 :
314 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
315 :
316 13 : call spec%disp%show("stats.time.perFuncCallAccepted")
317 13 : call spec%disp%show(stat%timer%clock / real(stat%numFunCallAccepted, RKC), format = format)
318 13 : call spec%disp%note%show("This is the average effective time cost of each accepted function call, in seconds.")
319 :
320 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
321 :
322 13 : call spec%disp%show("stats.time.perFuncCallAcceptedRejected")
323 13 : call spec%disp%show(stat%timer%clock / real(stat%numFunCallAcceptedRejected, RKC), format = format)
324 13 : call spec%disp%note%show("This is the average effective time cost of each accepted or rejected function call, in seconds.")
325 :
326 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
327 :
328 13 : call spec%disp%show("stats.time.perFuncCallAcceptedRejectedDelayed")
329 13 : call spec%disp%show(stat%timer%clock / real(stat%numFunCallAcceptedRejectedDelayed, RKC), format = format)
330 13 : call spec%disp%note%show("This is the average effective time cost of each accepted or rejected function call (including delayed-rejections, if any requested), in seconds.")
331 :
332 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
333 :
334 13 : call spec%disp%show("stats.time.perFuncCallAcceptedRejectedDelayedUnused")
335 13 : call spec%disp%show(stat%timer%clock / real(stat%numFunCallAcceptedRejectedDelayedUnused, RKC), format = format)
336 : spec%msg = "This is the average effective time cost of each accepted or rejected or unused function call (including delayed-rejections, if any requested), in seconds. &
337 : &This timing can be directly compared to the corresponding timing of other parallel runs of the same simulation but with different processor counts to assess the parallel scaling. &
338 13 : &A lower timing value implies a better parallel scaling and performance."
339 13 : call spec%disp%note%show(spec%msg)
340 :
341 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
342 :
343 13 : call spec%disp%show("stats.time.perInterProcessCommunication")
344 13 : call spec%disp%show(stat%avgCommPerFunCall, format = format)
345 13 : call spec%disp%note%show("This is the average time cost of parallel inter-process communications per used (accepted or rejected or delayed-rejection) function call, in seconds.")
346 :
347 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
348 :
349 13 : call spec%disp%show("stats.time.perFuncCall")
350 13 : call spec%disp%show(stat%avgTimePerFunCall, format = format)
351 13 : call spec%disp%note%show("This is the average pure time cost of each function call, in seconds.")
352 :
353 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
354 :
355 13 : call spec%disp%show("stats.parallelism.current.process.count")
356 13 : call spec%disp%show(spec%image%count, format = format)
357 13 : call spec%disp%note%show("This is the number of images/processes/threads used in this simulation.")
358 :
359 : #if CAF_ENABLED || MPI_ENABLED || OMP_ENABLED
360 : if (spec%image%count == 1_IK) then
361 : spec%msg = spec%method%val//SKC_" is used in parallel mode with only one processor. This can be computationally inefficient. &
362 : &Consider using the serial version of the code or provide more processes at runtime if it seems to be beneficial."
363 : call spec%disp%note%show(spec%msg, unit = output_unit)!, format = format, tmsize = 3_IK, bmsize = 1_IK)
364 : end if
365 : #endif
366 :
367 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
368 : ! Find individual image contributions and the Cyclic Geometric fit to the contributions.
369 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
370 :
371 : imageContribution_block: block
372 :
373 : real(RKC) :: successProbNormFac(2)
374 : integer(IK) :: pidSuccessLen, iell
375 : integer(IK), allocatable :: index(:), cntSuccess(:), pidSuccess(:)
376 :
377 13 : call setResized(index, spec%image%count)
378 13 : call setResized(cntSuccess, spec%image%count)
379 13 : call setResized(pidSuccess, spec%image%count)
380 13 : call setUnique(stat%cfc%processID, unique = pidSuccess, lenUnique = pidSuccessLen, count = cntSuccess)
381 :
382 13 : if (pidSuccessLen < spec%image%count) then
383 0 : pidSuccess(pidSuccessLen + 1 :) = getComplementRange(pidSuccess(1 : pidSuccessLen), start = 1_IK, stop = spec%image%count, step = 1_IK)
384 0 : cntSuccess(pidSuccessLen + 1 :) = 0_IK
385 : end if
386 13 : call setSorted(pidSuccess, index)
387 39 : pidSuccess(:) = pidSuccess(index)
388 39 : cntSuccess(:) = cntSuccess(index)
389 :
390 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
391 :
392 13 : call spec%disp%show("stats.parallelism.current.process.contribution.count")
393 65 : call spec%disp%show(css_type([character(15) :: "processID", "count"]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
394 26 : do iell = 1, spec%image%count
395 26 : write(spec%reportFile%unit, spec%reportFile%format%integer) iell, cntSuccess(iell)
396 : end do
397 13 : call spec%disp%skip(count = spec%disp%bmsize)
398 : spec%msg = SKC_"These are the contributions of individual processes to the construction of the output chain of the "//spec%method%val//SKC_" sampler. &
399 : &Essentially, they represent the total number of accepted states (useful contributions to the simulation) by the corresponding processor, &
400 13 : &starting from the first processor to the last. This information is mostly informative in parallel Fork-Join (singleChain) simulations."
401 13 : call spec%disp%note%show(spec%msg)
402 :
403 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
404 :
405 13 : call spec%disp%show("stats.parallelism.current.process.contribution.fit")
406 65 : call spec%disp%show(css_type([character(15) :: "successProb", "normFac"]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
407 13 : if (spec%parallelism%is%forkJoin) then
408 0 : pidSuccess = pack(pidSuccess, 0 < cntSuccess)
409 0 : cntSuccess = pack(cntSuccess, 0 < cntSuccess)
410 0 : err%occurred = isFailedGeomCyclicFit(pidSuccess, cntSuccess, spec%image%count, successProbNormFac(1), successProbNormFac(2))
411 0 : if (err%occurred) then
412 0 : call spec%disp%show(css_type([character(15) :: "NaN", "NaN"]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
413 0 : err%occurred = .false._LK
414 : else
415 0 : call spec%disp%show(successProbNormFac, format = spec%reportFile%format%allreal, tmsize = 0_IK)
416 : end if
417 : else
418 39 : successProbNormFac = [1._RKC, real(cntSuccess(size(cntSuccess, 1, IK)), RKC)]
419 13 : call spec%disp%show(successProbNormFac, format = spec%reportFile%format%allreal, tmsize = 0_IK)
420 : end if
421 13 : call spec%disp%skip(count = spec%disp%bmsize)
422 : spec%msg = SKC_"These are the two parameters of the Cyclic Geometric fit to the distribution of the processor contributions to the construction &
423 : &of the final output chain of visited states. The processor contributions are reported in the first column of the output chain file. &
424 : &The processor contribution frequencies are listed above. The fit has the following form: "//NL2// &
425 : SKC_" processConstribution(i) ="//NL1//&
426 : SKC_"successProb * normFac * (1 - successProb)^(i - 1) / (1 - (1 - successProb)^processCount)"//NL2// &
427 : SKC_"where `i` is the ID of the processor (starting from index `1`), `processCount` is the total number of &
428 : &processes used in the simulation, `successProb` is equivalent to an effective sampling efficiency computed &
429 13 : &from the contributions of individual processes to the output chain, and `normFac` is a normalization constant."
430 13 : call spec%disp%note%show(spec%msg)
431 :
432 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
433 :
434 : end block imageContribution_block
435 :
436 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
437 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% begin speedup compute %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
438 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
439 :
440 : blockParallelSpeedup: block
441 :
442 : integer(IK) :: scalingMaxLoc, iell
443 : integer(IK), parameter :: nscol = 5_IK
444 : integer(IK), allocatable :: numproc(:)
445 : real(RKC), allocatable :: scaling(:)
446 : real(RKC) :: scalingMaxVal
447 :
448 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
449 : ! Compute the effective efficiency from the processor contributions and the current strong scaling.
450 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
451 :
452 : ! \todo
453 : ! The specified epsilon values for `parSecTime` and `comSecTime` are points of weakness of the following call
454 : ! if the computers of the future civilization are roughly 1 billion times or more faster than the current 2020 technologies.
455 : ! Therefore, a more robust solution is required for cases where the entire simulation is a dry run of the old existing simulation.
456 : ! This situation is, however, such a rare occurrence that does not merit further investment in the current version of the library.
457 : call setForkJoinScaling ( conProb = real(stat%numFunCallAccepted, RKC) / real(stat%numFunCallAcceptedRejected, RKC) & ! LCOV_EXCL_LINE
458 : , seqSecTime = epsilon(1._RKC) & ! LCOV_EXCL_LINE time cost of the sequential section of the code, which is negligible here
459 : , comSecTime = real(stat%avgCommPerFunCall, RKC) / spec%image%count & ! LCOV_EXCL_LINE
460 : , parSecTime = real(stat%avgTimePerFunCall, RKC) & ! LCOV_EXCL_LINE
461 : , scalingMinLen = max(10_IK, spec%image%count * 2_IK) & ! LCOV_EXCL_LINE
462 : , scalingMaxLoc = scalingMaxLoc & ! LCOV_EXCL_LINE
463 : , scalingMaxVal = scalingMaxVal & ! LCOV_EXCL_LINE
464 : , numproc = numproc & ! LCOV_EXCL_LINE
465 : , scaling = scaling & ! LCOV_EXCL_LINE
466 13 : )
467 :
468 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
469 :
470 13 : call spec%disp%show("stats.parallelism.current.speedup")
471 13 : call spec%disp%show(scaling(spec%image%count), format = format)
472 : spec%msg = "This is the estimated maximum speedup gained via "//spec%parallelism%val// & ! LCOV_EXCL_LINE
473 13 : SKC_" parallelization model compared to serial mode given the current parallel communication overhead."
474 13 : call spec%disp%note%show(spec%msg)
475 :
476 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
477 :
478 13 : call spec%disp%show("stats.parallelism.optimal.process.count")
479 13 : call spec%disp%show(scalingMaxLoc, format = format)
480 : spec%msg = SKC_"This is the predicted optimal number of physical computing processes for "//spec%parallelism%val// & ! LCOV_EXCL_LINE
481 13 : SKC_" parallelization model, given the current sampling efficiency and parallel communication overhead."
482 13 : call spec%disp%note%show(spec%msg)
483 :
484 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
485 :
486 13 : spec%msg = ""
487 13 : call spec%disp%show("stats.parallelism.optimal.speedup")
488 13 : call spec%disp%show(scalingMaxVal, format = format)
489 13 : if (spec%parallelism%is%forkJoin .and. scaling(spec%image%count) < 1._RKC) then
490 : spec%msg = "The time cost of calling the user-provided objective function must be at least "//getStr(1._RKC / scaling(spec%image%count), SK_"(g0.6)")//&
491 : SKC_" times more (that is, ~"//getStr(10**6 * stat%avgTimePerFunCall / scaling(spec%image%count), SK_"(g0.6)")//" microseconds) to see any performance benefits from "//&
492 0 : spec%parallelism%val//SKC_" parallelization model for this simulation. "
493 0 : if (scalingMaxLoc == 1_IK) then
494 0 : spec%msg = spec%msg//SKC_"Consider simulating in serial mode in the future (if used within the same computing environment and with the same configuration as used here)."
495 : else
496 : spec%msg = spec%msg//SKC_"Consider simulating on "//getStr(scalingMaxLoc)//&
497 0 : SKC_" processors in the future (if used within the same computing environment and with the same configuration as used here)."
498 : end if
499 0 : if (.not. spec%outputSplashMode%is%silent) call spec%disp%note%show(spec%msg, unit = output_unit, tmsize = spec%disp%tmsize, bmsize = spec%disp%bmsize)
500 0 : spec%msg = NL1//spec%msg
501 : end if
502 : spec%msg = SKC_"This is the predicted optimal maximum speedup gained via "//spec%parallelism%val// & ! LCOV_EXCL_LINE
503 13 : SKC_" parallelization model, given the current sampling efficiency and parallel communication overhead."//spec%msg
504 13 : call spec%disp%note%show(spec%msg)
505 :
506 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
507 :
508 13 : call spec%disp%show("stats.parallelism.optimal.scaling.strong.speedup")
509 65 : call spec%disp%show(css_type([character(15) :: "processCount", "speedup"]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
510 416 : do iell = 1, size(scaling, 1, IK)
511 416 : write(spec%reportFile%unit, spec%reportFile%format%intreal) numproc(iell), scaling(iell)
512 : !call spec%disp%show(scaling(iell : min(iell + nscol - 1_IK, size(scaling, 1, IK))), format = scalingFormat, tmsize = 0_IK, bmsize = 0_IK)
513 : end do
514 13 : call spec%disp%skip(count = spec%disp%bmsize)
515 : spec%msg = SKC_"This is the predicted strong-scaling speedup behavior of the "//spec%parallelism%val//SKC_" parallelization model, &
516 13 : &given the current sampling efficiency and parallel communication overhead, for increasing numbers of processes, starting from a single process."
517 13 : call spec%disp%note%show(spec%msg)
518 :
519 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
520 : ! compute the absolute parallelism efficiency under any sampling efficiency.
521 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
522 :
523 : ! \todo
524 : ! The specified epsilon values for `parSecTime` and `comSecTime` are points of weakness of the following call
525 : ! if the computers of the future civilization are roughly 1 billion times or more faster than the current 2020 technologies.
526 : ! Therefore, a more robust solution is required for cases where the entire simulation is a dry run of the old existing simulation.
527 : ! This situation is, however, such a rare occurrence that does not merit further investment in the current version of the library.
528 : call setForkJoinScaling ( conProb = 0._RKC & ! LCOV_EXCL_LINE
529 : , seqSecTime = epsilon(1._RKC) & ! LCOV_EXCL_LINE time cost of the sequential section of the code, which is negligible here
530 : , comSecTime = real(stat%avgCommPerFunCall, RKC) / spec%image%count & ! LCOV_EXCL_LINE
531 : , parSecTime = real(stat%avgTimePerFunCall, RKC) & ! LCOV_EXCL_LINE
532 : , scalingMinLen = max(10_IK, spec%image%count * 2_IK) & ! LCOV_EXCL_LINE
533 : , scalingMaxVal = scalingMaxVal & ! LCOV_EXCL_LINE
534 : , scalingMaxLoc = scalingMaxLoc & ! LCOV_EXCL_LINE
535 : , numproc = numproc & ! LCOV_EXCL_LINE
536 : , scaling = scaling & ! LCOV_EXCL_LINE
537 13 : )
538 :
539 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
540 :
541 13 : call spec%disp%show("stats.parallelism.absolute.process.count")
542 13 : call spec%disp%show(scalingMaxLoc, format = format)
543 : spec%msg = "This is the predicted absolute number of physical computing processes for "//spec%parallelism%val// & ! LCOV_EXCL_LINE
544 13 : SKC_" parallelization model, under any sampling efficiency for this sampling problem, given the current parallel communication overhead."
545 13 : call spec%disp%note%show(spec%msg)
546 :
547 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
548 :
549 13 : call spec%disp%show("stats.parallelism.absolute.speedup")
550 13 : call spec%disp%show(scalingMaxVal, format = format)
551 : spec%msg = SKC_"This is the predicted absolute optimal maximum speedup gained via "//spec%parallelism%val//SKC_" parallelization model, under any sampling efficiency. &
552 : &This simulation will likely NOT benefit from any additional computing processors beyond the predicted absolute optimal number, "//getStr(scalingMaxLoc)//SKC_", in the above. &
553 : &This is true for any value of MCMC sampling efficiency given the current parallel communication overhead. Keep in mind that the predicted absolute optimal number of processors &
554 : &is just an estimate whose accuracy depends on many runtime factors, including the topology of the communication network being used, the number of processors per node, &
555 13 : &and the number of tasks to each processor or node."
556 13 : call spec%disp%note%show(spec%msg)
557 :
558 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
559 :
560 13 : call spec%disp%show("stats.parallelism.absolute.scaling.strong.speedup")
561 65 : call spec%disp%show(css_type([character(15) :: "processCount", "speedup"]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
562 416 : do iell = 1, size(scaling, 1, IK)
563 416 : write(spec%reportFile%unit, spec%reportFile%format%intreal) numproc(iell), scaling(iell)
564 : !call spec%disp%show(scaling(iell : min(iell + nscol - 1_IK, size(scaling, 1, IK))), format = scalingFormat, tmsize = 0_IK, bmsize = 0_IK)
565 : end do
566 13 : call spec%disp%skip(count = spec%disp%bmsize)
567 : spec%msg = SKC_"This is the predicted absolute strong-scaling speedup behavior of the "//spec%parallelism%val//&
568 : SKC_" parallelization model, under any MCMC sampling efficiency, given the current parallel communication overhead,&
569 13 : &for increasing numbers of processes, starting from a single process."
570 13 : call spec%disp%note%show(spec%msg)
571 :
572 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
573 :
574 : end block blockParallelSpeedup
575 :
576 :
577 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
578 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end speedup compute %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
579 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
580 :
581 13 : call spec%disp%show("stats.domain.logVolume")
582 13 : call spec%disp%show(spec%domain%logVol, format = format)
583 13 : call spec%disp%note%show("This is the natural logarithm of the volume of the "//spec%ndim%str//SKC_"-dimensional domain over which the objective function was defined.")
584 :
585 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
586 :
587 13 : SET_DRAMDISE(call spec%disp%show("stats.chain.verbose.logFunc.max.val"))
588 : SET_ParaNest(call spec%disp%show("stats.chain.uniques.logFunc.max.val"))
589 13 : call spec%disp%show(stat%chain%mode%val, format = format)
590 13 : call spec%disp%note%show("This is the maximum logFunc value (the maximum of the user-specified objective function).")
591 :
592 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
593 :
594 13 : SET_DRAMDISE(call spec%disp%show("stats.chain.verbose.logFunc.max.crd"))
595 : SET_ParaNest(call spec%disp%show("stats.chain.uniques.logFunc.max.crd"))
596 67 : call spec%disp%show(css_type(spec%domainAxisName%val), format = spec%reportFile%format%fixform, bmsize = 0_IK)
597 13 : call spec%disp%show(stat%chain%mode%crd, format = spec%reportFile%format%allreal)
598 13 : call spec%disp%note%show("This is the coordinates, within the domain of the user-specified objective function, where the maximum function value occurs.")
599 :
600 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
601 :
602 : #if ParaDRAM_ENABLED || ParaDISE_ENABLED
603 13 : call spec%disp%show("stats.chain.compact.logFunc.max.loc")
604 13 : call spec%disp%show(stat%chain%mode%loc, format = format)
605 13 : call spec%disp%note%show("This is the location of the first occurrence of the maximum logFunc in the compact chain.")
606 :
607 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
608 :
609 13 : call spec%disp%show("stats.chain.verbose.logFunc.max.loc")
610 104094 : call spec%disp%show(sum(stat%cfc%sampleWeight(1 : stat%chain%mode%loc - 1)) + 1_IK, format = format) ! stat%chain%mode%Loc%verbose
611 13 : call spec%disp%note%show("This is the location of the first occurrence of the maximum logFunc in the verbose (Markov) chain.")
612 :
613 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
614 :
615 13 : call spec%disp%show("stats.chain.compact.burnin.loc.likelihoodBased")
616 13 : call spec%disp%show(stat%burninLocMCMC%compact, format = format)
617 13 : call spec%disp%note%show("This is the burnin location in the compact chain, based on the occurrence likelihood.")
618 :
619 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
620 :
621 13 : call spec%disp%show("stats.chain.compact.burnin.loc.adaptationBased")
622 13 : call spec%disp%show(stat%burninLocDRAM%compact, format = format)
623 13 : call spec%disp%note%show("This is the burnin location in the compact chain, based on the value of burninAdaptationMeasure simulation specification.")
624 :
625 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
626 :
627 13 : call spec%disp%show("stats.chain.verbose.burnin.loc.likelihoodBased")
628 13 : call spec%disp%show(stat%burninLocMCMC%verbose, format = format)
629 13 : call spec%disp%note%show("This is the burnin location in the verbose (Markov) chain, based on the occurrence likelihood.")
630 :
631 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
632 :
633 13 : call spec%disp%show("stats.chain.verbose.burnin.loc.adaptationBased")
634 13 : call spec%disp%show(stat%burninLocDRAM%verbose, format = format)
635 13 : call spec%disp%note%show("This is the burnin location in the verbose (Markov) chain, based on the value of burninAdaptationMeasure simulation specification.")
636 :
637 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
638 :
639 : ! reset burninLocation to the maximum value
640 :
641 13 : if (stat%burninLocDRAM%compact > stat%burninLocMCMC%compact) then
642 0 : stat%burninLocMCMC%compact = stat%burninLocDRAM%compact
643 0 : stat%burninLocMCMC%verbose = stat%burninLocDRAM%verbose
644 : end if
645 :
646 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
647 : ! Compute the statistical properties of the MCMC chain
648 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
649 :
650 13 : if (spec%image%is%first .and. .not. spec%outputSplashMode%is%silent) then
651 13 : call spec%disp%note%show("Computing the statistical properties of the Markov chain...", unit = output_unit)
652 : end if
653 13 : call spec%disp%text%wrap(NL1//SKC_"The statistical properties of the Markov chain"//NL1)
654 :
655 : ! Compute the covariance and correlation upper-triangle matrices.
656 :
657 : !> \warning
658 : !> forrtl: severe (174): SIGSEGV, segmentation fault occurred
659 : !> Apparently, when the chain is too long (e.g., 300000), the stack size can lead to such an error.
660 : !> The stack size will have to be adjusted in such cases.
661 : !> \todo
662 : !> A robust fix to the above segmentation fault error required further digging into this problem.
663 : !> This could become a serious hard-to-resolve error at higher levels in dynamic languages, especially on Windows.
664 : !> update: The problem is now resolved by allocating on the heap.
665 :
666 13 : call setResized(stat%chain%avg, ndim)
667 39 : call setResized(stat%chain%cov, [ndim, ndim])
668 13 : call setMean(stat%chain%avg, stat%cfc%sampleState(:, stat%burninLocMCMC%compact : stat%cfc%nsam), 2_IK, stat%cfc%sampleWeight(stat%burninLocMCMC%compact : stat%cfc%nsam), stat%chain%size)
669 13 : call setCov(stat%chain%cov, lowDia, stat%chain%avg, stat%cfc%sampleState(:, stat%burninLocMCMC%compact : stat%cfc%nsam), 2_IK, stat%cfc%sampleWeight(stat%burninLocMCMC%compact : stat%cfc%nsam), stat%chain%size)
670 13 : call setMatCopy(stat%chain%cov, rdpack, stat%chain%cov, rdpack, lowDia, transHerm)
671 238 : stat%chain%cor = getCor(stat%chain%cov, subsetv = lowDia)
672 :
673 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
674 :
675 13 : call spec%disp%show("stats.chain.verbose.size.burninExcluded")
676 13 : call spec%disp%show(stat%chain%size, format = format)
677 13 : call spec%disp%note%show("This is the length of the verbose (Markov) Chain excluding burnin.")
678 :
679 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
680 :
681 13 : call spec%disp%show("stats.chain.verbose.avg")
682 67 : call spec%disp%show(css_type(spec%domainAxisName%val), format = spec%reportFile%format%fixform, bmsize = 0_IK)
683 13 : call spec%disp%show(stat%chain%avg, format = spec%reportFile%format%allreal)
684 13 : call spec%disp%note%show("This is the mean (avg) of the verbose (Markov) chain variables.")
685 :
686 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
687 :
688 13 : call spec%disp%show("stats.chain.verbose.std")
689 67 : call spec%disp%show(css_type(spec%domainAxisName%val), format = spec%reportFile%format%fixform, bmsize = 0_IK)
690 67 : call spec%disp%show([(sqrt(stat%chain%cov(idim, idim)), idim = 1, ndim)], format = spec%reportFile%format%allreal)
691 13 : call spec%disp%note%show("This is the standard deviation (std) of the verbose (Markov) chain variables.")
692 :
693 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
694 :
695 13 : call spec%disp%show("stats.chain.verbose.covmat")
696 120 : call spec%disp%show(css_type([character(len(spec%domainAxisName%val),SKC) :: SKC_"axis", spec%domainAxisName%val]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
697 13 : call spec%disp%skip(count = spec%disp%tmsize)
698 40 : do idim = 1, ndim
699 40 : write(spec%reportFile%unit, spec%reportFile%format%strreal) trim(adjustl(spec%domainAxisName%val(idim))), stat%chain%cov(1 : ndim, idim)
700 : end do
701 13 : call spec%disp%skip(count = spec%disp%bmsize)
702 13 : call spec%disp%note%show("This is the covariance matrix of the verbose (Markov) chain.")
703 :
704 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
705 :
706 13 : call spec%disp%show("stats.chain.verbose.cormat")
707 120 : call spec%disp%show(css_type([character(len(spec%domainAxisName%val),SKC) :: SKC_"axis", spec%domainAxisName%val]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
708 13 : call spec%disp%skip(count = spec%disp%tmsize)
709 40 : do idim = 1, ndim
710 40 : write(spec%reportFile%unit, spec%reportFile%format%strreal) trim(adjustl(spec%domainAxisName%val(idim))), stat%chain%cor(1 : ndim, idim)
711 : end do
712 13 : call spec%disp%skip(count = spec%disp%bmsize)
713 13 : call spec%disp%note%show("This is the correlation matrix of the verbose (Markov) chain.")
714 :
715 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
716 :
717 : ! Compute the chain quantiles.
718 : ! \bug Intel ifort bug with heap memory allocations:
719 : ! The Intel ifort cannot digest the following task for an input 2D sample, throwing `double free or corruption (out)`.
720 : ! The source of this error was traced back to the return point from `setExtrap()` within `getQuan()`.
721 : ! The runtime error seems to resolve if the quantile matrix is computed in a loop like the following.
722 :
723 39 : call setResized(stat%chain%quantile%quan, [size(stat%chain%quantile%prob, 1, IK), ndim])
724 : block
725 : real(RKC), allocatable :: sample(:)
726 40 : do idim = 1, ndim
727 1015026 : sample = stat%cfc%sampleState(idim, stat%burninLocMCMC%compact : stat%cfc%nsam)
728 283 : stat%chain%quantile%quan(:, idim) = getQuan(neimean, stat%chain%quantile%prob, sample, stat%cfc%sampleWeight(stat%burninLocMCMC%compact : stat%cfc%nsam), stat%chain%size)
729 : end do
730 : end block
731 : !stat%chain%quantile%quan = getQuan(neimean, stat%chain%quantile%prob, stat%cfc%sampleState(idim, stat%burninLocMCMC%compact : stat%cfc%nsam), 2_IK, stat%cfc%sampleWeight(stat%burninLocMCMC%compact : stat%cfc%nsam), stat%chain%size)
732 :
733 13 : call spec%disp%show("stats.chain.verbose.quantile")
734 120 : call spec%disp%show(css_type([character(len(spec%domainAxisName%val),SKC) :: SKC_"quantile", spec%domainAxisName%val]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
735 13 : call spec%disp%skip(count = spec%disp%tmsize)
736 130 : do iq = 1, size(stat%chain%quantile%prob, 1, IK)
737 373 : write(spec%reportFile%unit, spec%reportFile%format%allreal) stat%chain%quantile%prob(iq), (stat%chain%quantile%quan(iq, idim), idim = 1, ndim)
738 : end do
739 13 : call spec%disp%skip(count = spec%disp%bmsize)
740 13 : call spec%disp%note%show("This is the quantiles table of the variables of the verbose (Markov) chain.")
741 : #endif
742 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
743 : ! Generate the i.i.d. sample statistics and output file (if requested)
744 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
745 :
746 : ! report refined sample statistics, and generate output refined sample if requested.
747 :
748 13 : if (spec%image%is%first .and. .not. spec%outputSplashMode%is%silent) then
749 13 : call spec%disp%note%show("Computing the final refined i.i.d. sample size...", unit = output_unit)
750 : end if
751 : #if ParaDRAM_ENABLED || ParaDISE_ENABLED
752 : err = stat%sfc%getErrRefinement ( sampleWeight = stat%cfc%sampleWeight(stat%burninLocMCMC%compact :) & ! LCOV_EXCL_LINE
753 : , sampleLogFunc = stat%cfc%sampleLogFunc(stat%burninLocMCMC%compact :) & ! LCOV_EXCL_LINE
754 : , sampleState = stat%cfc%sampleState(:, stat%burninLocMCMC%compact :) & ! LCOV_EXCL_LINE
755 : , outputSampleRefinementCount = spec%outputSampleRefinementCount%val & ! LCOV_EXCL_LINE
756 : , outputSampleRefinementMethod = spec%outputSampleRefinementMethod%val & ! LCOV_EXCL_LINE
757 13 : )
758 :
759 : ! compute the maximum integrated autocorrelation times for each variable.
760 :
761 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
762 :
763 13 : call spec%disp%show("stats.chain.refined.act")
764 13 : call spec%disp%skip(count = spec%disp%tmsize)
765 40 : write(spec%reportFile%unit, spec%reportFile%format%fixform) "refinementStage", "sampleSize", "ACT_sampleLogFunc", (SKC_"ACT_"//trim(adjustl(spec%domainAxisName%val(idim))), idim = 1, ndim)
766 13 : block
767 : integer(IK) :: iref
768 : character(:,SKC), allocatable :: thisForm
769 : thisForm = SKC_"('"//spec%reportFile%indent//SKC_"',2(I"//spec%outputColumnWidth%max//SKC_",' '),*("//& ! LCOV_EXCL_LINE
770 13 : spec%real%ed//spec%outputColumnWidth%max//SKC_"."//spec%outputPrecision%str//spec%real%ex//SKC_",:,' '))"
771 13 : call spec%disp%skip(count = spec%disp%tmsize)
772 52 : do iref = 0, stat%sfc%nref
773 52 : write(spec%reportFile%unit, thisForm) iref, stat%sfc%sumw(iref), stat%sfc%act(0 : ndim, iref)
774 : end do
775 26 : call spec%disp%skip(count = spec%disp%bmsize)
776 : end block
777 13 : spec%msg = SKC_"This is the table of the Integrated Autocorrelation (ACT) of individual variables in the verbose (Markov) chain, at increasing stages of chain refinements."
778 : if (stat%sfc%nref == 0_IK) spec%msg = spec%msg//NL1// & ! LCOV_EXCL_LINE
779 : SKC_"The user-specified `outputSampleRefinementCount` ("//getStr(spec%outputSampleRefinementCount%val)//SKC_") &
780 0 : &is too small to ensure an accurate computation of the decorrelated i.i.d. effective sample size. No refinement of the Markov chain was performed."
781 13 : call spec%disp%note%show(spec%msg)
782 :
783 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
784 :
785 : ! Report the final Effective Sample Size (ESS) based on ACT.
786 :
787 13 : stat%ess = stat%sfc%sumw(stat%sfc%nref)
788 :
789 : #elif ParaNest_ENABLED
790 :
791 : ! Read the ParaNest output chain file contents.
792 :
793 : err = getErrChainRead(stat%cfc, file = spec%chainFile%file, spec%outputChainFileFormat%val, pre = stat%numFunCallAccepted)
794 : if (err%occurred) then
795 : err%msg = PROCEDURE_NAME//err%msg ! LCOV_EXCL_LINE
796 : exit blockLeaderPostProcessing ! LCOV_EXCL_LINE
797 : return ! LCOV_EXCL_LINE
798 : end if
799 :
800 : ! Compute the effective sample size.
801 :
802 : stat%ess = nint(sum(exp(stat%cfc%sampleLogWeight - maxval(stat%cfc%sampleLogWeight))), IK)
803 : #else
804 : #error "Unrecognized interface."
805 : #endif
806 13 : stat%sample%size = abs(spec%outputSampleSize%val)
807 13 : if (spec%outputSampleSize%val < 0_IK) stat%sample%size = stat%sample%size * stat%ess
808 :
809 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
810 :
811 13 : call spec%disp%show("stats.chain.refined.ess")
812 13 : call spec%disp%show(stat%ess, format = format)
813 13 : call spec%disp%note%show("This is the estimated Effective (i.i.d.) Sample Size (ESS) of the final refined chain.")
814 :
815 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
816 :
817 13 : call spec%disp%show("stats.chain.refined.efficiency.essOverAccepted")
818 13 : call spec%disp%show(real(stat%ess, RKC) / real(stat%numFunCallAccepted, RKC), format = format)
819 : spec%msg = SKC_"This is the effective sampling efficiency given the accepted function calls, that is, &
820 13 : &the final refined effective sample size (ESS) divided by the number of accepted function calls."
821 13 : call spec%disp%note%show(spec%msg)
822 :
823 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
824 :
825 13 : call spec%disp%show("stats.chain.refined.efficiency.essOverAcceptedRejected")
826 13 : call spec%disp%show(real(stat%ess, RKC) / real(stat%numFunCallAcceptedRejected, RKC), format = format)
827 : spec%msg = SKC_"This is the effective sampling efficiency given the accepted and rejected function calls, that is, &
828 13 : &the final refined effective sample size (ESS) divided by the number of (accepted + rejected) function calls."
829 13 : call spec%disp%note%show(spec%msg)
830 :
831 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
832 :
833 13 : call spec%disp%show("stats.chain.refined.efficiency.essOverAcceptedRejectedDelayed")
834 13 : call spec%disp%show(real(stat%ess, RKC) / real(stat%numFunCallAcceptedRejectedDelayed, RKC), format = format)
835 : spec%msg = SKC_"This is the effective sampling efficiency given the accepted, rejected, and delayed-rejection (if any requested) function calls, &
836 13 : &that is, the final refined effective sample size (ESS) divided by the number of (accepted + rejected + delayed-rejection) function calls."
837 13 : call spec%disp%note%show(spec%msg)
838 :
839 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
840 :
841 13 : call spec%disp%show("stats.chain.refined.efficiency.essOverAcceptedRejectedDelayedUnused")
842 13 : call spec%disp%show(real(stat%ess, RKC) / real(stat%numFunCallAcceptedRejectedDelayedUnused, RKC), format = format)
843 : spec%msg = SKC_"This is the effective sampling efficiency given the accepted, rejected, delayed-rejection (if any requested), and unused function calls, &
844 13 : &(in parallel simulations), that is, the final refined effective sample size (ESS) divided by the number of (accepted + rejected + delayed-rejection + unused) function calls."
845 13 : call spec%disp%note%show(spec%msg)
846 :
847 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
848 :
849 : !end associate blockEffectiveSampleSize
850 :
851 : ! Generate the output refined sample if requested.
852 :
853 13 : blockSampleFileGeneration: if (spec%outputSampleSize%val == 0_IK) then
854 :
855 0 : call spec%disp%note%show(SKC_"Skipping the generation of the refined chain and the output "//spec%sampleFile%kind//SKC_" file, as requested by the user...")
856 :
857 : else blockSampleFileGeneration
858 :
859 : ! report to the report-file(s)
860 :
861 13 : call spec%disp%note%show(SKC_"Generating the output "//spec%sampleFile%kind//SKC_" file:"//NL1//spec%sampleFile%file)
862 13 : if (spec%image%is%first .and. .not. spec%outputSplashMode%is%silent) then
863 : ! print the message for the generating the output sample file on the first image
864 13 : call spec%disp%note%show(SKC_"Generating the output "//spec%sampleFile%kind//SKC_" file:"//NL1//spec%sampleFile%file, unit = output_unit, bmsize = 0_IK)
865 : ! print the message for the generating the output sample file on the rest of the images in order.
866 : #if CAF_ENABLED || MPI_ENABLED || OMP_ENABLED
867 : if (spec%parallelism%is%multiChain) then
868 : block
869 : integer(IK) :: imageID
870 : do imageID = 2, spec%image%count
871 : call spec%disp%note%show(getReplaced(spec%sampleFile%file, SKC_"pid1", SKC_"pid"//getStr(imageID)), unit = output_unit, tmsize = 0_IK, bmsize = 0_IK)
872 : end do
873 : end block
874 : end if
875 : #endif
876 13 : call spec%disp%skip(unit = output_unit)
877 : end if
878 :
879 : ! Begin sample file generation.
880 :
881 133 : stat%sfc%colname = [css_type(chainFileColName(size(chainFileColName, 1, IK))), css_type(spec%domainAxisName%val)]
882 13 : stat%sfc%header = getStr(stat%sfc%colname, format = spec%sampleFile%format%header)
883 : #if ParaDISE_ENABLED || ParaDRAM_ENABLED
884 13 : if (spec%outputSampleSize%val /= -1_IK) then
885 : ! Regenerate the refined sample, this time with the user-specified sample size.
886 : block
887 : integer(IK), allocatable :: cumSumWeight(:), unifrnd(:)
888 : integer(IK) :: isam, iloc
889 7 : stat%sfc%nref = 0_IK
890 7 : call setRebound(stat%sfc%nsam, 0_IK, 0_IK)
891 7 : call setRebound(stat%sfc%sumw, 0_IK, 0_IK)
892 7 : stat%sfc%nsam(stat%sfc%nref) = stat%numFunCallAccepted - stat%burninLocMCMC%compact + 1_IK
893 7 : call setResized(cumSumWeight, stat%sfc%nsam(stat%sfc%nref))
894 7 : call setCumSum(cumSumWeight, stat%cfc%sampleWeight(stat%burninLocMCMC%compact :))
895 14 : unifrnd = getUnifRand(spec%rng, 1_IK, cumSumWeight(stat%sfc%nsam(stat%sfc%nref)), stat%sample%size)
896 21 : call setRebound(stat%sfc%sampleLogFuncState, [0_IK, 1_IK], [ndim, stat%sfc%nsam(stat%sfc%nref)])
897 70007 : stat%sfc%sampleLogFuncState(1 :, :) = stat%cfc%sampleState(:, stat%burninLocMCMC%compact :)
898 35007 : stat%sfc%sampleLogFuncState(0, :) = stat%cfc%sampleLogFunc(stat%burninLocMCMC%compact :)
899 35014 : stat%sfc%sampleWeight = getFilled(0_IK, stat%sfc%nsam(stat%sfc%nref))
900 17507 : loopOverSample: do isam = 1, stat%sample%size
901 43487256 : do iloc = 1, size(cumSumWeight, 1, IK)
902 43487256 : if (cumSumWeight(iloc) < unifrnd(isam)) cycle
903 17500 : stat%sfc%sampleWeight(iloc) = stat%sfc%sampleWeight(iloc) + 1_IK
904 43469756 : cycle loopOverSample
905 : end do
906 7 : error stop "This is a miracle! Please report it to the developers at: https://github.com/cdslaborg/paramonte"
907 : end do loopOverSample
908 35007 : stat%sfc%sumw(stat%sfc%nref) = sum(stat%sfc%sampleWeight)
909 : end block
910 : end if
911 : #elif ParaNest_ENABLED
912 : stat%sfc%nref = 0_IK
913 : call setRebound(stat%sfc%nsam, 0_IK, 0_IK)
914 : call setRebound(stat%sfc%sumw, 0_IK, 0_IK)
915 : stat%sfc%sampleWeight = nint(getReweight(exp(stat%cfc%sampleLogWeight - getLogSumExp(stat%cfc%sampleLogWeight)), 1._RKC / stat%sample%size), IK)
916 : stat%sfc%nsam(stat%sfc%nref) = size(stat%sfc%sampleWeight, 1, IK)
917 : call setRebound(stat%sfc%sampleLogFuncState, [0_IK, 1_IK], [ndim, stat%sfc%nsam(stat%sfc%nref)])
918 : stat%sfc%sumw(stat%sfc%nref) = sum(stat%sfc%sampleWeight)
919 : #else
920 : #error "Unrecognized interface."
921 : #endif
922 : ! open the output sample file and write the sample.
923 :
924 : sfc_block: block
925 : integer(IK) :: isam, iwei
926 13 : spec%sampleFile%iomsg = SKC_""
927 13 : spec%sampleFile%unit = getFileUnit() ! for some unknown reason, if newunit is used, GFortran opens the file as an internal file
928 13 : open(unit = spec%sampleFile%unit, file = spec%sampleFile%file, form = spec%sampleFile%form, status = spec%sampleFile%status, position = spec%sampleFile%position, iostat = spec%sampleFile%iostat, iomsg = spec%sampleFile%iomsg SHARED)
929 13 : if (spec%sampleFile%iostat /= 0_IK) exit sfc_block
930 13 : write(spec%sampleFile%unit, spec%sampleFile%format%header, iostat = spec%sampleFile%iostat, iomsg = spec%sampleFile%iomsg) stat%sfc%header
931 13 : if (spec%sampleFile%iostat /= 0_IK) exit sfc_block
932 104884 : do isam = 1, stat%sfc%nsam(stat%sfc%nref)
933 193879 : do iwei = 1, stat%sfc%sampleWeight(isam)
934 88995 : write(spec%sampleFile%unit, spec%sampleFile%format%rows, iostat = spec%sampleFile%iostat, iomsg = spec%sampleFile%iomsg) stat%sfc%sampleLogFuncState(:, isam)
935 193866 : if (spec%sampleFile%iostat /= 0_IK) exit sfc_block
936 : end do
937 : end do
938 13 : close(spec%sampleFile%unit, iostat = spec%sampleFile%iostat, iomsg = spec%sampleFile%iomsg)
939 : if (spec%sampleFile%iostat /= 0_IK) exit sfc_block
940 : end block sfc_block
941 13 : err%occurred = spec%sampleFile%iostat /= 0_IK
942 13 : if (err%occurred) then
943 0 : err%stat = spec%sampleFile%iostat
944 0 : err%msg = spec%sampleFile%iomsg
945 0 : exit blockLeaderPostProcessing
946 : end if
947 :
948 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
949 : ! Compute the statistical properties of the refined sample
950 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
951 :
952 13 : if (spec%image%is%first .and. .not. spec%outputSplashMode%is%silent) then
953 13 : call spec%disp%note%show(SKC_"Computing the statistical properties of the final output sample...", unit = output_unit)
954 : end if
955 13 : call spec%disp%text%wrap(NL1//SKC_"The statistical properties of the final output sample"//NL1)
956 : #if ParaDRAM_ENABLED || ParaDISE_ENABLED
957 52 : CHECK_ASSERTION(__LINE__, stat%sample%size == stat%sfc%sumw(stat%sfc%nref), \
958 : SK_"@getErrSampling(): The condition `stat%sample%size == stat%sfc%sumw(stat%sfc%nref)` must hold. stat%sample%size, stat%sfc%sumw(stat%sfc%nref), stat%sfc%nref = "//\
959 : getStr([stat%sample%size, stat%sfc%sumw(stat%sfc%nref), stat%sfc%nref]))
960 : #endif
961 : ! Compute the covariance and correlation upper-triangle matrices.
962 :
963 13 : call setRebound(stat%sample%avg, 0_IK, ndim)
964 39 : call setRebound(stat%sample%cov, [0_IK, 0_IK], [ndim, ndim])
965 13 : call setMean(stat%sample%avg, stat%sfc%sampleLogFuncState(:, 1 : stat%sfc%nsam(stat%sfc%nref)), 2_IK, stat%sfc%sampleWeight(1 : stat%sfc%nsam(stat%sfc%nref)), stat%sfc%sumw(stat%sfc%nref))
966 13 : if (ndim < stat%sfc%nsam(stat%sfc%nref)) then
967 12 : call setCov(stat%sample%cov, lowDia, stat%sample%avg, stat%sfc%sampleLogFuncState(:, 1 : stat%sfc%nsam(stat%sfc%nref)), 2_IK, stat%sfc%sampleWeight(1:stat%sfc%nsam(stat%sfc%nref)), stat%sfc%sumw(stat%sfc%nref))
968 12 : call setMatCopy(stat%sample%cov, rdpack, stat%sample%cov, rdpack, lowDia, transHerm)
969 372 : stat%sample%cor = getCor(stat%sample%cov, subsetv = lowDia)
970 : else
971 16 : stat%sample%cor = getMatInit([ndim, ndim], uppLowDia, 0._RKC, 0._RKC, 0._RKC)
972 16 : stat%sample%cov = getMatInit([ndim, ndim], uppLowDia, 0._RKC, 0._RKC, 0._RKC)
973 : end if
974 :
975 : ! Report the refined chain statistics.
976 :
977 : !formatStr = "(1A"//spec%outputColumnWidth%max//SKC_",*(E"//spec%outputColumnWidth%max//SKC_"."//spec%outputPrecision%str//SKC_"))"
978 :
979 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
980 :
981 13 : call spec%disp%show("stats.chain.refined.length")
982 13 : call spec%disp%show(stat%sample%size, format = format)
983 13 : spec%msg = "This is the final output refined sample size."
984 13 : if (spec%outputSampleSize%val /= -1_IK) then
985 7 : if (stat%sample%size < stat%ess) then
986 : spec%msg = spec%msg//NL1//SKC_"The user-requested sample size ("//getStr(stat%sample%size)//SKC_") is smaller &
987 : &than the potentially-optimal i.i.d. sample size ("//getStr(stat%ess)//SKC_"). The output sample &
988 : &contains i.i.d. samples, however, the sample-size could have been larger if it had been set to the optimal size. &
989 0 : &To get the optimal size in the future runs, set outputSampleSize = -1, or drop it from the input list."
990 7 : elseif (stat%ess < stat%sample%size) then
991 : spec%msg = spec%msg//NL1//SKC_"The user-requested sample size ("//getStr(stat%sample%size)//SKC_") is larger than &
992 : &the potentially-optimal i.i.d. sample size ("//getStr(stat%ess)//SKC_"). The resulting sample likely contains &
993 : &duplicates and is not independently and identically distributed (i.i.d.). To get the optimal &
994 7 : &size in the future runs, set outputSampleSize = -1, or drop it from the input list."
995 : else ! LCOV_EXCL_LINE
996 : spec%msg = spec%msg//NL1//SKC_"How lucky that could be! The user-requested sample size ("//getStr(stat%sample%size)//& ! LCOV_EXCL_LINE
997 : SKC_") is equal to the potentially-optimal i.i.d. sample size determined by the "//spec%method%val//SKC_" sampler." ! LCOV_EXCL_LINE
998 : end if
999 : end if
1000 13 : call spec%disp%note%show(spec%msg)
1001 :
1002 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1003 :
1004 13 : call spec%disp%show("stats.chain.refined.avg")
1005 67 : call spec%disp%show(css_type(spec%domainAxisName%val), format = spec%reportFile%format%fixform, bmsize = 0_IK)
1006 13 : call spec%disp%show(stat%sample%avg(1:), format = spec%reportFile%format%allreal)
1007 13 : call spec%disp%note%show("This is the mean (avg) of the final output refined sample.")
1008 :
1009 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1010 :
1011 13 : call spec%disp%show("stats.chain.refined.std")
1012 67 : call spec%disp%show(css_type(spec%domainAxisName%val), format = spec%reportFile%format%fixform, bmsize = 0_IK)
1013 67 : call spec%disp%show([(sqrt(stat%sample%cov(idim, idim)), idim = 1, ndim)], format = spec%reportFile%format%allreal)
1014 13 : call spec%disp%note%show("This is the standard deviation (std) of the final output refined sample.")
1015 :
1016 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1017 :
1018 13 : call spec%disp%show("stats.chain.refined.covmat")
1019 120 : call spec%disp%show(css_type([character(len(spec%domainAxisName%val),SKC) :: SKC_"axis", spec%domainAxisName%val]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
1020 13 : call spec%disp%skip(count = spec%disp%tmsize)
1021 40 : do idim = 1, ndim
1022 40 : write(spec%reportFile%unit, spec%reportFile%format%strreal) trim(adjustl(spec%domainAxisName%val(idim))), stat%sample%cov(1 : ndim, idim)
1023 : end do
1024 13 : call spec%disp%skip(count = spec%disp%bmsize)
1025 13 : call spec%disp%note%show("This is the covariance matrix of the final output refined sample.")
1026 13 : if (.not. isMatClass(stat%sample%cov, posdefmat)) call spec%disp%warn%show("The sample covariance matrix is not positive-definite. sample size = "//getStr(stat%sfc%nsam(stat%sfc%nref)))
1027 :
1028 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1029 :
1030 13 : call spec%disp%show("stats.chain.refined.cormat")
1031 120 : call spec%disp%show(css_type([character(len(spec%domainAxisName%val),SKC) :: SKC_"axis", spec%domainAxisName%val]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
1032 13 : call spec%disp%skip(count = spec%disp%tmsize)
1033 40 : do idim = 1, ndim
1034 40 : write(spec%reportFile%unit, spec%reportFile%format%strreal) trim(adjustl(spec%domainAxisName%val(idim))), stat%sample%cor(1 : ndim, idim)
1035 : end do
1036 13 : call spec%disp%skip(count = spec%disp%bmsize)
1037 13 : call spec%disp%note%show("This is the correlation matrix of the final output refined sample.")
1038 13 : if (.not. isMatClass(stat%sample%cor, posdefmat)) call spec%disp%warn%show("The sample correlation matrix is not positive-definite. sample size = "//getStr(stat%sfc%nsam(stat%sfc%nref)))
1039 :
1040 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1041 :
1042 : ! Compute the sample quantiles.
1043 : ! \bug Intel ifort bug with heap memory allocations:
1044 : ! The Intel ifort cannot digest the following task for an input 2D sample, throwing `double free or corruption (out)`.
1045 : ! The source of this error was traced back to the return point from `setExtrap()` within `getQuan()`.
1046 : ! The runtime error seems to resolve if the quantile matrix is computed in a loop like the following.
1047 :
1048 39 : call setResized(stat%sample%quantile%quan, [size(stat%sample%quantile%prob, 1, IK), ndim])
1049 : block
1050 : real(RKC), allocatable :: sample(:)
1051 40 : do idim = 1, ndim
1052 314252 : sample = stat%sfc%sampleLogFuncState(idim, :)
1053 283 : stat%sample%quantile%quan(:, idim) = getQuan(neimean, stat%sample%quantile%prob, sample, stat%sfc%sampleWeight, stat%sfc%sumw(stat%sfc%nref))
1054 : end do
1055 : end block
1056 : !stat%sample%quantile%quan = getQuan(neimean, stat%sample%quantile%prob, stat%sfc%sampleLogFuncState, 2_IK, stat%sfc%sampleWeight, stat%sfc%sumw(stat%sfc%nref))
1057 :
1058 13 : call spec%disp%show("stats.chain.refined.quantile")
1059 120 : call spec%disp%show(css_type([character(len(spec%domainAxisName%val),SKC) :: "quantile", spec%domainAxisName%val]), format = spec%reportFile%format%fixform, bmsize = 0_IK)
1060 13 : call spec%disp%skip(count = spec%disp%tmsize)
1061 130 : do iq = 1, size(stat%sample%quantile%prob, 1, IK)
1062 373 : write(spec%reportFile%unit, spec%reportFile%format%allreal) stat%sample%quantile%prob(iq), (stat%sample%quantile%quan(iq, idim), idim = 1, ndim)
1063 : end do
1064 13 : call spec%disp%skip(count = spec%disp%bmsize)
1065 13 : call spec%disp%note%show("This is the quantiles table of the variables of the final output refined sample.")
1066 :
1067 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1068 : ! Begin inter-chain convergence test in multiChain parallelization mode
1069 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1070 :
1071 : #if CAF_ENABLED || MPI_ENABLED || OMP_ENABLED
1072 : blockInterChainConvergence: if (spec%parallelism%is%multiChain .and. spec%image%count > 1_IK) then
1073 :
1074 : call spec%disp%note%show("Computing the inter-chain convergence probabilities...")
1075 : if (spec%image%is%first .and. .not. spec%outputSplashMode%is%silent) then
1076 : call spec%disp%note%show("Computing the inter-chain convergence probabilities...", unit = output_unit, bmsize = 0_IK)
1077 : end if
1078 : call spec%image%sync()
1079 :
1080 : ! compute and report the KS convergence probabilities for all images.
1081 :
1082 : multiChainConvergenceTest: block
1083 :
1084 : real(RKC), allocatable :: sampleLogFuncState1(:,:)
1085 : real(RKC), allocatable :: sampleLogFuncState2(:,:)
1086 : integer(IK) :: imageID, idimMinProbKS, pidMinProbKS
1087 : character(:, SK), allocatable :: sampleFilePath
1088 : real(RKC), allocatable :: probKS(:)
1089 : real(RKC) :: minProbKS
1090 :
1091 : minProbKS = 2._RKC !huge(minProbKS)
1092 : call setResized(probKS, ndim + 1_IK)
1093 : !call setRebound(probKS, 0_IK, ndim)
1094 :
1095 : ! sort the refined chain on the current image.
1096 :
1097 : err%stat = getErrTableRead(spec%sampleFile%file, sampleLogFuncState1, trans, sep = spec%outputSeparator%val, roff = 1_IK)
1098 : err%occurred = err%stat /= 0_IK
1099 : if (err%occurred) then
1100 : err%msg = PROCEDURE_NAME//SKC_": "//err%msg ! LCOV_EXCL_LINE
1101 : exit blockLeaderPostProcessing ! LCOV_EXCL_LINE
1102 : end if
1103 : do idim = 1, ndim + 1
1104 : call setSorted(sampleLogFuncState1(:, idim))
1105 : end do
1106 :
1107 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1108 :
1109 : call spec%disp%show("stats.chain.refined.kstest.prob")
1110 : call spec%disp%show([css_type("processID"), stat%sfc%colname], format = spec%reportFile%format%fixform, bmsize = 0_IK)
1111 :
1112 : do imageID = 1, spec%image%count
1113 :
1114 : if (spec%image%id /= imageID) then
1115 :
1116 : ! read the refined chain on the other image.
1117 : sampleFilePath = getReplaced(spec%sampleFile%file, spec%sampleFile%suffix, getReplaced(spec%sampleFile%suffix, SKC_"_pid"//getStr(spec%image%id), SKC_"_pid"//getStr(imageID)))
1118 : err%stat = getErrTableRead(sampleFilePath, sampleLogFuncState2, trans, sep = spec%outputSeparator%val, roff = 1_IK)
1119 : err%occurred = err%stat /= 0_IK
1120 : if (err%occurred) then
1121 : err%msg = PROCEDURE_NAME//SKC_": "//err%msg ! LCOV_EXCL_LINE
1122 : exit blockLeaderPostProcessing ! LCOV_EXCL_LINE
1123 : end if
1124 :
1125 : do idim = 1, ndim + 1
1126 : ! sort the refined chain on the other image.
1127 : call setSorted(sampleLogFuncState2(:, idim))
1128 : ! compute the inter-chain KS probability table.
1129 : probKS(idim) = getProbKS(getDisKolm(sampleLogFuncState1(:, idim), sampleLogFuncState2(:, idim), ascending), size(sampleLogFuncState1, 1, IK), size(sampleLogFuncState2, 1, IK))
1130 : if (minProbKS <= probKS(idim)) cycle
1131 : minProbKS = probKS(idim)
1132 : pidMinProbKS = imageID
1133 : idimMinProbKS = idim
1134 : end do
1135 :
1136 : ! write the inter-chain KS probability table row
1137 :
1138 : write(spec%reportFile%unit, spec%reportFile%format%intreal) imageID, probKS
1139 :
1140 : end if
1141 :
1142 : end do
1143 :
1144 : call spec%disp%skip(count = spec%disp%bmsize)
1145 : spec%msg = SKC_"This is the table of pairwise inter-chain Kolmogorov-Smirnov (KS) convergence (similarity) probabilities. &
1146 : &Higher KS probabilities are better, presenting less evidence for a lack of convergence of all chains to the same target density function."
1147 : call spec%disp%note%show(spec%msg)
1148 :
1149 : ! write the smallest KS probabilities
1150 :
1151 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1152 :
1153 : call spec%disp%show("stats.chain.refined.kstest.prob.min")
1154 : call spec%disp%show(minProbKS, format = format)
1155 : spec%msg = SKC_"This is the smallest KS-test probability for the inter-chain sampling convergence, which has happened between "//&
1156 : stat%sfc%colname(idimMinProbKS)%val//SKC_" on the chains generated by processes "//getStr(spec%image%id)//SKC_" and "//getStr(pidMinProbKS)//SKC_"."
1157 : call spec%disp%note%show(spec%msg)
1158 :
1159 : ! Report the smallest KS probabilities on stdout.
1160 :
1161 : if (.not. spec%outputSplashMode%is%silent) then
1162 : if (spec%image%is%first) then
1163 : call spec%disp%note%show("The smallest KS probabilities for the inter-chain sampling convergence (higher is better):", unit = output_unit, bmsize = 0_IK)!, tmsize = 2_IK
1164 : end if
1165 : do imageID = 1, spec%image%count
1166 : if (spec%image%id == imageID) then
1167 : spec%msg = getStr(minProbKS)//SKC_" for "//stat%sfc%colname(idimMinProbKS)%val//&
1168 : SKC_" on the chains generated by processes "//getStr(spec%image%id)//SKC_" and "//getStr(pidMinProbKS)//SKC_"."
1169 : call spec%disp%note%show(spec%msg, unit = output_unit, tmsize = 0_IK, bmsize = 0_IK)
1170 : end if
1171 : flush(output_unit)
1172 : call spec%image%sync()
1173 : end do
1174 : if (spec%image%is%first) then
1175 : call execute_command_line("", cmdstat = err%stat)
1176 : if (1 < spec%disp%bmsize) call spec%disp%skip(unit = output_unit, count = spec%disp%bmsize - 1)
1177 : end if
1178 : end if
1179 :
1180 : end block multiChainConvergenceTest
1181 :
1182 : call spec%image%sync()
1183 :
1184 : end if blockInterChainConvergence
1185 : #endif
1186 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1187 : ! End inter-chain convergence test in multiChain parallelization mode
1188 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1189 :
1190 : end if blockSampleFileGeneration
1191 :
1192 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1193 : ! End of generating the i.i.d. sample statistics and output file (if requested)
1194 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1195 :
1196 : ! Mission accomplished.
1197 :
1198 13 : call spec%disp%note%show("Mission Accomplished.", tmsize = 2_IK * spec%disp%note%tmsize)!, tmsize = 3_IK, bmsize = 1_IK
1199 13 : if (spec%reportFile%unit /= output_unit .and. spec%image%is%first .and. .not. spec%outputSplashMode%is%silent) then
1200 13 : flush(output_unit)
1201 13 : call execute_command_line("", cmdstat = err%stat)
1202 13 : call spec%disp%note%show("Mission Accomplished.", unit = output_unit, tmsize = 1_IK)!, tmsize = 1_IK, bmsize = 3_IK)
1203 : end if
1204 13 : close(spec%progressFile%unit)
1205 13 : close(spec%reportFile%unit)
1206 39 : err%msg = SK_""
1207 :
1208 : end associate
1209 :
1210 : end if blockLeaderPostProcessing
1211 :
1212 : ! A global sync is essential for parallel applications.
1213 :
1214 : SET_CAFMPI(call spec%image%sync())
1215 : SET_CAFMPI(if (spec%parallelismMpiFinalizeEnabled%val) call spec%image%finalize())
1216 : #else
1217 : !%%%%%%%%%%%%%%%%%%%%%%%%
1218 : #error "Unrecognized interface."
1219 : !%%%%%%%%%%%%%%%%%%%%%%%%
1220 : #endif
1221 : #undef RETURN_IF_FAILED
1222 : #undef SET_DRAMDISE
1223 : #undef SET_ParaNest
1224 : #undef SET_CAFMPI
1225 : #undef SET_SERIAL
1226 : #undef SET_DEBUG
1227 : #undef SET_CAF
1228 : #undef SET_MPI
1229 : #undef SET_OMP
1230 : #undef COINDEX
1231 : #undef TRACE
|
