Program snippet that shows how to define the integrand subroutine (named fxn here):
EXTERNAL fxn
[ main program ]
SUBROUTINE fxn(x, k, wgt, aux, f)
IMPLICIT NONE
DOUBLE PRECISION x(5), wgt, aux(1), f(1)
INTEGER k(1)
[...]
END
The array x contains the random numbers for the continuous dimensions,
the array k contains the random numbers for the discrete dimensions,
the array aux contains the auxilliary random numbers. The VEGAS weight (For
each iteration the sum of the weights is normalized to 1.)
is passed in variable wgt, so that it can be passed on to routines that
fill histograms, for example. The subroutine fxn is expected to fill the
array f with the function values that correspond to the arguments
x(1), x(2), ..., k(1), k(2), ... and aux(1), aux(2), ...
The value of the first function is assigned to f(1). The integration is optimized
for this function according to the VEGAS algorithm. Secondary functions can simultaneously be
integrated and their values are expected to be returned as f(2), f(3), etc.
First this subroutine has to be called:
CALL dvegas_init(cdim, nsdim, cnbin, ddim, adim, nvals, fn, fxn)
which sets the following (internal) INTEGER variables:
cdim: number of continuous dimensions
nsdim: number of non-separable continuous dimensions (first nsdim dimensions)
cnbin: number of bins for each continuous dimension
ddim: number of discrete dimensions
adim: number of auxilliary dimensions (no adaptation)
nvals: array: number of values for each discrete dimension (see below)
fn: number of functions to integrate, (weights are adapted for first function)
fxn: integrand subroutine (this is a subroutine identifier, not an INTEGER variable)
The default integration/summation range is
Random numbers in (0, 1) are provided to facilitate the addition of auxilliary dimensions that are not VEGAS driven.
The parameters ddim or adim can be set to zero if no discrete
or auxilliary dimensions are needed.
Optional initialization calls:
CALL dvegas_init_parallel(ncpu)
specify number of CPUs for multiprocessor machines (ncpu <= 4, see section 5.1)
CALL dvegas_init_output()
save status (init. parameters and weights) to file dvegas.dat-new after next dvegas call
CALL dvegas_init_files(outputfile, inputfile)
restore previously saved status from file inputfile and save new status to file outputfile
After dVegas has been initialized you can make multiple calls to dvegas:
CALL dvegas(ncall, iterats, init, integral, error, redchi2)
ncall: INTEGER*8 (!), request ncall shots per iteration INTEGER, request iterats iterations, init is explained below DOUBLE PRECISION arrays, pass accumulated results, errors and "reduced chi^2" = chi^2/(nr of estimates-1)
The parameter init has the same meaning as in the original VEGAS program:
init = 0: fresh start (no prior information)
init = 1: use previous adaptation info, but discard shot data
init = 2: use previous adaptation info and build on previous shots data
PROGRAM TESTDVEGAS
IMPLICIT NONE
INTEGER nvals(1)
INTEGER*8 ncalls
DOUBLE PRECISION integral(1), error(1), redchi2(1)
CHARACTER*20 outputfile, inputfile
EXTERNAL fxn2
nvals(1) = 100
ncalls = 100000
CALL dvegas_init(5, 0, 50, 1, 0, nvals, 1, fxn2)
* 2 iterations with init = 0, no initial adaptation
CALL dvegas(ncalls, 2, 0, integral, error, redchi2)
* result is 199.0
WRITE(*,*) 'Result:'
WRITE(*,*) 'integral = ', integral(1)
WRITE(*,*) 'error = ', error(1)
WRITE(*,*) 'chi2/itn = ', redchi2(1)
WRITE(*,*) ' '
* now 1 iteration with init = 1, use previous adaptation
* save grid
CALL dvegas_init_output()
CALL dvegas(ncalls, 1, 1, integral, error, redchi2)
* now another run after restoring adaptation from file
outputfile = 'run2.out'
inputfile = 'run2.in'
CALL dvegas_init_files(outputfile, inputfile)
CALL dvegas(ncalls, 1, 1, integral, error, redchi2)
END
SUBROUTINE fxn2(x, k, wgt, aux, f)
IMPLICIT NONE
DOUBLE PRECISION x(5), wgt, aux(1), f(1)
INTEGER k(1)
IF (k(1) .LT. 99) THEN
f(1) = 1 * x(1)*x(2)*x(3)*x(4)*x(5)* 2.0**5
ELSE IF (k(1) .EQ. 99) THEN
f(1) = 100 * x(1)*x(2)*x(3)*x(4)*x(5)* 2.0**5
END IF
END
CALL dvegas_init(cdim, nsdim, cnbin, ddim, adim, nvals, fn)
which sets the following (internal) INTEGER variables:
cdim: number of continuous dimensions
nsdim: number of non-separable continuous dimensions (first nsdim dimensions)
cnbin: number of bins for each continuous dimension
ddim: number of discrete dimensions
adim: number of auxilliary dimensions (no adaptation)
nvals: array: number of values for each discrete dimension (see below)
fn: number of functions to integrate, (weights are adapted for first function)
The default integration/summation range is
Random numbers in (0, 1) are provided to facilitate the addition of auxilliary dimensions that are not VEGAS driven.
The parameters ddim or adim can be set to zero if no discrete
or auxilliary dimensions are needed.
To specify filenames for subsequent dvegas_save and dvegas_restore calls, use
CALL dvegas_init_files(outputfile, inputfile)
The defaults are dvegas.dat-new and dvegas.dat, respectively.
After dVegas has been initialized one can build iteration blocks as follows:
DO itn = 1, iterats
[...]
CALL dvegas_before_iteration(ncall, init)
DO i = 1, ncall
CALL dvegas_get(x, k, wgt, aux)
f(1) = fxn(x, k, aux)
CALL dvegas_put(f)
END DO
CALL dvegas_after_iteration()
END DO
CALL dvegas_result(integral, error, redchi2)
[or CALL dvegas_acc_result(integral, error, redchi2)]
[...]
CALL dvegas_end()
Subroutine argument types:
ncall: INTEGER*8 (!)
init: INTEGER, explained below
x, aux, f: DOUBLE PRECISION arrays, random variables and function values, explained below
k: INTEGER array, random variables, explained below
wgt: DOUBLE PRECISION, VEGAS weight
integral, error, redchi2: DOUBLE PRECISION arrays, pass accumulated results, errors and "reduced chi^2" = chi^2/(nr of estimates-1)
The parameter init has the same meaning as in the original VEGAS program:
init = 0: fresh start (no prior information)
init = 1: use previous adaptation info, but discard shot data
init = 2: use previous adaptation info and build on previous shots data
The get and put subroutines retrieve the random variables and return the function values:
CALL dvegas_get(x, k, wgt, aux)
CALL dvegas_put(f)
The array x contains the random numbers for the continuous dimensions,
the array k contains the random numbers for the discrete dimensions,
the array aux contains the auxilliary random numbers. The VEGAS weight (For
each iteration the sum of the weights is normalized to 1.)
is passed in variable wgt, so that it can be used, for example, for
routines that fill histograms. The array f is expected to be filled
with the function values that correspond to the arguments
x(1), x(2), ..., k(1), k(2), ... and aux(1), aux(2), ...
The value of the first function is assigned to f(1). The integration is optimized
for this function according to the VEGAS algorithm. Secondary functions can simultaneously be
integrated and their values are expected to be passed as f(2), f(3), etc.
Persistence: The status of dVegas can be saved and restored later with these two calls:
CALL dvegas_save()
CALL dvegas_restore()
dvegas_save saves the dVegas status (init. parameters and weights) to file dvegas.dat-new.
dvegas_restore restores a previously saved status from file dvegas.dat. Other filenames
can be specified with dvegas_init_files, see above.
PROGRAM TESTDVEGAS
IMPLICIT NONE
INTEGER nvals(1)
INTEGER*8 ncalls
DOUBLE PRECISION integral(1), error(1), redchi2(1)
CHARACTER*20 outputfile, inputfile
INTEGER i, init, itn
DOUBLE PRECISION fxn2
DOUBLE PRECISION x(5), wgt, aux(1), f(1)
INTEGER k(1)
nvals(1) = 100
ncalls = 100000
CALL dvegas_init(5, 0, 50, 1, 0, nvals, 1)
* 2 iterations with init = 0, no initial adaptation
DO itn = 1, 2
init=1
IF (i .EQ. 1) init=0
CALL dvegas_before_iteration(ncalls, init)
DO i = 1, ncalls
CALL dvegas_get(x, k, wgt, aux)
f(1) = fxn2(x, k)
CALL dvegas_put(f)
END DO
CALL dvegas_after_iteration()
END DO
CALL dvegas_result(integral, error, redchi2)
* CALL dvegas_acc_result(integral, error, redchi2)
* result is 199.0
WRITE(*,*) 'Result:'
WRITE(*,*) 'integral = ', integral(1)
WRITE(*,*) 'error = ', error(1)
WRITE(*,*) 'chi2/itn = ', redchi2(1)
WRITE(*,*) ' '
* now 1 iteration with init = 1, use previous adaptation
* save grid
CALL dvegas_before_iteration(ncalls, 1)
DO i = 1, ncalls
CALL dvegas_get(x, k, wgt, aux)
f(1) = fxn2(x, k)
CALL dvegas_put(f)
END DO
CALL dvegas_after_iteration()
CALL dvegas_save()
* now another run after restoring adaptation from file
outputfile = 'run2.out'
inputfile = 'run2.in'
CALL dvegas_init_files(outputfile, inputfile)
CALL dvegas_restore()
CALL dvegas_before_iteration(ncalls, 1)
DO i = 1, ncalls
CALL dvegas_get(x, k, wgt, aux)
f(1) = fxn2(x, k)
CALL dvegas_put(f)
END DO
CALL dvegas_after_iteration()
CALL dvegas_save()
CALL dvegas_end()
END
DOUBLE PRECISION FUNCTION fxn2(x, k)
IMPLICIT NONE
DOUBLE PRECISION x(5)
INTEGER k(1)
IF (k(1) .LT. 99) THEN
fxn2 = 1 * x(1)*x(2)*x(3)*x(4)*x(5)* 2.0**5
ELSE IF (k(1) .EQ. 99) THEN
fxn2 = 100 * x(1)*x(2)*x(3)*x(4)*x(5)* 2.0**5
END IF
END