A multi-dimensional integration package for a vector processor

Computer Physics Communications 48 (1988) 353—365 North-Holland, Amsterdam

353

A MULTI-DIMENSIONAL INTEGRA11ON PACKAGE FOR A VECTOR PROCESSOR S. KAWABATA National Laboratory for High Energy Physics (KEK), Tsukuba, Ibaraki, Japan

and T. KANEKO Meiji-Gakuin University, Totsuka-ku, Yokohama 244, Japan Received 15 September 1987

A program package, VBASES, a supercomputer version of multi-dimensional integration package BASES, has been developed. It was applied to the calculation of cross section for the radiative Bhabha process. The execution time of the integration was shortened by a factor of 80 in comparison with that by a scalar computer.

PROGRAM SUMMARY Title of program: VBASES

No. of lines in combined program and test deck: 2474

Catalogue number: ABBU

Keywords: multidimensional numerical integration, vector processor, supercomputer, BASES

Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue) Computer for which the program is designed and others on which it is operable: FACOM, HITAC, CRAY Computer: HITAC-S810; Installation: National Laboratory for High Energy Physics (KEK), Tsukuba, Ibaraki, Japan Operating system: VOS3 (HITAC) Programming language used: FORTRAN77 High speed storage required: 250 Kwords No. of bits vi a word: 32 Peripherals used: terminal or card reader for input, magnetic disc and line printer for output

Nature of physical problem In high energy physics the numerical integration package is used to obtain the cross section for the physical process. When the physical process is complicated, the integrand may become a very long FORTRAN source code. In such a case, the integration itself requires too much CPU time for practical use. Method of solution We have developed a program package VBASES, a supercomputer version of multi-dimensional integration package BASES. By using YBASES, we could shorten the CPU time by a factor of 80 for the integration for the radiative Bhabha scattering. Typical running time The running time depends especially on the length of the integrand subprogram which gives the differential cross section for the process. If we take the process e~e —+ (Z°) —~ vi’y as an example, the integration takes 1.37 s on the HITAC S810 Model 10.

OO1O-4655/88/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

354

S. Kawabata, T. Kaneko

/ A multi-dimensional integration package

LONG WRITE-UP 1. Introduction In high energy physics, when physical results from experiments are compared with theoretical prediction, it is necessary to calculate the cross section of the physical process. A numerical calculation of cross section is made by the following way. First we write Feynman amplitudes of the process and take the square of the sum of them. Traces of gamma matrices in square of amplitudes are dealt with by a symbolic manipulation system REDUCE [1]. REDUCE generates FORTRAN source code of the differential cross section. This differential cross section is integrated over the phase space by BASES [2]. This method is powerful to obtain the cross section for many physical processes [3]. However, a large amount of CPU time is necessary for the numerical integration, when a very lengthy code is generated by REDUCE. There are two reasons for such a large CPU time. The first reason is that the integrand has a huge number of arithmetic operations and is calculated at many sampling points in the phase space. When some calculations are performed repeatedly, a vector computer wields its power. Since, in our problem, this is just the case, it is possible to shorten the execution time of the integration by using the vector computer. In order to make program fast with the vector computer, one has to write a vectorizable program which is interpreted into vec-

package SPROC [4] divides the source code into small subroutines and rewrites it into a vectorized form which compilers are able to optimize easily. In section 4, we present results of the performance test of VBASES on several vector cornputers as a conclusion.

2. Vectorization method of integration We describe briefly the algorithm of the multidimensional integration used in BASES [2]. Let us consider one-dimensional integration for simplicity, 1

I

J f(x) dx.

=

—

A cumulative estimate S of I is given by i

~2

~=

a2 a=1

(4)

Ga

and 1

—

,~,

—

~

1 a

where a2 is the variance of the estimate S in the ath iteration and n is the number of iterations. S is calculated by a —

1

—

N

N~] 1

tor operations by a compiler. A vectorization method of our integration algorithm is described in section 2. We have developed a multi-dimensional integration package VBASES for that purpose. It requires users to prepare subprograms USERIN and VBFNCT for the initialization and integrand, respectively. Their specifications are given in the section 3. Another reason for the large CPU time is the problem of FORTRAN compilers. For a very long source code, many FORTRAN compilers give up high level optimization because it requires large work space and large CPU time. In some cases the size of source code generated by REDUCE may exceed the capacity of the compilers. A program

(3)

0

Sa

•~

f( x,)

6

PaXj’

where pa(x) is a probability density on the integration volume and Ncaii is the number of sampling points in one iteration. The density pa(X) is modified so as to minimize a~for each iteration of the estimation. In the first estimation, uniformly distributed random points are chosen (i.e., in determining Si). The information about the behavior of f( x) obtained in this first sampling is used to define a new density Pa(x) which reduces a~in the next iteration of the estimation for S2. After each subsequent iteration, the density PaCX) is again refined for the next one. The right hand side of eq. (6) can be easily vectorized, since the

S. Kawabata, T. Kaneko

/ A multi-dimensional integration package

calculation of the function f is performed at ~a1l sampling points which are independent of each other. The integration volume is divided into small hypercubes. A set of sampling points for an iteration is constructed by selecting several points from each hypercube. The number of sampling points of one iteration ~1~a11 is the product of the number of hypercubes NHCUBE and the number of sampling points in a hypercube NTRIAL. We take NHCUBE as the size of the vectorized DO-loop. We show the structures of BASES and VBASES in figs. la and b, respectively. In BASES, the calculation of the integrand FUNC, which returns the result for a single sampling point, is repeated by three DO-loops. On the other hand, one of the DO-loops, whose length is NHCUBE, does not appear in fig. ib, and the subroutine VBFNCT returns NHCUBE results of the calculation of the

355

The performance of VBASES depends directly on the performance of the subroutine VBFNCT. There is a useful technique to make a vectorized program faster. In calculating the integrand, some sampling points do not satisfy the kinematical condition and the integrand is set to zero. If many such points exist, it is more efficient to omit then from the calculation. For this purpose we take the gathering and scattering method. A performance test showed that this method made our program run faster by about factor of two. For vectorization of a program, one should note the following points. Consider an addition operator A ~ which operates on two n-dimensional vectors x, y and returns an n-dimensional vector z whose ith element has the value of the sum of the ith elements of x and y: Z

=A~,(x,y),

where z.=x 1+y1,

integrand. The subroutine VBFNCT, which is prepared by users, should calculate the integrand NHCUBE times by a vectorized code.

a

start CALL USER1N initialization of BASES DO [TB

=

DO NP

1. (number of iterations) =

DO NT

1, NHCUBE =

The parallel operation for each element does not always give the correct result even though the code is vectorized. In the calculation of x1 = a + x,2, for example, x3 should be determined before x5 = a + x3 is calculated. Thus there should be no mixing between input and output data is a vectorized program.

numbers for one point

3. How to use VBASES

calculate the tntegrand FUNC

We take the process e~e (Z°) viy as an example throughout this paper. This example is used in the description of BASES [2]. Since the usage of YBASES is almost the same as that of

ENDDO

—~

ENDDO esttmate the tntegratton ENDDO

b

start CALL USERIN initialization of VBASES =

(2)

1, NTRIAL

sample a set of random

DO [TB

for i = 1,..., n.

1, (number of tterattons)

DO NT 1, NTRIAL sample a set of random numbers for NHCUBE potnts CALL VBFNCT ENDDO estimate the integration ENDDO

Fig. 1. (a) Structure of BASES; (b) structure of VBASES.

—*

BASES, its detailed description is omitted here. The structure of VBASES is shown in fig. 2. In order to use VBASES, one should prepare the main program VMAIN, subroutines VBFNCT and USERIN, which are enclosed with ((and)) in fig. 2. VBMAIN called by VMAIN controls the job flow. To steer the job, VBMAIN reads several parameters from unit 5. They are the loop control parameters Li and L2, print option flag NPRINT, job control flag IFLAG and CPU time limit CTIME of the job. At the beginning of the integration job, IFLAG should be set equal to zero.

356
S. Kawabata, T. Kaneko VMAIN

/ A multi-dimensional integration package

>>

subroutine is identical to that of BASES. We show an example in fig. 3. One has to set at least values of parameters XL, XU, NDIM, NTRIAL in the common block /BASEi/. Arrays XL and XU specify the lower and upper limits of the integration volume in each axis. NDIM should be equal to the dimensionality of the integration and set equal to the number MDIM, defined in subsection 3.3. The reason why we have two variable NDIM and MDIM for the dimensionality is not to change the specification in USERIN from that in BASES. NTRIAL is the number of sampling points in a hypercube and its default value is two. One can set here the values of kinematical parameters which are used in the calculation of the integrand. Variables CPnn’s and Cnnn’s appearing in fig. 3 are .

H---- VBMAIN F---

vwri~ia I---

BUINIT USERIN

>>

H-- BSREAD F—--- VBASES VBINTO H-- VBTIME I

---

~ VB~NDMvB~ F----

VBTRNS

F--<< VBFNCT XHFILL I--- DHFILL I

F--- VBTERM

~ii= H—--— XHPLOT XHRNGE F--- XHORDR ‘---

F

I

---

XHSCLE

DHPLOT

H-- VBPGRD H-- BSWR!T

powers and products of constants, which are generated by SPROC for optimization. One can obtam histograms and scatter plots as optional outputs in the same way as in BASES. Their initializations should be made here by calling subroutines XHINIT and DHINIT for histogram and scatter plot, respectively.

Fig. 2. Structure of VBASES.

3.2. Subroutine VBFNCT If the integration is not finished in the CPU time limit, the job is terminated after all information is saved in a file of the unit number 23. To continue the integration by the subsequent job, IFLAG is to be set to the number which is given at the end of the previous job. In the following two subsections, we describe specification of subroutines USERIN and VBFNCT. While USERIN is not necessarily required to be vectorized, VBFNCT is. It is almost impossible to rewrite the integrand to a vectorized code, however, especially when a very lengthy code is generated by REDUCE. Therefore a program package SPROC [4] translates output of REDUCE into a vectorized code and produces VBFNCT and USERIN as well. We used this package to make these subroutines in our example. 3.1. Subroutine USERIN Here were give oniy some comments on subroutine USERIN, since the specification of this

Since the calculation of VBFNCT is the most time consuming part in VBASES, it must be vectorized. An example of VBFNCT is shown in fig. 4. In FORTRAN vector operations are expressed as inner DO-loops which satisfy several conditions, e.g. not including subroutines or external functions, no dependence among input and output data which breaks parallelism of vector operations, etc. These conditions depend on the machine. VBASES generates a set of sampling points in the MDIM-dimensional volume and passes them to VBFNCT through an array X(NHCUBE, MDIM) in the common block /VBVECT/. Parameters NHCUBE and MDIM are the number of hypercubes and the dimensionally of the integration volume, respectively. The values of the integrand corresponding to sampling points must be returned through the array FX(NHCUBE). Arrays PWnn and Tnnnn are powers and products of variables, which are generated by SPROC in order to optimize the original code. The final

S. Kawabata, I

Kaneko / A multi-dimensional integration package

351

C

—

0.

C

~0.

—

.0

-~

0 — 11.0

C 0 141 0 11:0 0.00 .4-.— 5.

a. U .a --‘CC

OlaIn C ——~ — .400* .-......0 04004:. 1104*04

04 4: • I. O

040040t.0..l000000 • .4 441 (450 0 0 0 0 0 0 0 0 00 000 00000000 444 - 4:. 0. . . 0. 4:. 0. 0 0 0 0 00 0 U * CO Cl) .44441 0 U Cl) 4440 00 U .4 0 U 4444440 ..I.-(14104.C... * * * (‘5 (. 04 44104(44(4144.104 000044. 141 (‘5 0 0 0 0 0 0 .C441000000000000000000000000000040 ,. = = a * 0 * 0.0. 0. (1. 0. 0. 0. 0. 0. 0. 0. 0. 4:. 4:. 4:. — 0 0 0 0 0 U044:4CO04N04UUUUU0000000000U0000U00000UUCO (I

4(00

.4.4 0 0 . .00 a 0 .4 0 0 U U 0 40 t- * * .4 0 0

. 144 * ‘0 0

0 I .04 04 * * (14 t . 0 (.41

0 0 0 0. 0.. 0 0

tIIIlIIIIIOIIIt0000flhIOIlllIINlllllItlIIIlflflU(IuI

..100401n4044000.40fl0In0t..O0 0.C1410401n 40t..0O.0 (44 0 In (0 44.. 04 0 0 0 0 0 0 0 0 000.4..I ..l.4..I .-. .0 .40040004040404 444 (44 (‘5 000000000000000000000000000000000000 0.0.0.0.0.0.00000000000000000000000000000 U000000000U00000000000U000UUOUUUUUUU

—

C C -~ .0 0

.4.4)l.)& IC IC IC 0

.11 *4 *

—

- - . 4:. .411..I -

C 0 .

4: 0 0 .4: 0 0 0.

— — — (.41.14.41

.4 .0 0 —~ •1 -.

—

LU

* 10)00 .~ .3 .4 ..I..3..~ -CCC 000

0 0

.0 *4

=

~0 (lIZ 0404

U

0 C 0

*4 *

o

0. C —

.4 ‘41—

0 44400000fl

fl

In IflInIn

In In

flInInlO

(0 fl4040’0400

0(4-

—

401’-

V

—

~ 44. (.4

(.444. B

N

Z

0

0

0.0(44

40 — 0.

004

04 04 U 0(40 0 14 - =

U

000 000 000

~

01

H II C H I: H 0 011 II — H H 04 :4 04:3 I: (/4 II 0:: II H *4 II (I 04 0 0 — 0 11.1* 0 0 11 0 0 II II 04 fl II 04 11 II 0 (I (/1 0 IIII O1

0. U

0.-I 141(44 0000 0000

. 0 0 ~0

In 0 4:. U

In 145 .40 0.4(44(40 0000 0000 000U

. 04 0 (.4

.4 0 4:.

U

0

0(4000 0 .4(44444 0 0 0 0 0000 U 000 (4: .4 Cl I~00.004%. 0 0 0 0 0000 000U

00 . . (4-444 fl~.. .4 ‘.-04 1114 0.. Z Ia (12

U

. 04 C 0 N Cd,

= — — —

0 ZO - 5< — = 0 4.11 .4*( N — : 0404 0 Ic U .4 U

-

04 C 1411. *

=4

0 COO — — U Ia U 0 04 C *

UUU

‘aooo

n 0 0. U

-

000U

141 04 4:1 0 04

C .1101 Dl1/4 04 (/5 04 40

* 0 (—I

C C 40004

0 U

0040 — 00 Ia 04= 0. = 04 04CC —00

0 0 04 04 0 U

(4* (4400110 0 0 .0.-.(.J S.11. *000 0000 0 0 U 00 U 0 3 .0- a 0 .0 0 (4.. 145 040.00.00...r.... 04 044.) 0.0 0 0 0 0 0 0* 0 04 0000 • .0 0 * 04 0 04 = 0 0 0 0

00

<~ .4 - =.0 0 0 01.4 4004 (41 0 00 114 04 00 .0 . (‘40 04 — (04 0.0 C C 1.4 1.~ .14< 00

(44

z 04 (14 0 *

14..-. 0441 (/5 * * * 0 C

(/5 = 4: 4:. .1.4 (41 .14 .4 * In N .0 0404 0* 04 0 04 * I 14 4 N . U 0 0. 0 (*4 14 * 0 0 3... — Ia Or 04 ...--OC04 *0004 <.4143 — - * U (11 In 3 In In - .0 ~. * * 0(11 115 — Z (0 04 0 0 * Cl 0 t~(1- - — A. 04 1111 = 04 (/5 = (44 0. (‘4. :14 (.4 0 II U H II * H 111(11 II (I C (I (I II 4441 04— a 0 40 *

I. 4: 0 4: 0 5. C 0. 1/4 4:4 1/4 C 04 04 II .4 0

C C 0.044404 040Ia<04014C>< 04CO044013314.1404140404000II.

0 — 114 0 . (‘4 In 0 0 04

(I

II

0

in 0 0 4 0 . a . * In 0 II — — 00

0

03 — * = It: 04 (C 0-

II II

0

—a0—04004= 0 .4 *104 = U II’ ~ —C ~

0

— — 40 0 CC 04 0 * (
—

to to 00 0 0 —

c~ , 0 0 0 0 040

II II

(‘4 * 0 4: * C (I4 U H (I 0 (‘4 00 11.0. 00

358

S. Kawabata, I Kaneko

/ A multi-dimensional integration package

04004 I0 00.4—004 000*004 000.4 I — O .4 U 004 * (144*0.04.4 0.0.
— ... — — - *00— a — — * 0.0.-. 0...U000 (*_t.I :0 * 140 .4~0 *SCOCC-*0 040...0.1 — 0.101 0*..0*.__

~ 11 U 0 0 11 — 0 (0*U0.~..0 0

~ 0 OCCPJ(-. 0 0000.0 0 ( %. 0 (‘1 — (.4

C 44

-000 .4 H II 0 0 — — — — — (0.41 0 0 0 111 0 0

114 0 * 0

—

1(1 * — 0 (.1

~C

(1) 0 Z (-I 04 0 U

—00+0_a 0001.1.1(010 U0—0* .-. -( 0 — I Oct ..(0I0+ 0000+0l .4 0.00l11~*4.( ~.U*l.41...~1100 fl*0+.*004:. *CCU04O..— 0. 0 I 0 . 0 — (44 0000001. * 0.4.0 ~00 * 00000414*— 0_al 00O(44 0.*UP4O0400 0.. 4. 000 *0. —*1(0*UO~U

In * -.000 0.0 — I 000 0 U 00004 * (-.1.. * ..3...OIC-*O04 IC * =.~..1*1.4la0

__100000440 0.1/4 a : ..e I 00 0*. * U — * (4400. 0...~. * * 000 U Ia.-0&.-~U*+ 1< 1/4 01 (41 — I 0(0 — *001040(44*

U

(14 0 0 0 ((4 04 *4 * *4 1. C (4

.404.0In0104 0 0(0010040*1 0.0 * — I 001.1 04 I H. - (44 0 U ~ 0 40 — 0 — a U — 0 014:~+*U

— .4 0 0 I (14

.4 (41 * 11

4

*0000*1C —000.0_a 00000010 *00... (-. 0 U .)l44.:;:(. *001(.4 .00*10(0 W00..~....* .40004 444 .-.U I......0a_.l. (‘I (H * — 04 .4 —. 004.0.00 00 U_0_0O 000041.U .0001 03 I 0.-.04.* ...0....O.-.0t4

..(.tM0.00*( (/3 .4 0 * 0 0 II 0 41. — 00 (I 0(41.4. 4.0 .4 U 0 * 0 (0..0400 0 — .4 U 0 0 *

0 0. U * (C N

1f~:;~

;~~:~I

— — — 0.0440-. 0)0)5 — 0 * * 40 IaIn,0* IC — (0 =~0~

— —

04 14 — .0 1/3 — — * 4:4 — . — *1 — SC IC

*4

— 040 0— 00 UC

0

Z((.

— I.

.0*1 ..(I.. U I.. SI CC I.. I II

04

.*~U0CI0-.a~

00 .4 (0 + — U (0 I U II 0(44004 .4 0111 .4 —0 .
—

—

(44 0 0 1.1 0 0

—( 0 C -l Cl 0

0 II (11 (11 — — 0.4 04

0 .o 11 U

— — (.1 (4. (C 04 0 0 II. 0 0 U

—

0. 0

C

(4. 0

— (C

C —

0. 0

— (*4

(11 — 1<

(45 0

SC 4*

.41 — 15

SC IC

14 4: 0 5. *1 11 4. II

.0

(~4

* .4

— Ia Ia — 0.

— Ia .4 (4.

— Ia .4 l~4 (4.

0

1<

SC

0

4:4 C I .4

Ia .4 C U

Ia .4 C U

Ia Ia C U

.4

H Dl 0 .11 U 0 4:.

ILl

1 3)

04 0 0 (.1 0 (04 04 0(14

* .0 *4 *

E

0

40 5<

S. Kawabata, I. Kaneko

/ A multi-dimensional integration package

result should be stored into the array FX. Histograms and scatter plots are filled here by calling subroutines XHFILL and DHFILL, respectively. The arguments of these routines have the same meaning as in BASES, but they should be arrays in VBASES. 3.3. Common block VB VECT As mentioned in the previous subsection, the integration variables and resultant values of the integrand at all sampling points are contained in the common block /VBVECT/. Because of this, the program size is mainly determined by the size of these arrays. In order to keep the program size as small as possible, we must define these size parameters by the PARAMETER statement. The sizes of arrays are determined by the dimensionality of the integration MDIM. The number of hypercubes NHCUBE is given in table 1 as a function of MDIM in the same way as in BASES. One must change the values of these parameters in the PARAMETER statement of the subroutines VBINIT, VBITRA, VBRNDM, VRAN, VRAND, VBTRNS, VBTERM, XHFILL and DHFILL, which use common block /VBVECT/. The definitions of Ng and Nd in table 1 are given in ref. [2]. 3.4. Machine dependent subprograms There are machine dependent subprograms in VBASES. Before using VBASES on the user’s machine, it is necessary to check whether the following routines are suited for the machine, Table 1 Relation between MDIM (NDIM) and NHCUBE MDIM

NHCUBE

359

(a) Subroutines VBRNDM, VRAND and VRAN Subroutine VBRNDM generates NHCUBEMDIM uniform random numbers between 0 and 1 and returns them through array X(NHCUBE, MDIM). Subroutines VRAND and VRAN are auxiliary ones for VBRNDM. Since the method of random number generation uses the length of the word of the machine and the data format of floating point numbers, one should be careful to use subroutine VBRNDM on user’s machine. (b) Subroutine VBTIME Subroutine VBTIME measures CPU time by using system clock routine. When the value of an argument IFLG is 0, the CPU timer should be started. If IFLG has other values, this subroutines should return used CPU time in second by the argument TIME.

4. Conclusion We tested performance of our system by calculating the cross section of the e ±e e + e y reaction in electro-weak theory, because the FORTRAN code of its integrand is long enough to test the performance of VBASES. The REDUCE output of the differencial coross section has 2655 lines including 49441 arithmetic operations. After being preprocessed by SPROC, the source code is divided into 25 subroutines. The number of arithmetic operations is reduced to 14188 by the optimization. The optimization is performed with the constraint of available main memory less than 13 Mbytes. The execution times on vector computers are shown as follows: — —

Computer

CPU time

—

Ratio of acceleration

Original FACOM M382 3192.63 1.00 Vector FACOM VP5O 115.05 27.75 FACOM VP100 75.20 42.46 FACOM VP200 38.52 82.88

Ng

Nd

1 2

25 25

50 50

25 625

6 7 8

4 3 2

48 48 50

4096 2187 256

where CPU time corresponds to the GO step. For comparison, we also tested BASES on a scalar computer FACOM M382 with the same integrand

I

whose code isThe the ratio original form of the output source by REDUCE. of acceleration is

360

S. Kawabata, I Kaneko

/ A multi.dimensional integration package

the relative CPU time on each machine to that on FACOM M382. The same program was tested on the vector computers HITAC S810-10, S810-20 and CRAY-XMP. We obtained expected performance on each vector computer according to its hardware.

Acknowledgement We thank Professor Y. Shimizu for valuable discussions and continuous encouragement through this work. Thanks are also due to Dr. K. Tobimatsu for helping us make the test programs. We are indebted to Professors R. Kajikawa and Y. Oyanagi, Dr. Y. Fukui and members of KEK computer center for their hospitality to test our program on vector computers.

Appendix. Program structure and its components The program structure of the integration package VBASES is shown in fig. 2. The program components, which concern users, are briefly described in the following. The routines marked by * * should be prepared by the user. MAIN Routine VMAIN** Purpose: To call subroutine VBMAIN. Subroutine VBMAIN Purpose: To control the flow of program package VBASES according to the job control flag IFLAG and to set the VBASES parameters to default values. Subroutine USERIN * * Purpose: To initialize parameters in a common block /BASE1/ and for calculation of the integrand and to define histograms by XHINIT and DHINIT. Comment: This routine is called by VBMAIN at the beginning of the integration job. Subroutine VBRNDM Purpose: Interface routine between VBASES and random number generator VRAN.

Comment: the generated NHCUBE * NDIM random numbers are stored in the array X(NHCUBE,NDIM) of the common block /VBVECT/. When the random number generator VRAN is not suited for the user’s machine, it should be replaced by a suitable one. Subroutine VBTRNS Purpose: To transform the uniform random numbers into the integration variables between XL(i) and XU(i), where i = 1 NDIM. Comment: The following arrays in the common block /VBVECT/ are changed. X(i, j) Transformed random numbers between XL and XU are given as outputs, where the distribution of random number X(i, j) is to be matched to the probability density p(x~). WGT( i) Resultant probability density p (x) at the sampling point X(i). =

=

Subroutine VBFNCT * * Purpose: To calculate the integrand. Comment: This routine must be prepared by user. If histograms are required, histogram-filling routines XHFILL and DHFILL should be called here. The resultant values of the integrand should be filled in the array FX( i), which corresponds to the independent variables X(i, j= 1, NDIM). Subroutine XHINIT(ID,XMIN,XMAX,NBIN, TITLE) Purpose: To initialize a histogram. Comment: This routine should be called in USERIN, if a histogram is required. ID Histogram identification number. XMIN = Lower boundary of the variable. XMAX Upper boundary of the variable. NBIN = Number of bins of the histogram (max. 50). TITLE Title of the histogram, which is a character string up to 62 characters. =

=

=

=

Subroutine D H IN IT(ID ,XM IN ,XM AX, NXBIN,YMIN,YMAX,NYBIN,TITLE) Purpose: To initialize a scatter plot. Comment: This routine should be called in USERIN, if a scatter plot is required.

S. Kawabata, I Kaneko

/ A multi-dimensional integration package

ID = Scatter plot identification number. XMIN Lower boundary for the horizontal axis. XMAX Upper boundary for the horizontal axis. NXBIN = Number of bins for the horizontal axis (max. 50). YMIN Lower boundary for the vertical axis. YMAX Upper boundary for the vertical axis. NYBIN Number of bins for the vertical axis (max. 50). TITLE Title of the plot, which is a character string up to 62 characters. =

=

=

=

=

=

=

=

Subroutine XHFILL (ID, X, FX) Purpose: To fill the IDth histogram. Comment: This routine should be called in VBFNCT. Before calling this routine, the routine XHINIT should be called with the same ID number. ID Histogram identification number. X(i) Variable array. FX(i) Weight array of the points X(i), e.g. values of the integrand at X(i). =

=

=

361

XU( i): Upper boundary of i th integration variable. NDIM: Dimension of the integration (up to 10). NOCUB: Maximum number of hypercubes (default and maximum numbers are 6000 and 17000, respectively). NTRIAL: Number of sampling points per hypercube (default is 2). ACC1: Required accuracy in the grid defining step (default value is 0.2%). ACC2: Required accuracy in the accumulation step (default value is 0.01%). ITMX1: Maximum number of iterations in the grid defining step (default number is 15). ITMX2: Maximum number of iterations in the accumulation step (default number is 100). At least XL(i), XU(i) and NDIM should be given in the routine USERIN. The functions of all program components of VBASES are summarized below. A more detailed description is found in ref. [2]. Components, prepared by user himself

Subroutine DHFILL (ID, X, Y, FX) Purpose: To fill the IDth scatter plot. Comment: This routine should be called in VBFNCT program. Before calling this routine, the routine DHINIT should be called with the same ID number. ID Plot identification number. X( i) Variable array for the horizontal axis. Y(i) Variable array for the vertical axis. FX(i) Weight array of the points (X(i), Y(i)), e.g. values of the integrand at (X(i), Y(i)). =

VMAIN

: Main routine in which VBMAIN is called. USERIN : Initialization of user parameters. VBFNCT : Calculation of the integrand. Components, only for VBASES

= =

=

In order to prepare the routine USERIN, it is necessary to know the contents of the common blocks /BASE1/ and /BASE2/. Some of these parameters for VBASES should be set by user and transferred to the routine VBASES. The contents of /BASE1/ and /BASE2/ are as follows: COMMON/BASE1/ XL(10),XU(10),NDIM,NOCUB,NTRIAL COMMON/BASE1/ ACC1, ACC2,ITMX1,ITMX2 XL( i): Lower boundary of i th integration variable.

VBMAIN VBTIME VBASES VBINTO VBINIT VBITRA VBRNDM

: : : : : : :

VRAN VBTRNS

: :

VBTERM : VBGDEF : VBPGRD : XHFILL DHFILL

: :

Steering routine of VBASES. Interface to system clock routine. Main integration package. Initialization subprogram no. 1. Initialization subprogram no. 2. Control part of the iteration loop. Interface to random number generator. Random number generator. To transform the random numbers into integration variables. Loop termination for one iteration. re-definition of grids. To print the grids for all integration variables. To fill histogram. To fill scatter plot.

362

5. Kawabata, I Kaneko

/ A multi-dimensional integration package

Components, in common use with BASES

References

BSUSRI

[1] A.C. Beam, REDUCE User’s Manual, version 3.2 (The Rand Corporation, 1985). [2] S. Kawabata, Comput. Phys. Comniun. 41 (1986) 127. [3] K. Tobiniatsu and Y. Shimizu, Prog. Theor. Phys. 74(1985)

: To check the VBASES parameters set in USERIN. BSREAD : To read the VBASES data from unit 23. BSWRIT : To write the VBASES data on unit 23. BHINIT : To reset histograms and scatter plots (for loop option). BHRSET : To reset histograms and scatter plots (before accumulation). BHPLOT : Control part for printing histograms and plots. XHINIT : Initialization of histogram. DHINIT : Initialization of scatter plot. XHPLOT : To print histograms. DHPLOT To print scatter plots.

567. J. Fujimoto and M. Igarashi, Prog. Theor. Phys. 74 (1985) 791. M. Katuya, J. Morishita, T. Munebisa and Y. Shimizu; Prog. Theor. Phys. 75(1986) 92. K. Tobimatsu and Y. Shimizu, Prog. Theor. Phys. 75 (1986)

~os. M. Igarashi and N. Nakazawa, Nuci. Phys. B288 (1987) 301. J. Fujimoto, M. Igamashi and Y. Shimizu, Prog. Theor. Phys. 77 (1987) 118. S. Kuroda, T. Kamitarn, K. Tobimatsu, S. Kawabata and Y. Shimizu, Comput. Phys. Commun. 48 (1988) 335. [4] S. Kawabata and T. Kaneko, to be published.

S. Kawabata, I Kaneko / A multi-dimensional integration package +1.4

+

01 CI /451* I 01* ID II

I I I I

I * ID

I I

1141* —1*4 04/4:4:

I

OI***

I

0/4*4:4:

I

I

*******

01

o o

I

I

I **********4********

I

I I I I /

I I

**********4*************************** ******* *4* * *4* ********* **************** *4*4*4:4: *4* ***************************** * **** ** ** * * *****#**************** ******4 *4: * * *4* *4: ***********************

.or_I.J.oo.n0004awIo..0e,oI..~Ia,11rl4oorn.N.-oo.0o..oIn-t~~~N.--o I O0~0401.l0010100044.~ 1O(0(01040(O4OIflIflUI01 I (~10114000000000000000000000000000000i

U 0 C

I 0*4:4:4:41 I

I ItSIflIII11101145 (1101 I

I 000000000000000000000000000000000000000

I

I (001(110401.0(004040111 -. .- (I1~0 54501(45400404 ~ (~ (11(00*4(45l0*0140001010 I 0011400404*00.01-01., 0(11 ~ (4410.~ 001l4-01~OS ~040*0 ,- 040101* I OOO 04O01(nO4(flW101040401040404W4W’Pfl04010404004*11fl~1C’l04040104040400404

I I I

14 ,000000000000000000000000000000000000000:

144 I 000000000000000000000000000000000000000 I C I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I WWWWWUWWWWWUWUWWWWWWWWWI4JUUIQISJWUUWWI4JUWUIWW I — I .— 0* (1) ((1 (4I111 III 44.4 (14 0 ,.I,—,—~004.— 154 (‘I 0 ((I 1- -4 1I (44 (14 (1) .44 .1 .._ 1(1 .1*11 (‘4 (14 .— 154 ‘0 I I *44*4.4.4.4*4*4 I 00 0. U CI I II/l

*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4*4

000000000000

000

0000

000000

00000

I 000

0000

I 115*4t.t.t.e*4.**4.t.*.*.*4**.*.**4.**4 I 000000000000000000000000000000000000000 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I UUJWWWUWWUUWUUWUUWUWILJLIJUUISJISIUWIIJWISJLLJUWUIULUWUW I 0001 04’0*4040.15J(...O0100I’JPflIfl0...*4Wl..0..(00144.040(0(C01*4C,l0’(0 I (11Pl.OJ0fllI0(140’40.t4.0*t.~ 04fll4-001001’.F..(0’OIflUCIfl*4.*04040404154 I . 0’*P~ .0(4101.4: -t -t -.4. .4:04040404 (‘1 040404(1104010100101010101010101010101010101

/

I

I 000000000000000000000000000000000000000 141 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I .4 I ULUWUIS4UWWUWWUUWIAJLIJWUWWWUUISJLUWI4JIUWWIUWWUUUWUU O I 0101 11%0’0O5C(J0I1104(l10 0404*4Cdr..’0.*OIflOCllst’0t..Olfl,...0.04 04 I .-I5J.O04N.e.401004’0.0400 0404.000040.WI.00.-t0.*0.00..t’O..r-N.--o* (42 I (00’004010140010100010000000..00 C140041h1(0111*4.*.t..*St*4 CI00400,..00010101.-OJ01010100lN0l44JC1lI%JCII15JCIIOJ15IC’4I’JOI

I

I I I I I I I

I II I I I I

:~:~: I 000000000000000000000000000000000000000

I

*11 *

I 000000000000000000000000000000000000000

:

z ~

I

~

~

. I 000000000000000000000000000000000000000 ‘I-I I- I 000000000000000000000000000000000000000

~

I I I

0Q0000000000000002200000000000020000000

O I ‘0010010*4~O04,00,..tr..00440040010004(0040101C~*41*.0401010..*1.ISA

E I 001114(1404

* ..1.4010151111t1(0’0’0P...

—

I 0000000000000000000000000

—

I I

0104.4:04(0 I~.00.0.-.

II. ~—r—*0001010404000.-.-.~

I -~0101010404(44

~

1’4(14 .4 (Cl ‘0 I~.O0. O,-0104*4In ‘0?’. 0010,(‘4 (14 .4: 1(4 .0 1-004 010101010101010101151 (11 ((1 (14 (11 (1404(44 (44((1 In

I I

I I

363

364

5. Kawabata, I Kaneko

/ A multi-dimensional integration package

+

I I 0/ 01 1111 1 0/ I I * I I I I I I 0 I CI *4 I I 01 I * I I * I I I (4. I .0 I 0 I 0 0 I *1041 I 001 o + I t I C I I I O I I IS I 00 I 4(40 I (‘4 I I ISO I O + C I I .4 I I I I I I 0 I 0 I — I I 0 * I I I I I I I I I I I 0 I .1 0+ (4* 01 .I O O C o 44 0 .0 0. * 0 0. o —

+

0 0 0 0 0 0 0 000 0 000 0 000 0 000 0 0000 0 0000 0 0000 0 00000 00 000 * 00 00 0000 * 00 000 0000 * 00 000000000 0 * 00 0000 00000 0 * 000 0000000000*000 0000000000 * 000 0000000000 * 000 000000000040000 000o000ooo*oooo 00000 00000 * 0 000 00000 00 000 * 0000 0000000000 4: 00000 000000000 0 * 0000 0 0000 00000 0 * 0000 00 0000000000*000000 00000 00 000 * 0000000 00000 00 000 * 0000000 0000000000 * 00400000 0000000000 * 0000000 oooooooooo.oooooooo * 0000 0 0 000 * 00000000 0 *eeooooooo*ooooooooo * 000000000 * 0000000000 * 0 0 0 0 0 0 0 0 0 4 -0 0 0 0 0 0 0 0 0 0 * 0000 0000 0 * 00000000000 *000000000*00000000000 * 000000000 4: 000000000000 * 000000000 * 000000000000 * 000000000 * 0000 0000000000 * 000000000 * 00000000000000 * 00000000 4: * 00000000000000 * 00000000 * * 000000000000000 * 00000000 * * * 00000000000000 * 00000000 * * * 000000000000000 * 00000000 4 * * 0000000000000000 * 00000000 4 * * 0000000000000000 *00000000***000000000000000000 *oooooooo***oooooooooooooooooe 4: 00000000 * * * 000000000000000000 * 00000000 * * * 0000000000000000000 * 00000000 * * * 0000000000000000000 * 00000000 * * * 000000000000000000000 * 00000000 * * * 0000000000000000000000 * 00000000 * * * 00400000000000000000000 * 00000000 * * * 00000000000000000000000 * 004000000 * * * 0000000000000000000000 * 00000000 * * * 00040000000000000040000000 * 00000000 * * * 0000000040000000000000000 * 0000000 * * * * 000000000000000000000000 * 0000000 * * * * 0000000000000000000000000 * 0000000 * * * * 0000000000000000000000000 * 0000000 * * * * 0000000000000000000000000 * 00400000 * * * * 0000000000000000000000000 * 000000 0 * * * * * 0000000000000000000000000 * 000000 * * * * * * 0000000000000000000000000 * * 00000 * * * * * * 00000000000000000000000000 * * 0000 * * * * * * * 00000000000000000000000000 * * * 000 * * * * * * * 00000000000000000000000000 * * * * * * * * * 4 * * * 00000000000000000000000000 * * * * * * * * * * * * * * 0000000000000000000000000 * * * * * * * 4: * * * * * * 0000000000000000000000000 * * * * * * * * * * * * * * 4: * 00000000000000000000000 ***************************************0

I OIM04IM*4..t*4440 04*40404,*. ?*.*.‘tInlfl*%0P1C4N.04.115*4’00. 011 .H10104004010’0000*0104*00*401000010001010400.0’00.-.-In0000 C I 1.-15. ‘001151 *4 000.0..0 ~* 1*. 0. .410101010104. ~1I I’J 0.1

I I 10 /0 I. /0 I I + I I I I 0 I I I 114 I 0 10 / — / / 0 I * 0 I I I >~ I 144 .0 10 I I / 112 CI I 044 I 0 15 1.-U I lOt I C + -— I I 0 I I ~ ~ I 0 C) I I In I U I 0 0 I 0. I . 40 I •~ I 0*( I I I I I I I I I I

04 0

01 0 .4.4 0 I 144 0 C

I I * 0 I I I I I I Ill I 0 I I I U I 0 I 0 I (. +0 4(41

N...0

I /0) I C

0.

151 (141

/0 III)

II

II

01 01 I +

10 10 I

-

I

I .-.-.-0101010101010404040404040404 .-* .4*4*4*4*4*4 ‘44 040101 t I 0000000 0000000000000000000000000000 ..I 1111111 111111111111111111? 1/1111111 I UUUUIUIUU IIJUWII.IUUWUUUUUUJIUUUWISJWU1WUUWUUUU (UI *01N.0.040010.*4’01-010*0404010)N.0.*00..N.’0N..0401.-0040’00*4041-010* El 040010101I-04*00401t(J.0,040.0In01*004N.01*40,’004*4r-040*0*04040*4 Oil .0400104N.01,-0*4N.0404N.*4H,0’0N.*400100’001’0040011-04N.1510000)*4’._IO040104InIn04*4’0.-01.t01,O014-N..t04.-4-4-p.01-t0401I-.-,-01.0~11II1I(4C4l.-Ifl(fl.-OI.-

I /

t I’. I IS E 01 I I I

I00000000000100000000000000000000000000000001 tI

It

041

~d0 0 44

(4 I ,r.J*4’0OI04Ifl.. 114 pC .4: .0 15- 00..0104 *4 ((4 P.. 00.001(14*4145 I 0.00104.*01N.00*1410,-*4N.O04~~O0

=

O I 004 2 I t ..-fll04(n*0 01CC .410 *

0 . P..

44.. 00401-04

.4: 01(0001

0

I

>

1.

04 *4 (44 f~..00..004 *4 1.- 00,01 .4 04*000400104 -2 HO 1. 001 4- I 4: 010101010101010404 ((4 0404040404 *4 *4 *4 .2 .14 .44 -4: *41(4 0 I 2 (410 01.4 +

S. Kawabata, I Kaneko

/ A multi-dimensional integration package

365

(4* ‘I.o 0

o

o (U

I~+ + — CII.’.’ 0.-I ((I...JI ‘ I fll*.l 1(101 0.1 0. . + I...

0 .0

I I

O 0

I I +

C 0 44 O .0 0. 0. o

I (4. 01

I I

O C O

I I I

C o 4.1 O .0 0.

.

•

111.4

I.. o I-

0.1-I .01 031..0 I 0.04040

0 —I,’ (‘4 (‘4 II.~ 0 °0 44..00 H. 4: .14.-I

, ..

..

.

.

.... . .

...

..

.. .

.

.

.

.

.

.

.

.. .

.

.

(4-04 lOin 0. (0.

.

.

+

I I

.

‘‘

.

44 o

. .

.

+ I I I 1 (4*1*4 I * ‘0 *44 0404(14 00. I 0.0)04010 Ol*401,-.~,. 04 1(404. 044+1. 11 14201

4o ‘I-

.

4+. (4: I..0 (.41(4

I +

.

‘

:‘::~..‘‘

(

.~ ..

‘

...

. .

.

..

.

.‘ .

:

.

:

I I I I + I I I

444.4410 010101004 -.44 ‘0 * I 01010104040. I 0......-I’4*4I 0 .-I1JI 154 — lEo 1.4154

)-

11-0. 1-IO~ —150 040104040)0 I 0 (4-

000

o

44

44

E 04.4 ~0 0.1~4: 00. 0.

0. .404 ‘0 1-01 21 o