% Latex file % New efficient programs to calculate general recoupling coefficients % Part II: Evaluation of a summation formula % V. Fack, S.N. Pitre, J.Van der Jeugt % Comp. Phys. Commun. 86 (1995), 105-122. % \documentclass[12pt]{article} %\usepackage{a4} %\setlength{\textwidth}{16cm} %\setlength{\textheight}{22cm} %\setlength{\evensidemargin}{0cm} %\setlength{\oddsidemargin}{0cm} \newcommand{\sixj}[6]{\left\{ \begin{array}{ccc} #1 & #2 & #3 \\ #4 & #5 & #6 \end{array} \right\}} \newcommand{\ninej}[9]{\left\{ \begin{array}{ccc} #1 & #2 & #3 \\ #4 & #5 & #6 \\ #7 & #8 & #9 \end{array} \right\}} \newcommand{\recoup}[2]{\langle \: #1 \: | \: #2 \: \rangle} \newcommand{\bra}[1]{\langle \: #1} \newcommand{\ket}[1]{| \: #1 \: \rangle} \newcommand{\kmin}{k^{\rm (min)}} \newcommand{\kmax}{k^{\rm (max)}} \begin{document} \addtolength{\baselineskip}{0.2\baselineskip} % % TITLE % ----- % \title{New efficient programs to calculate general recoupling coefficients \\ Part II~: Evaluation of a summation formula} \author{V.~Fack, S.N.~Pitre, J.~Van der Jeugt\thanks{Senior Research Associate N.F.W.O.~(National Fund for Scientific Research, Belgium)} \\ \small Department of Applied Mathematics and Computer Science, \\ \small Universiteit Gent, Krijgslaan 281--S9, B--9000 Gent, Belgium \\ \small E-mail~: Veerle.Fack@rug.ac.be} \date{Published in~: Comput.\ Phys.\ Commun.\ {\bf 86} (1995) 105--122} \maketitle \vfill \begin{abstract} In~\cite{FPV} we described a program {\tt NJFORMULA} to generate a summation formula in terms of 6-$j$ coefficients for a general angular momentum recoupling coefficient. Here we present programs for the second part of the problem, i.e.\ the numerical computation of a general recoupling coefficient by evaluating the formula generated by the program {\tt NJFORMULA}. We describe a sequential version {\tt NJSUMMATION}, written in C, which runs on a range of computers, as well as a parallel version {\tt NJSUMPAR}, also written in (Parallel) C, which runs on an arbitrary network of transputers. \end{abstract} \vfill % % PROGRAM SUMMARY % --------------- % \pagebreak {\large \bf PROGRAM SUMMARY} \begin{description} \item[Title of program~:] {\tt NJSUMMATION}, {\tt NJSUMPAR} \item[Catalogue number~:] \item[Program obtainable from~:] CPC Program Library, Queen's University of Belfast, \\ N.~Ireland (see application form in this issue). \item[Computers used (Operating systems)~:] 386/486-based PCs (MS-DOS, Linux~\cite{Linux}), Sun Sparc (Unix), T800 transputer system (MS-DOS). \item[Programming language used~:] C. \\ (Compilers~: Turbo C++~\cite{TC}, GNU CC~\cite{GnuC}, SPARCompiler C, 3L Parallel C~\cite{ParC}) \item[No.\ of lines in sequential program (module + 2 example programs)~:] 450 \item[No.\ of lines in parallel program~:] 500 \item[Keywords~:] atomic structure, nuclear structure, scattering, general recoupling coefficient, angular momentum, Racah coefficient, 3$n$-$j$ coefficient, recursive techniques, arbitrary nested loops, parallel processing, transputer, farming technique. \item[Nature of physical problem~:] A general recoupling coefficient for an arbitrary number of (integer or half-integer) angular momenta can be expressed as a formula consisting of products of 6-$j$ coefficients summed over a certain number of variables. Such a formula can be generated using the program {\tt NJFORMULA}~\cite{FPV}. The present programs perform the evaluation of a generated formula for given data of angular momenta. \item[Method of solution~:] A summation formula for a general recoupling coefficient can contain an arbitrary number of summation variables, which is determined by the binary trees involved in the search process when generating the formula~\cite{FPV}. To develop a general algorithm for the evaluation of such a formula, we use recursive programming techniques which allow to implement arbitrary nested loops. Apart from the sequential implementation {\tt NJSUMMATION} of this algorithm, we also present a parallel version {\tt NJSUMPAR}, which can be executed on a network of transputers. The latter uses a farming technique, where a master processor distributes work packets to several worker processors who calculate products of 6-$j$ coefficients. \item[Typical running time~:] For the sequential program the running time is determined mainly by the number of 6-$j$ coefficients computed in the evaluation process. Typical examples, involving upto ten thousand 6-$j$ coefficients, take less than 0.25 seconds on the Sun Sparc, less than 0.5 seconds on the 486-based PC, and less than 50 seconds on a single T800 transputer. More complicated examples, involving upto one million 6-$j$ coefficients, take less than 20 seconds on the Sun Sparc and less than 40 seconds on the 486-based PC. For the parallel program we focus our attention on the advantage factor of the parallel version compared with the sequential version. On a network of 4 transputers, the parallel program is 3 to 4 times faster than the corresponding sequential version. \end{description} % % INTRODUCTION % ------------- % \pagebreak {\large \bf LONG WRITE-UP} \bigskip \section{Introduction} % The need for angular momentum coupling coefficients arises in atomic and nuclear structure or scattering calculations. The problem of the calculation of a general recoupling coefficient for an arbitrary number of integer and half integer angular momenta has been studied by several authors~\cite{Bur,Bar,Lim} and is usually treated in two parts. In~\cite{FPV} we presented a program {\tt NJFORMULA}, which is a new approach to the first part of the problem, i.e.\ expressing the recoupling coefficient as a (multiple) sum over products of Racah coefficients multiplied by phase factors and square root factors. The algorithm described in~\cite{FPV} makes use of the method of binary tree transformations, which was first introduced by Burke~\cite{Bur} in the program {\tt NJSYM}. Another approach is followed by Bar-Shalom and Klapisch~\cite{Bar} (and later also by Lima~\cite{Lim}), who calculate angular momentum recoupling coefficients by means of more advanced graphical methods as developed by Yutsis, Levinson and Vanagas~\cite{Yut} and explained in many books~\cite{Bie,Bri,Jud,Elb}. This has given rise to a fast program {\tt NJGRAF}~\cite{Bar}, yielding efficient summation formulae (i.e.\ with a minimal number of summation variables) in most cases. However in~\cite{FPV} we showed that the simpler method of binary tree transformations can also be used to develop a fast program generating efficient formulae. Both programs {\tt NJSYM} and {\tt NJGRAF} also provide a subroutine {\tt GENSUM}, for the second part of the problem, i.e.\ the numerical evaluation of a general recoupling coefficient for which a summation formula has been generated. The program {\tt NJGRAF} implements a modified version of {\tt GENSUM}, which is more efficient than the version in {\tt NJSYM}. As remarked by Bar-Shalom and Klapisch, in the original {\tt GENSUM} of Burke, the limits of the values of the summation variables are not calculated accurately, but tests are performed on all the arguments of the 6-$j$ coefficients in order to check whether the latter should be computed or not. In their version of {\tt GENSUM}, Bar-Shalom and Klapisch build a matrix expressing the connections between variables and perform a sort of diagonalisation in order to find the most stringent limitations to each variable. These modifications are responsible for a noticeable part of the total gain of time for {\tt NJGRAF} over {\tt NJSYM}, which is upto 2 orders of magnitude. In the present paper we describe our approach to this second part of the problem. The procedure {\tt evaluate\_formula} is the equivalent of {\tt GENSUM} in the programs {\tt NJSYM} and {\tt NJGRAF}, and uses a completely different method viz recursion to calculate the sum over several variables. A general algorithm for the evaluation of a summation formula with an arbitrary number of summation variables is developed, in which nested loops of arbitrary depth are `simulated' by means of recursive programming techniques. The fact that each step in this recursive algorithm basically handles only one summation variable, also allows to formulate a clear algorithm for the determination of the range of a summation variable. We also discuss the possibilities for parallelisation of the numerical evaluation of a summation formula for a general recoupling coefficient. The use of parallel computers for numerical computation is of recent origin and it is conceivable that parallel systems will become easily accessible for scientists in various areas, thereby motivating this aspect of our program. The parallel implementation discussed here runs on a network of transputers, with an arbitrary configuration. A farming technique is used, with a master processor distributing work (consisting of the calculation of products of 6-$j$ coefficients) to worker processors. It is a rather straightforward extension of a parallel program we developed for the computation of 9-$j$ coefficients~\cite{FVR}. The structure of the present paper is as follows. In section~\ref{sec-meth} we discuss how a recursive algorithm can be built for the evaluation of a general summation formula, as well as how we determine the range of a summation variable using the triangle conditions. Section~\ref{sec-pro} describes the (sequential) program {\tt NJSUMMATION} and some implementation details concerning the datastructures used and the evaluation algorithm. The parallel implementation {\tt NJSUMPAR} of the program is described in section~\ref{sec-par}, where we show how the task can be broken up into subtasks, using the farming technique, and discuss some specific problems arising in this case. Finally in section~\ref{sec-dis} some results obtained by our programs are discussed and compared with results obtained by the program {\tt NJGRAF}~\cite{Bar}. % % METHOD OF SOLUTION % ------------------ % \section{Method of solution} \label{sec-meth} % % EVALUATION OF GENERAL SUMMATION FORMULA % \subsection{Evaluation of a general summation formula} The evaluation of an $n$-fold summation of the form \begin{eqnarray} S & = & \sum_{k_1} \sum_{k_2} \ldots \sum_{k_n} P_{k_1,\ldots,k_n}, \label{eq-sum} \end{eqnarray} where $P_{k_1,\ldots,k_n}$ stands for a product of 6-$j$ coefficients (including phase factors and square root factors) in terms of the summation variables $k_1,\ldots,k_n$ (and of angular momenta $j_i$), consists basically of an algorithm with $n$ nested for-loops~: {\sl \begin{tabbing} ww\=ww\=ww\=ww\=ww\=ww\= \kill \\ \> determine $\kmin_1$ and $\kmax_1$ \\ \> for $k_1$ from $\kmin_1$ to $\kmax_1$ do \\ \> \> determine $\kmin_2$ and $\kmax_2$ \\ \> \> for $k_2$ from $\kmin_2$ to $\kmax_2$ do \\ \> \> \> \ldots \\ \> \> \> \> determine $\kmin_n$ and $\kmax_n$ \\ \> \> \> \> for $k_n$ from $\kmin_n$ to $\kmax_n$ do \\ \> \> \> \> \> compute the next term $P_{k_1,\ldots,k_n}$ of the summation \\ \> \> \> \> \> add this to the temporary result for $S$ \end{tabbing}} For a constant value of $n$ it is easy to write a procedure evaluating an $n$-fold summation. However it would be unwise to write a separate evaluation procedure for every possible value of $n$, all of which will be different but also very alike. The obvious means is to write only one general evaluation procedure, that takes a variable $n$ as a parameter, for which the actual value will be given at runtime. This can be done by making use of recursive programming techniques, which allow to `simulate' nested for-loops of arbitrary depth~\cite{Rohl}. To that purpose, note that the actions in every for-loop in the above algorithm are very similar~: for every summation variable $k_i$ we have to determine the range $\kmin_i \to \kmax_i$ and then perform a for-loop for $k_i$ over this range, in which each step requires similar actions for $k_{i+1}$. This idea is expressed in the following recursive algorithm~{\tt UNOPT}, which describes how to evaluate a summation $S_i$ over the variables $k_i,\ldots,k_n$, with the current values for $k_1,\ldots,k_{i-1}$, for $1\leq i\leq n$~: {\sl \begin{tabbing} ww\=ww\=ww\=ww\=ww\=ww\=ww\= \kill \\ \> determine $\kmin_i$ and $\kmax_i$ \\ \> for $k_i$ from $\kmin_i$ to $\kmax_i$ do \\ \> \> evaluate $S_{i+1}$ over $k_{i+1},\ldots,k_n$, with current $k_1,\ldots, k_{i-1}$ and $k_i$ (recursively) \\ \> \> add this $S_{i+1}$ (for current $k_i$) to the temporary result for $S_i$ \end{tabbing}} The terminating case of the recursion is reached when $i=n+1$; at that point the actual computation of the product $P_{k_1,\ldots,k_n}$ of the 6-$j$ coefficients for the current values of the summation variables is performed. The complete sum $S$ is obtained by calling this algorithm to evaluate the summation $S_1$ over the variables $k_1,\ldots,k_n$. Using the Horner scheme of polynomial evaluation~\cite{Rao81}, the above algorithm can be optimised by `pulling out' independent factors for every summation. Denoting by $P_{k_i}$ the product of 6-$j$ coefficients (and the corresponding phase factors and square root factors) that contain only the summation variable $k_i$ (together with some of the summation variables $k_1,\ldots,k_{i-1}$ and some of the angular momenta), then the summation~(\ref{eq-sum}) can be written as~: \begin{eqnarray} \sum_{k_1} \sum_{k_2} \ldots \sum_{k_n} P_{k_1,\ldots,k_n} & = & \sum_{k_1} P_{k_1} \sum_{k_2} P_{k_2} \ldots \sum_{k_n} P_{k_n}. \end{eqnarray} Taking this into account, we obtain the following algorithm {\tt OPT}, to evaluate a summation $S_i$ over the variables $k_i,\ldots,k_n$, with current values for $k_1,\ldots,k_{i-1}$, for $1\leq i\leq n$~: {\sl \begin{tabbing} ww\=ww\=ww\=ww\=ww\=ww\=ww\= \kill \\ \> determine $\kmin_i$ and $\kmax_i$ \\ \> for $k_i$ from $\kmin_i$ to $\kmax_i$ do \\ \> \> compute $P_{k_i}$ \\ \> \> evaluate $S_{i+1}$ over $k_{i+1},\ldots,k_n$, for current $k_1,\ldots,k_{i-1}$ and $k_i$ (recursively) \\ \> \> add $P_{k_i} \times S_{i+1}$ (for current $k_i$) to the temporary result for $S_i$. \end{tabbing}} The terminating case of the recursion, when $i=n+1$, needs no computations at all and simply returns 1 as the result for $S_{n+1}$. % % RANGE OF A SUMMATION VARIABLE % \subsection{Determining the range of a summation variable $k_i$} As described by Burke~\cite{Bur} and Bar-Shalom and Klapisch~\cite{Bar}, the range of values taken by the summation variables can be decided by examining the triangular conditions placed on the different values in the 6-$j$ coefficients appearing in the product of the summation. Since the recursive algorithm, described in the previous section, basically considers only one summation variable $k_i$ in each algorithm step, we can develop a straightforward algorithm for the determination of the range of this variable. A 6-$j$ coefficient of the following general form \begin{equation} \sixj{j_1}{j_2}{j_3}{l_1}{l_2}{l_3} \end{equation} is trivially zero unless the following four triangular conditions~\cite{Edm} are satisfied~: \begin{equation} \triangle(j_1 j_2 j_3), \triangle(j_1 l_2 l_3), \triangle(l_1 j_2 l_3), \triangle(l_1 l_2 j_3). \end{equation} A triangular condition $\triangle(t_1 t_2 t_3)$ means that \begin{eqnarray} \triangle(t_1 t_2 t_3) & = & \left[ (t_1+t_2-t_3)!(t_1-t_2+t_3)!(-t_1+t_2+t_3)! \over{(t_1+t_2+t_3+1)!} \right]^{1\over2} \end{eqnarray} should be non-zero or the following inequality should be satisfied~: \begin{eqnarray} | t_1-t_3 | \leq t_2 \leq t_1+t_3. \end{eqnarray} For each 6-$j$ coefficient appearing in the products of the summation, these triangular inequalities should be checked before actually computing the 6-$j$ coefficient, to avoid unnecessary computations. However, some of the 6-$j$ coefficients will also involve one or more summation variables. In that case the triangular inequalities can be used to determine the range of the summation variables, in the following way. At each step $i$ in the recursive algorithm, described above, we are left with determining the range of the summation variable $k_i$. At this point we know, from the way in which the summation is evaluated, that the preceding summation variables $k_1,\ldots,k_{i-1}$ have already been given a current numerical value. To find $\kmin_i$ and $\kmax_i$ we have to consider all possible triangular inequalities appearing in the 6-$j$ coefficients involving $k_i$, which will be of the form \[ | t - t' | \leq k_i \leq t + t', \] where $t$ and $t'$ stand for some $j_x$ coefficient or $k_x$ variable. Conditions involving also some of the summation variables $k_{i+1},\ldots,k_n$ as a value for $t$ or $t'$, can be ignored at this point, since they cannot give any more information about the range of $k_i$. But for conditions involving only $j_x$ angular momentum values or $k_1,\ldots,k_{i-1}$ summation variables, the values of $t$ and $t'$ will be specific numerical values, which does give some information about the range of $k_i$. In particular, the value of $\kmin_i$ will be the maximum of all appropriate values of $ | t - t' | $, while the value of $\kmax_i$ will be the minimum of all appropriate values of $ t + t' $. % \subsubsection*{Example} As an example we show how to obtain the range of the summation variables $k_1$ and $k_2$ in the following summation formula~: \begin{eqnarray} & & \sum_{k_1} \sum_{k_2} (-1)^{k_1+k_2+j_2+j_8} (2k_1+1) (2k_2+1) \sixj{j_1}{j_2}{j_6}{j_8}{j_9}{k_1} \sixj{j_7}{j_3}{j_8}{j_4}{k_1}{j_2} \nonumber \\ & & \mbox{} \qquad \qquad \times \sixj{j_7}{j_4}{k_1}{k_2}{j_1}{j_9} \sixj{j_5}{j_4}{j_7}{k_2}{j_9}{j_3} \sixj{j_7}{j_1}{k_2}{j_4}{j_2}{j_8}. \label{eq-form} \end{eqnarray} In the first recursive algorithm step we have to determine the range of the summation variable $k_1$. For the first 6-$j$ coefficient in the formula the following four triangular conditions should be satisfied~: \[ \triangle(j_1 j_2 j_6), \triangle(j_1 j_9 k_1), \triangle(j_8 j_2 k_1), \triangle(j_8 j_9 j_6), \] only two of which involve $k_1$. So we already have two triangular inequalities which help to decide the range of $k_1$~: \[ |j_1-j_9| \leq k_1 \leq j_1+j_9 \qquad {\rm and} \qquad |j_8-j_2| \leq k_1 \leq j_8+j_2. \] In a similar way we obtain two more triangular inequalities from the second 6-$j$ coefficient in the formula~: \[ |j_7-j_2| \leq k_1 \leq j_7+j_2 \qquad {\rm and} \qquad |j_4-j_8| \leq k_1 \leq j_4+j_8. \] The third 6-$j$ coefficient in the formula gives rise to the following triangular inequalities~: \[ \triangle(j_7 j_4 k_1), \triangle(j_7 j_1 j_9), \triangle(k_2 j_4 j_9), \triangle(k_2 j_1 k_1), \] of which we only have to consider those involving $k_1$ and can ignore those involving $k_2$ at this point. This gives one more condition to decide on the range of $k_1$~: \[ |j_7-j_4| \leq k_1 \leq j_7+j_4. \] The other two 6-$j$ coefficients in the formula involve only $k_2$ and will give no further information in determining the range of $k_1$. So finally the values of $\kmin_1$ and $\kmax_1$ are decided~: \begin{eqnarray*} \kmin_1 & = & \max (|j_1-j_9|,|j_8-j_2|,|j_7-j_2|,|j_4-j_8|,|j_7-j_4|), \nonumber \\ \kmax_1 & = & \min ( j_1+j_9 , j_8+j_2 , j_7+j_2 , j_4+j_8 , j_7+j_4 ). \end{eqnarray*} For every value of $k_1$ in this range we have to compute the inner summation over $k_2$, for which we also have to determine the range first. The triangular inequalities involved arise first of all from the third 6-$j$ coefficient in the formula~: \[ |j_4-j_9| \leq k_2 \leq j_4+j_9 \qquad {\rm and} \qquad |j_1-k_1| \leq k_2 \leq j_1+k_1. \] Note that $k_1$ has a specific numerical value at this point, and so it can be used to determine the range of $k_2$. Also from the last two 6-$j$ coefficients in the formula some triangular conditions can be obtained, finally giving rise to the following values for $\kmin_2$ and $\kmax_2$~: \begin{eqnarray*} \kmin_2 & = & \max\ (|j_4-j_9|,|j_1-k_1|,|j_4-j_3|,|j_9-j_7|,|j_7-j_1|,|j_4-j_2|), \nonumber \\ \kmax_2 & = & \min\ ( j_4+j_9 , j_1+k_1 , j_4+j_3 , j_9+j_7 , j_7+j_1 , j_4+j_2 ). \end{eqnarray*} % % CALCULATION OF 6-J COEFFICIENTS % \subsection{Calculation of 6-$j$ coefficients} After each variable of summation has been replaced by its current value, we have an expression consisting of a product of 6-$j$ coefficients which can be computed numerically. For the computation of a 6-$j$ coefficient we use the classical single sum expression due to Racah~\cite{Rac}~: \begin{eqnarray} \lefteqn{\sixj{j_1}{j_2}{j_3}{l_1}{l_2}{l_3} \quad = \quad (-1)^{j_1+j_2+l_1+l_2} \Delta(j_1j_2j_3)\Delta(l_1l_2j_3) \Delta(l_1j_2l_3)\Delta(j_1l_2l_3)} \nonumber \\ & & \times \sum_{k=\kmin}^{\kmax} \frac{(-1)^k (j_1+j_2+l_1+l_2+1-k)!}{k!(k-\alpha_1)!(k-\alpha_2)! (\beta_1-k)!(\beta_2-k)!(\beta_3-k)!(\beta_4-k)!}, \label{eq-6j} \end{eqnarray} with \[ \alpha_1=j_1+l_1-j_3-l_3\,, \qquad \alpha_2=j_2+l_2-j_3-l_3\,, \] \[ \beta_1=j_1+j_2-j_3\,, \qquad \beta_2=l_1+l_2-j_3\,, \qquad \beta_3=j_1+l_2-l_3\,, \qquad \beta_4=l_1+j_2-l_3\,, \] and \[ \kmin = \max(\alpha_1,\alpha_2,0)\,, \qquad \kmax = \min(\beta_1,\beta_2,\beta_3,\beta_4)\,. \] When implementing the computation of 6-$j$ coefficients, a look-up table for the factorials is provided, as suggested by Scott and Hibbert~\cite{SH}. This is useful since factorials are often needed more than once in the computations. As described by Rao {\sl et al.}~\cite{RRC}, the logarithms of the factorials are stored instead of the actual values themselves, for the sake of accuracy. % % PROGRAM DESCRIPTION % ------------------- % \section{Program description} \label{sec-pro} The program {\tt NJSUMMATION} presented here, is a module to be used together with the program {\tt NJFORMULA}, described in an earlier paper~\cite{FPV}. While the program {\tt NJFORMULA} provides routines for the {\em generation} of a summation formula for a general recoupling coefficient, the present program {\tt NJSUMMATION} implements a function {\tt evaluate\_formula} for the {\em evaluation} of such a summation formula for given values of the angular momenta. The module also contains some auxiliary routines, such as a procedure {\tt read\_jvals} to read actual values for the angular momenta, a function {\tt trian} to check triangle conditions imposed on the angular momenta, a procedure {\tt order} to order the 6-$j$ coefficients in a formula w.r.t.\ the summation variables, and a procedure {\tt comp\_fac} to build a look-up table for the factorials needed in the computation of 6-$j$ coefficients, for which a routine {\tt comp\_6j} is supplied. In this section we describe the implementation as well as some special features concerning the functionality of these routines. We also show how to use them in a simple main program and how to combine the modules {\tt NJFORMULA} and {\tt NJSUMMATION}. % % DATASTRUCTURES % \subsection{Datastructures} The datastructure used in this program to represent a summation formula is essentially the same as in the program {\tt NJFORMULA}. Again the index $i$ of a $j_i$ coefficient is used to identify the $j_i$ coefficient, while a negative index $-i$ is used to identify a summation variable $k_i$ of the summation formula. For the actual values of the angular momenta, which are integer or half-integer values, the double value is stored, which allows to work with integers. (We would like to remark that the procedure {\tt read\_jvals} expects a value $2j_i+1$ to be given as input, this is done for reasons of compatibility with the program {\tt NJGRAF}~\cite{Bar}.) A summation formula is represented by a record structure {\tt FORMULA}, which contains first of all the same fields as the structure described in~\cite{FPV}~: The fields {\tt nrjs}, {\tt nrks}, {\tt nrsixjs} contain resp.\ the number of $j_i$ coefficients in the recoupling coefficient, the number of summation variables $k_i$ and the number of 6-$j$ coefficients in the products of the formula. Fields {\tt jsigns}, resp.\ {\tt ksigns}, and {\tt jsqrts}, resp.\ {\tt ksqrts}, are arrays with an entry for every $j_i$ coefficient, resp.\ $k_i$ summation variable, such that \begin{quote} {\tt jsigns[i]}, resp.\ {\tt ksigns[i]} = power of $(-1)^{j_i}$, resp.\ $(-1)^{k_i}$, in the formula, \\ {\tt jsqrts[i]}, resp.\ {\tt ksqrts[i]} = power of $\sqrt{(2 j_i + 1)}$, resp.\ $\sqrt{(2 k_i + 1)}$. \end{quote} The field {\tt sixjs} is an array for which every entry is a representation of one of the 6-$j$ coefficients of the formula, as an array with 6 elements which are $j_i$ coefficients or $k_i$ summation variables, represented by an index $i$ or $-i$. For the purposes of the present program some extra fields have been added to this structure {\tt FORMULA}~: A field {\tt js} of the array type will be used to store the actual values for the $j_i$ angular momenta during the evaluation of the formula. The field {\tt asixjs} is a copy of the array {\tt sixjs} except that it contains the actual $j_i$ and $k_i$ values wherever possible; e.g.\ if $j_1=2$, $j_6=2$, $j_8=4$, $j_9=2$, $k_1=4$ and {\tt sixjs[1]}=$(1,-1,6,8,9,-2)$, then {\tt asixjs[1]}=$(2,4,2,4,2,-2)$. In the following we will assume that the 6-$j$ coefficients in the formula are ordered in increasing order w.r.t.\ the summation variables, i.e.\ first the 6-$j$ coefficients with no summation variables, then the 6-$j$ coefficients with only $k_1$, then the 6-$j$ coefficients with only $k_1$ and $k_2$, and so on. This ordering will be necessary for implementing the optimised algorithm~{\tt OPT}, which `pulls out' 6-$j$ coefficients independent of the summation variables $k_{i+1},\ldots,k_n$ at each step $i$ in the algorithm. For that purpose an additional field {\tt kp} of the array type is added to the datastructure {\tt FORMULA}, to keep track of how many 6-$j$ coefficients can be pulled out at different levels of the evaluation. An entry {\tt kp[i]} will contain the number of 6-$j$ coefficients that contain only the summation variables $k_1,\ldots,k_i$, where the entry {\tt kp[0]} signifies the number of 6-$j$ coefficients containing no summation variables; e.g.\ for the summation~(\ref{eq-form}), which is already ordered according to our conventions, we will have {\tt kp[0]} = 0, {\tt kp[1]} = 2 and {\tt kp[2]} = 5. % % EVALUATING A SUMMATION FORMULA % \subsection{Evaluating a summation formula} The function {\tt evaluate\_formula} performs the evaluation of a given summation formula, assuming that the structure of the formula as well as the actual $j_i$ angular momentum values have been filled in properly. It returns a real (double precision) value which results from this evaluation. Before starting the actual evaluation, it first performs some initial computations and checks. The triangle conditions are checked for all the 6-$j$ coefficients for the triads that contain only angular momenta (by means of the function {\tt trian}); if these are not satisfied the recoupling coefficient is trivially zero. Next it is checked if the product of 6-$j$ coefficients is ordered as described above and if the {\tt kp} array has already been set up properly; if necessary this ordering and building of {\tt kp} is done here, by means of a procedure {\tt order}. After that the array {\tt asixjs} is initialised, as a copy of the array {\tt sixjs}, replacing all positive indices $i$ by the corresponding $j_i$ value, while keeping the negative indices as they are (they correspond to summation variables $k_i$ and will be replaced later during the recursive evaluation process). Then the product of the 6-$j$ coefficients containing only $j_i$ angular momentum values, together with the corresponding square root factors, is computed. For the calculation of a 6-$j$ coefficient we have written a function {\tt comp\_6j}, which implements formula~(\ref{eq-6j}) as described above (we assume a look-up table for the factorials has been built initially by the program, by means of the procedure {\tt comp\_fac}). Finally, if the formula does not contain a summation (i.e.\ if \mbox{{\tt nrks} = 0}), this value (multiplied with the appropriate sign factor) is returned as the result of the function, otherwise the recursive evaluation of the summation is started by calling the function {\tt r\_eval\_formula}. The recursive function {\tt r\_eval\_formula} takes as parameters the formula to be evaluated and a value indicating the step $i$ in the algorithm, which is the index $i$ indicating for which variable $k_i$ the summation will be evaluated next. An additional parameter is a sign factor, which is a integer power of $-1$, with which the obtained sum has to be multiplied. Initially the function is called for step $i=1$, with the sign factor being the powers of $-1$ containing only angular momenta. The result obtained from this initial call, multiplied with the product of the square root factors and the product of 6-$j$ coefficients containing only angular momenta (which have been computed during the above initialisation), gives the final result for the recoupling coefficient. In each call to {\tt r\_eval\_formula} the range of a summation variable $k_i$ is determined as described previously. Different ranges are obtained for the same summation variable depending on its position in different 6-$j$ coefficients, and the final range is decided by keeping the maximum of all the lower limits and the minimum of all the upper limits. Next a loop is performed from this $\kmin_i$ to $\kmax_i$, which replaces every occurrence of the index $-i$ in the 6-$j$ coefficients with the current value of $k_i$. Then the product of 6-$j$ coefficients containing $k_i$ (but no variables $k_{i+1},\ldots,k_n$) can be computed, which means we `pull out' those 6-$j$ coefficients which are independent of those further summation variables. It is at this stage that the information in the array {\tt kp} is used. Finally the rest of the summation, i.e.\ over the summation variables $k_{i+1},\ldots,k_n$, for the current value of $k_i$, must be evaluated. This is done by a recursive call to the function {\tt r\_eval\_formula}, for step $i+1$, with the new sign factor being the old sign factor plus the current sign factor for $k_i$. This result is multiplied with the current square root factor for $k_i$ and the product of 6-$j$ coefficients just computed, and this is added to a temporary value for the current summation. This procedure is continued recursively till all the summation variables have been replaced by their values, i.e.\ until a call of {\tt r\_eval\_formula} for step $n+1$. The terminating case of the recursion needs no computations and simply returns the appropriate sign $-1$ or $+1$, depending on the sign factor. % % MAIN PROGRAM % \subsection{Main program} The following simple main program can be used to calculate a general recoupling coefficient and to evaluate it for different sets of actual angular momentum values. The input required here is compatible with the input style of the program {\tt NJGRAF}~\cite{Bar}. A header file {\tt njsummat.h} must be included, since it contains a set of declarations which are needed by the program. The program starts by generating a summation formula for a general recoupling coefficient, using the procedure {\tt read\_njsym} to read a recoupling coefficient, the function {\tt generate\_formula} to generate the formula and the procedure {\tt write\_formula} to write the obtained formula. The functionality of these routines is described in detail in~\cite{FPV}. Then the program builds a look-up table for the factorials, using the procedure {\tt comp\_fac}. Next it performs a loop, in each step reading a set of actual angular momentum values (using the procedure {\tt read\_jvals}), evaluating the obtained formula for this set of data (using the function {\tt evaluate\_formula}), and finally writing the obtained value for the recoupling coefficient. We also keep track of the number of 6-$j$ coefficients actually computed during the evaluation of a recoupling coefficient, for which purpose a global variable {\tt count} is updated automatically. As mentioned earlier, the procedure {\tt read\_jvals} expects a value $2j_i+1$ for an angular momentum $j_i$. Moreover, after reading the values for the angular momenta, it also reads a line with the same number of {\tt T} or {\tt F} characters as the number of angular momenta, but these characters are not stored. All this is done for reasons of compatibility of the input with the program {\tt NJGRAF}. \begin{verbatim} #include "njsummat.h" main () { int nrjs, next; NODE *bra, *ket; DELTAS delta; FORMULA *form; double recup; read_njsym (&nrjs, &bra, &ket, &delta); form = generate_formula (nrjs, bra, ket); write_formula (form); comp_fac (); do { read_jvals (form); recup = evaluate_formula (form); printf ("recup = %e\n", recup); printf ("number of 6-j's computed = %d\n", count); scanf ("%d", &next); } while (next != 0); } \end{verbatim} An example of an input file for this program is the following~: \begin{verbatim} 18 6 1 1 2 8 3 4 9 8 9 10 6 7 11 5 11 12 10 12 13 2 4 14 6 14 15 1 7 16 3 16 17 15 17 18 5 18 13 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 T T T T T T T T T T T T T T T T T T 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 T T T T T T T T T T T T T T T T T T 1 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 T T T T T T T T T T T T T T T T T T 0 \end{verbatim} % % INDEPENDENT SUMS % \subsection{The case of several independent sums} For some recoupling coefficients it happens that the summation variables in the formula can be partitioned into disjoint sets, such that the summation breaks up into several products of independent summations. Such independent factors can, and should, be evaluated separately, since this will result in a smaller number of 6-$j$ coefficients actually being computed. Bar-Shalom and Klapisch~\cite{Bar} handle this case by separating the obtained formula into products of sums, in between the generation and the evaluation process. So far we have not taken this possibility into account when developing our algorithms, which means that the case of several independent sums is not yet handled in the most efficient way. However, the essential ideas of our algorithms still hold and we will show here how our routines can be used to handle this case appropriately. It is easily seen that, in our approach, such independent sums arise only when for the recoupling coefficient (given as a two coupling trees~\cite{FPV}) there exist some pairs of angular momenta $(j_{i_1},j_{i_2})$ such that the coupled node $j_{i_1}$ in the bra tree and the coupled node $j_{i_2}$ in the ket tree are essentially the same, in the sense that they have the same set of leaf nodes. Indeed, in that case the transformation of the subtree $j_{i_1}$ in the bra tree to the subtree $j_{i_2} \equiv j_{i_1}$ in the ket tree will be independent of the rest of transformations needed in the coupling trees, and the subformula obtained will be an independent factor in the final summation formula. So from the beginning, the original problem, i.e.\ a recoupling coefficient resulting in a summation formula with several independent sums, can be split into several subproblems, i.e.\ recoupling coefficients for which the summation formula contains only one independent sum. For each of these subproblems the routines we have supplied so far can be used efficiently. The only routine which we have to add here is a procedure {\tt split}, which splits the two coupling trees into a set of pairs of coupling trees. For this purpose it uses the information returned by the input routines in~\cite{FPV}, which is in the form of a list of delta functions $\delta(j_{i_1},j_{i_2})$, indicating the pairs of equivalent angular momenta $j_{i_1}$ and $j_{i_2}$. The procedure {\tt split} returns an array with the independent bra subtrees as well as an array with the independent ket subtrees, obtained from a given bra and ket tree, which have to be supplied at call-time in the zeroth element of the corresponding arrays. Also the number of independent pairs of subtrees (apart from the main pair of trees) obtained is returned. The above main program can easily be extended in such a way that, after reading a recoupling coefficient, it calls the procedure {\tt split} when necessary. Next it has to generate each independent summation subformula separately, and finally it has to evaluate the product of these subformulae for a range of given $j_i$ values. To distribute the actual $j_i$ values over all the independent subformulae, we supply a routine {\tt copy\_jvals}. Below we give the listing of this main program. \begin{verbatim} #include "njsummat.h" main () { int i, nrjs, nrtrees, next; NODE *bra[10], *ket[10]; DELTAS delta; FORMULA *form[10]; double recup; read_njsym (&nrjs, &bra[0], &ket[0], &delta); if (delta.n > 0) split (delta, &nrtrees, bra, ket); for (i=0; i<=nrtrees; i++) { form[i] = generate_formula (nrjs, bra[i], ket[i]); write_formula (form[i]); } comp_fac (); do { read_jvals (form[0]); recup = evaluate_formula (form[0]); for (i=1; i<=nrtrees; i++) { copy_jvals (form[0], form[i]); recup *= evaluate_formula (form[i]); } printf ("recup = %e\n", recup); printf ("number of 6-j's computed = %d\n", count); scanf ("%d", &next); } while (next != 0); } \end{verbatim} % % THE PARALLEL VERSION % -------------------- % \section{The parallel version} \label{sec-par} Apart from the sequential program {\tt NJSUMMATION} described above, we also present a parallel program {\tt NJSUMPAR} to calculate a general recoupling coefficient. This parallel version is developed on a network of transputers, using a so-called farming technique, which makes the program completely independent of the actual transputer configuration. The method used here is similar to the method we used earlier when developing a parallel program to calculate 9-$j$ coefficients~\cite{FVR}. % % INTRODUCTION % \subsection{Transputers and parallel programming} For our purposes a transputer can be viewed as a processor that can run program code independently and that has four links, through which it can communicate with other processors in the network. Our specific system consists of 20 transputers and is hosted by four PC-AT microcomputers. For that purpose, four of the transputers have host interface capability; they are T800/20 MHz with 4Mb of RAM. The other 16 transputers are T800/20 MHz with 1Mb RAM. The PC-AT hosts run server software providing access to the PC-AT for system services such as file I/O, keyboard input and screen output. Other software (such as compilers, applications, etc.) runs on the transputer nodes. The parallel program {\tt NJSUMPAR} is written using the stand-alone Parallel C compiler from 3L~\cite{ParC}, which is basically a standard C compiler extended with several features supporting the concurrency offered by a transputer system. Parallel C uses the common abstract model of parallel processing in transputer systems, which is based on the idea of {\em communicating sequential processes}. A complete application is viewed as a collection of concurrently active sequential processes ({\em tasks\/}), which can communicate with each other over {\em channels}, each channel connecting one process to one other process. It is important to note that such channel communications are synchronised~: a process wishing to send a message over a channel is forced to wait until the receiving process reads the message. In that way, tasks can be treated as the atomic building blocks for parallel programs, being software `black boxes' connected together by channels. One possible type of application in Parallel C is a so-called {\em farming application}, which can run on an {\em arbitrary\/} network of transputers. A problem is suited for such an implementation when its solution consists in applying the same technique several times to different data. In the processor farm technique an application consists of one {\em master\/} task, organising the work, and any number of anonymous {\em worker\/} tasks, all performing similar parts of the work (thus all running the same code). Sending work packets by the master to the workers or returning results by the workers to the master, is done transparently by {\em routing software\/} which is supplied with the Parallel C compiler. Another important concept, supported by Parallel C, is that of a task being multi-threaded. This means that a task can contain any number of concurrent processes (called {\em threads\/}) running on the same processor, each of which is independently executing the code of the task. Each thread has its own private stack of data but shares the rest of its data with all the other threads in the same task. Threads within a task can communicate with each other via shared memory. The software construct of semaphores, which is commonly used in order to prevent threads from interfering with each other while operating on shared data, is available. Alternatively, internal channels can be used to synchronize the threads' operations and transmit data between them by passing messages. % % MASTER TASK % \subsection{The master task of {\tt NJSUMPAR}} The task {\tt MASTER} of the farming application {\tt NJSUMPAR} runs on the master transputer, which is a transputer with host interface capability. It is responsible for organising the work, which means first of all breaking up the problem into independent work packets which are to be processed by the workers. In this case a work packet will consist of the computation of one term in the summation, which means the computation of a product of 6-$j$ coefficients multiplied by a constant factor. Thus the master program has all the units of the sequential program {\tt NJSUMMATION} except for the routines involved with the actual computation of 6-$j$ coefficients. It contains the recursive subroutine {\tt evaluate\_formula} which does the job of finding the ranges of different summation variables and then performs nested for-loops over these ranges recursively as described before. The major difference here is that once we have a product of 6-$j$ coefficients with all the actual values filled in, its value has to be computed by a worker task. Another difference in the parallel version is the fact that we use the unoptimised algorithm~{\tt UNOPT} instead of the optimised algorithm~{\tt OPT}. This is done because pulling out independent products of 6-$j$ coefficients for every summation, introduces a sort of data-dependency between the recursion steps, which makes it quite non-trivial to parallelise the optimised algorithm. However, trying to do so will be a topic of our further research, since optimising the algorithm in this way yields a sizeable reduction in the computation time (as is shown in section~\ref{sec-dis}). The master task also takes care of sending the work packets to the workers, for which it uses the Parallel C procedure {\tt net\_send}. This will send a work packet to a free worker, i.e.\ a worker who is ready to receive a work packet. If no workers are free, the {\tt net\_send} procedure will `sleep' and the sending process will be blocked until a worker becomes free. Finally the master also has to collect the results from the workers, which have to be added in order to obtain the final result to the problem. For this, it uses the Parallel C procedure {\tt net\_receive}, which is also blocking, i.e.\ the process `sleeps' until one of the workers is ready to send a result. Because of the fact that the {\tt net\_send} and {\tt net\_receive} procedures are blocking, the master task must be a two-threaded task, in which the operations of sending work packets to the workers and receiving results from the workers have to be performed in parallel. Otherwise the program would end in a `deadlock' situation, when the master is waiting for a worker to become free, but no worker can return its result (and become free) because the master is not ready to receive this result. Thus there is a send thread, which takes care of building the work packets (in a recursive algorithm as described above) and of sending these work packets to the workers (which happens at the lowest level in the recursion). There is also a receive thread performing a loop, which receives the result packets from the workers and adds them to obtain the final recoupling coefficient. An additional problem here is the fact that initially it is not known how many work packets will be sent by the master and thus how many packets the receive thread must receive. To solve this problem we can use semaphores in the following way. One semaphore is used as a sort of counter indicating how many packets the send thread has sent~: it is updated by the send thread each time it sends a packet and it is consulted by the receive thread to check if there are unreceived packets. To indicate that all packets have been sent, another semaphore is set by the send thread, which is checked at the appropriate time by the receive thread. The receive thread can finish its loop when there are no unreceived packets and the send thread has sent all its packets. % % WORKER TASK % \subsection{The worker task of {\tt NJSUMPAR}} The task {\tt WORKER} runs on several slave transputers and performs most of the actual computation. It basically consists of a simple sequential loop in which it accepts work packets from the master, performs the processing for which it is programmed on the data of the received work packet and sends the result back to the master. Since the essential work done by the {\tt WORKER} tasks is the computation of products of 6-$j$ coefficients, each {\tt WORKER} task starts by building a look-up table for the factorials, which contains the logarithms of factorials (as described earlier). Next it performs a `loop forever', computing products of 6-$j$ coefficients, in the following way. Using the Parallel C procedure {\tt net\_receive}, it reads from `the network' a one-record message from the master, which contains a constant factor and a list of actual 6-$j$ coefficients that have to be computed. It then computes the product of the received 6-$j$ coefficients, using the procedure {\tt comp\_6j} as described above. Finally it sends the computed result, multiplied by the given constant factor, back over the network to the master task, using the Parallel C procedure {\tt net\_send}. % % RESULTS AND DISCUSSION % ---------------------- % \section{Results and Discussion} \label{sec-dis} The program {\tt NJSUMMATION} is written in C and can be run on any system providing a C (or C++) compiler. We have tested the program on a range of machines (such as 386/486-based PCs running MS-DOS and Linux~\cite{Linux}, a Sun Sparc running Unix, and a single T800 transputer hosted by a PC running MS-DOS) using several compilers (such as Turbo~C++~\cite{TC}, Gnu~CC~\cite{GnuC}, SPARCompiler C, and 3L Parallel~C~\cite{ParC}). The parallel version {\tt NJSUMPAR} is written in Parallel C and runs on an arbitrary network of T800 transputers. To compare the results obtained by our program {\tt NJSUMMATION} with the results obtained by the program {\tt NJGRAF}~\cite{Bar}, we have compiled the program {\tt NJGRAF} on the Sun Sparc using the Sparc Fortran-77 compiler and on the transputer system using the 3L Parallel Fortran compiler~\cite{ParF}. The test cases for which we want to discuss our results, are a selection of those described in~\cite{FPV}, corresponding to the following recoupling coefficients~: {\small \begin{itemize} \item[($G_4$)] $\bra{((j_1,j_2)j_{12},(j_3,(j_4,(j_5,((j_6,j_7)j_{13},(j_8,(j_9,(j_{10},j_{11})j_{14})j_{15})j_{16})j_{17})j_{18})j_{19})j_{20})j_{21}}$ \\ $\ket{(j_5,(j_7,((j_6,(j_{11},(j_9,(j_{10},j_8)j_{22})j_{23})j_{24})j_{25},(j_1,(j_3,(j_2,j_4)j_{26})j_{27})j_{28})j_{29})j_{30})j_{21}}$ \item[($F_3$)] $\bra{(((j_1,j_2)j_7,(j_3,j_4)j_8)j_9,(j_5,j_6)j_{10})j_{11}}$ $\ket{((j_3,j_6)j_{12},((j_2,j_4)j_{13},(j_1,j_5)j_{14})j_{15})j_{11}}$ \item[($F_4$)] $\bra{(((j_1,j_2)j_7,(j_3,j_4)j_8)j_9,(j_5,j_6)j_{10})j_{11}}$ $\ket{((j_1,j_6)j_{12},((j_3,j_5)j_{13},(j_2,j_4)j_{14})j_{15})j_{11}}$ \item[($F_5$)] $\bra{(((j_1,j_2)j_8,(j_3,j_4)j_9)j_{10},((j_5,j_6)j_{11},j_7)j_{12})j_{13}}$ \\ $\ket{((j_1,j_7)j_{14},(((j_3,j_5)j_{15},(j_2,j_4)j_{16})j_{17},j_6)j_{18})j_{13}}$ \item[($F_6$)] $\bra{(((j_1,j_2)j_8,(j_3,j_4)j_9)j_{10},(j_5,(j_6,j_7)j_{11})j_{12})j_{13}}$ \\ $\ket{(j_5,((j_6,(j_2,j_4)j_{14})j_{15},(j_3,(j_1,j_7)j_{16})j_{17})j_{18})j_{13}}$ \item[($F_7$)] $\bra{((j_1,(j_2,j_3)j_8)j_9,((j_4,j_5)j_{10},(j_6,j_7)j_{11})j_{12})j_{13}}$ \\ $\ket{(((j_1,j_4)j_{14},j_6),j_{15},((j_5,j_2)j_{16},(j_7,j_3)j_{17})j_{18})j_{13}}$ \item[($F_8$)] $\bra{(((j_1,j_2)j_9,j_3)j_{10},(((j_4,j_5)j_{11},j_6)j_{12},(j_7,j_8)j_{13})j_{14})j_{15}}$ \\ $\ket{(((j_1,j_4)j_{16},j_7)j_{17},(((j_2,j_5)j_{18},(j_8,j_3)j_{19})j_{20},j_6)j_{21})j_{15}}$ \item[($F_9$)] $\bra{(((j_1,(j_2,j_3)j_{11})j_{12},((j_4,j_5)j_{13},j_6)j_{14})j_{15},(((j_7,j_8)j_{16},j_9)j_{17},j_{10})j_{18})j_{19}}$ \\ $\ket{(((j_2,j_4)j_{20},j_7)j_{21},((((j_1,j_8)j_{22},(j_9,j_5)j_{23})j_{24},j_{10})j_{25},(j_6,j_3)j_{26})j_{27})j_{19}}$ \end{itemize} } Of these recoupling coefficients, only ($G_4$) gives rise to a summation formula with several independent sums; this is a case with two independent sums. So in order to test the version of the program handling the case of several independent sums, we have constructed some test cases by combining a pair of simpler recoupling coefficients (as the ones given above) in a form of {\sl `direct product'}. This is obtained by coupling the bra (resp.\ ket) trees of both recoupling coefficients of the pair to form the bra (resp.\ ket) tree of the combined recoupling coefficient. The following combinations have been made here~: {\small \begin{itemize} \item[] ($F_{10}$) = ($F_3$) $\times$ ($F_4$) \item[] ($F_{11}$) = ($F_4$) $\times$ ($F_7$) \item[] ($F_{12}$) = ($F_7$) $\times$ ($F_8$) \end{itemize} } In table~1 we compare the results obtained by the sequential program {\tt NJSUMMATION} with the results of the program {\tt NJGRAF}. For each test case ($F_3$)--($F_9$) we have calculated the value of the recoupling coefficient for several sets $(2j_i+1)$ values, all of which are chosen to be equal. Both programs do indeed obtain the same value. Since running times (given in seconds, obtained on a Sun Sparc computer) are very small for most cases, we also compare the number of 6-$j$ coefficients actually computed during the evaluation of the recoupling coefficient. Since none of these examples contain several independent sums, the non-splitting version of {\tt NJSUMMATION} has been used to generate these results, but of course the splitting version yields exactly the same results. For those cases where both programs have obtained a similar formula, i.e.\ with the same number of summation variables and the same number of 6-$j$ coefficients in the product of the formula, such as ($F_3$) and ($F_6$), it is clear from the table that both programs compute the same number of 6-$j$ coefficients, so both evaluation algorithms are equally efficient. For test cases where {\tt NJSUMMATION} uses a more efficient summation formula than {\tt NJGRAF}, such as ($F_4$) and ($F_5$), table~1 confirms the idea that this will yield a more efficient evaluation of the formula, with less 6-$j$ coefficients to be computed. Moreover we noticed that for larger $j_i$ values the program {\tt NJGRAF} crashes during the evaluation, giving the error message {\sl `Not a Number'\/} as a result; it is not clear to us why such a crash happens in {\tt NJGRAF}. We also give some results for cases in which {\tt NJGRAF} does not generate a summation formula, such as ($F_7$), ($F_8$) and ($F_9$). Also for these cases the program {\tt NJSUMMATION} is able to evaluate the formula in a reasonable time, such as 25 seconds on the Sun Sparc to compute 1.8 million 6-$j$ coefficients in a seven-fold summation ($F_8$). Table~2 gives results obtained for recoupling coefficients for which independent sums arise, i.e.\ ($G_4$), ($F_{10}$), ($F_{11}$) and ($F_{12}$). Again for each test case both the number of 6-$j$ coefficients computed and the running time are given; running times are now obtained on a 486-based PC running Linux. It is obvious that taking the independent sums into account results in a more efficient evaluation process, with fewer 6-$j$ coefficients to be computed. For large summations the advantage factor is even of several orders of magnitude. For the case ($G_4$) {\tt NJGRAF} too gives results which we can compare with ours. As shown in table~2, the number of 6-$j$ coefficients computed by {\tt NJGRAF} is even less than by our splitting version of {\tt NJSUMMATION}. Still it is of the same order of magnitude. The reason for this difference lies in the fact that, although both {\tt NJGRAF} and {\tt NJSUMMATION} obtain a summation formula with the same number of summation variables $k_i$, these formulae are different. And for the case of all equal $j_i$ values the formula obtained by {\tt NJGRAF} turns out to be evaluated somewhat more efficiently than the formula obtained by {\tt NJSUMMATION}. For the other examples {\tt NJGRAF} gives no results to compare with ours. Only for the case ($F_{10}$), {\tt NJGRAF} can generate a summation formula, which is less efficient than the formula generated by {\tt NJSUMMATION}. However this formula can only be evaluated for the smallest set of $j_i$ values; for larger values {\tt NJGRAF} crashes again (as already described above). For the other recoupling coefficients, ($F_{11}$) and ($F_{12}$), {\tt NJGRAF} does not generate a summation formula. Another interesting feature of the splitting version versus the non-splitting version is the fact that the splitting version is also more efficient for the {\em generation\/} of the formula. E.g.\ for ($F_{11}$) generating the formula by the non-splitting version takes 0.31 seconds, while the splitting version takes only 0.05 seconds; for ($F_{12}$) we get 2.83 seconds versus 0.22 seconds (on the PC). This can easily be explained by the fact that the splitting version also reduces the generation process to the generation of several subformulae. Such a subgeneration involves a search process on a coupling subtree, which is much smaller than the complete coupling tree. The complexity of such a subsearch is several orders of magnitude smaller than the complexity of a search on the complete tree. Table~3 shows a comparison between the parallel version {\tt NJSUMPAR} and the sequential version {\tt NJSUMMATION}. Since the parallel version implements the unoptimised algorithm {\tt UNOPT}, we compare it with a sequential implementation of this unoptimised algorithm, as well as with the optimal sequential version {\tt NJSUMMATION} which implements the algorithm {\tt OPT}. To give an idea of the difference in efficiency between these two algorithms, we also list the number of 6-$j$ coefficients actually computed during the evaluation, for both algorithms. For small summations, such as ($F_4$) and ($F_5$), the difference in the number of 6-$j$ coefficients computed is about a factor of 1.5 to 2, while for larger summations, such as ($F_7$), ($F_8$) and ($F_9$), this is already a factor of about 3. This fact is also clear from the running times (obtained on a single transputer) for the optimised and unoptimised sequential version. When comparing the parallel version (run on a network of 4 transputers, consisting of 1 master and 3 slaves) with the unoptimised sequential version, it is clear that we get a speed-up factor of about 3 on the average, which is a good efficiency. But since optimising the sequential version also gives a speed-up factor, the parallel version gives only a speed-up factor of 2 for smaller summations and a very small speed-up for larger summations, when compared with the optimised sequential version. Therefore we conclude that a parallel implementation using the unoptimised algorithm does not gain much more than simply optimising the sequential algorithm. The remaining question is then~: can the optimised evaluation algorithm be parallelised? This problem is quite non-trivial, since pulling out products of 6-$j$ coefficients introduces a sort of data-dependency between the different recursion steps. However, constructing techniques similar to parallel backtracking, we hope to develop an optimised parallel version resulting in a similar speed-up w.r.t.\ the optimised sequential version. % % ACKNOWLEDGEMENTS % \section*{Acknowledgements} It is a pleasure to thank Prof.\ K.\ Srinivasa Rao (IMSC, Madras, India) for stimulating discussions. This work was partly supported by the EEC (contract no.\ CI1*--CT92--0101). % % BIBLIOGRAPHY % \begin{thebibliography}{99} \bibitem{FPV} V.\ Fack, S.N.\ Pitre and J.\ Van der Jeugt, {\em New efficient programs to calculate general recoupling coefficients, part I~: Generation of a summation formula\/}, Comput.\ Phys.\ Commun.\ {\bf 83} (1994) 275. \bibitem{Bur} P.G.\ Burke, Comput.\ Phys.\ Commun.\ {\bf 1} (1970) 241. \bibitem{Bar} A.\ Bar-Shalom and M.\ Klapisch, Comput.\ Phys.\ Commun.\ {\bf 50} (1988) 375. \bibitem{Lim} P.M.\ Lima, Comput.\ Phys.\ Commun.\ {\bf 66} (1991) 89. \bibitem{Yut} A.P.\ Yutsis, I.B.\ Levinson and V.V.\ Vanagas, The Theory of Angular Momentum (Israel Program for Scientific Translation, Jerusalem, 1962). \bibitem{Bie} L.C.\ Biedenharn and J.D.\ Louck, The Racah-Wigner algebra in Quantum Theory, Encyclopedia of Mathematics and its Applications, Vol.~9, ed.\ G.-C.\ Rota (Addison Wesley, Reading, MA, 1981). \bibitem{Bri} D.M.\ Brink and G.R.\ Satchler, Theory of Angular Momentum (Clarendon Press, Oxford, 1968). \bibitem{Jud} B.R.\ Judd, Operator Techniques in Atomic Spectroscopy (McGraw-Hill, New York, 1963). \bibitem{Elb} E.\ El Baz and B.\ Castel, Graphical Methods of Spin Algebras (Marcel Dekker, New York, 1972). \bibitem{FVR} V.\ Fack, J.\ Van der Jeugt and K.\ Srinivasa Rao, Comput.\ Phys.\ Commun.\ {\bf 71} (1992) 285. \bibitem{Rohl} J.S.\ Rohl, Recursion via Pascal (Cambridge University Press, Cambridge, 1984). \bibitem{Rao81} K.\ Srinivasa Rao, Comput.\ Phys.\ Commun.\ {\bf 22} (1981) 227. \bibitem{Edm} A.\ Edmonds, Angular momentum in quantum mechanics (Princeton University Press, Princeton, 1957). \bibitem{Rac} G.\ Racah, Phys.\ Rev.\ {\bf 62} (1942) 438. \bibitem{SH} N.S.\ Scott and A.\ Hibbert, Comput.\ Phys.\ Commun.\ {\bf 28} (1982) 189. \bibitem{RRC} K.~Srinivasa Rao, V.~Rajeswari and Charles B.~Chiu, Comput.\ Phys.\ Commun.\ {\bf 56} (1989) 231. \bibitem{Linux} Linux (Unix clone for 386/486-based PCs) version 1.0; publicly available via anonymous ftp to {\tt nic.funet.fi} in directory {\tt /pub/OS/Linux}. \bibitem{TC} Turbo C++ (version 1.01) User Guide, Borland International Inc.\ (1990). \bibitem{GnuC} GNU CC version 2.5, Free Software Foundation, Cambridge, MA, USA; publicly available via anonymous ftp to {\tt prep.ai.mit.edu} in directory {\tt /pub/gnu}. \bibitem{ParC} Parallel C (version 2.2.2) User Guide, 3L Ltd., UK (1991). \bibitem{ParF} Parallel Fortran (version 2.1.3) User Guide, 3L Ltd., UK (1990). \end{thebibliography} % % TABLE 1 % \begin{table} \begin{center} \caption{Comparison of results obtained by the sequential version {\tt NJSUMMATION} and the program {\tt NJGRAF}. For each test case we list the results obtained for several sets of all equal $(2j_i+1)$ values, which are shown in column~2, while column~3 shows the corresponding value of the recoupling coefficient. For both programs the first column gives the number of summation variables $k_i$ and the number of 6-$j$ coefficients in the formula used, while the second column shows the number of 6-$j$ coefficients actually computed during the evaluation process, and the third column contains the running time in seconds on a Sun Sparc computer.} \bigskip \bigskip {\small \begin{tabular}{|c|r|r|c|r|r|c|r|r|} \hline & & &\multicolumn{3}{c|}{\tt NJSUMMATION} & \multicolumn{3}{c|}{\tt NJGRAF} \\ \cline{4-9} & $j$'s & value & formula & \#6-$j$'s & time & formula & \#6-$j$'s & time \\ \hline ($F_3$) & 3 & 6.250 000E$-2$ & 2 -- 6 & 27 & 0.00 & 2 -- 6 & 27 & 0.00 \\ & 5 & $-$3.240 837E$-2$ & & 63 & 0.00 & & 63 & 0.00 \\ & 7 & $-$7.638 889E$-2$ & & 114 & 0.01 & & 114 & 0.01 \\ & 9 & $-$1.570 528E$-2$ & & 180 & 0.01 & & 180 & 0.01 \\ & 11 & 1.166 757E$-2$ & & 261 & 0.01 & & 261 & 0.01 \\ \hline ($F_4$) & 9 & $-$1.570 528E$-$2 & 2 -- 6 & 222 & 0.01 & 3 -- 7 & 1137 & 0.07 \\ & 11 & $-$1.166 757E$-$2 & & 326 & 0.02 & & 2022 & 0.10 \\ & 13 & 4.549 656E$-$2 & & 450 & 0.02 & & 3280 & 0.13 \\ & 15 & $-$5.691 117E$-$3 & & 594 & 0.02 & & 4977 & 0.21 \\ & 25 & 1.332 245E$-$2 & & 1 614 & 0.04 & & -- & -- \\ \hline ($F_5$) & 5 & $-$3.781 349E$-$2 & 2 -- 7 & 68 & 0.01 & 4 -- 9 & 894 & 0.04 \\ & 7 & 7.060 185E$-$2 & & 121 & 0.01 & & 3101 & 0.13 \\ & 9 & $-$2.024 054E$-$2 & & 189 & 0.02 & & -- & -- \\ & 15 & $-$1.487 217E$-$2 & & 483 & 0.03 & & -- & -- \\ \hline ($F_6$) & 7 & 6.944 444E$-3$ & 2 -- 7 & 139 & 0.01 & 2 -- 7 & 139 & 0.01 \\ & 9 & $-$4.415 162E$-3$ & & 223 & 0.02 & & 223 & 0.02 \\ & 11 & $-$1.759 274E$-3$ & & 327 & 0.02 & & 327 & 0.02 \\ & 17 & $-$1.380 111E$-4$ & & 759 & 0.03 & & 759 & 0.02 \\ \hline ($F_7$) & 3 & 1.250 000E$-$1 & 4 -- 9 & 150 & 0.02 & -- & -- & -- \\ & 5 & 6.537 455E$-$3 & & 868 & 0.03 & & -- & -- \\ & 7 & $-$9.259 259E$-3$ & & 2 985 & 0.06 & & -- & -- \\ & 9 & 3.567 917E$-3$ & & 7 701 & 0.11 & & -- & -- \\ & 11 & 1.213 792E$-$2 & & 16 594 & 0.25 & & -- & -- \\ & 31 & $-$3.657375E$-$4 & & 952 689 & 19.35 & & -- & -- \\ \hline ($F_8$) & 5 & $-$2.286 994E$-$4 & 7 -- 13 & 37 212 & 0.69 & -- & -- & -- \\ & 7 & 2.700 617E$-$3 & & 340 368 & 5.15 & & -- & -- \\ & 9 & $-$1.299 748E$-$3 & & 1 846 922 & 25.47 & & -- & -- \\ \hline ($F_9$) & 3 & 3.906 250E$-$3 & 9 -- 17 & 7 482 & 0.62 & -- & -- & -- \\ & 5 & $-$2.027 983E$-$3 & & 480 752 & 7.07 & & -- & -- \\ & 7 & $-$7.019 462E$-$4 & & 8 535 793 & 109.86 & & -- & -- \\ \hline \end{tabular} } \end{center} \end{table} % % TABLE 2 % \begin{table} \begin{center} \caption{For some recoupling coefficients in which independent sums arise, we compare results obtained by the non-splitting and the splitting version of the sequential program {\tt NJSUMMATION}, as well as results of the program {\tt NJGRAF}. For each test case we give results obtained for several sets of all equal $(2j_i+1)$ values, which are shown in column~2, while column~3 lists the corresponding value of the recoupling coefficient. For both programs the first column gives the number of summation variables $k_i$ and the number of 6-$j$ coefficients in the formula used. For {\tt NJSUMMATION} the second and third column list the number of 6-$j$ coefficients computed and the running time, for the non-splitting version, while the fourth and fifth column list the same data for the splitting version. The same data is listed for {\tt NJGRAF}. Running times are given in seconds and are obtained on a 486-based PC running Linux.} \bigskip \bigskip {\small \begin{tabular}{|c|r|r|c|r|r|r|r|c|r|r|} \hline & & &\multicolumn{5}{c|}{\tt NJSUMMATION} & \multicolumn{3}{c|}{\tt NJGRAF} \\ \cline{4-11} & & & & \multicolumn{2}{c|}{non-splitting} & \multicolumn{2}{c|}{splitting} & & & \\ \cline{5-8} & $j$'s & value & formula & \#6-$j$'s & time & \#6-$j$'s & time & formula & \#6-$j$'s & time \\ \hline ($G_4$) & 3 & 2.241 607E$-$17 & 3 -- 11 & 117 & 0.01 & 45 & 0.01 & 3 -- 11 & 36 & 0.00 \\ & 5 & $-$1.125 161E$-$04 & & 463 & 0.03 & 103 & 0.01 & & 76 & 0.00 \\ & 7 & 7.822 823E$-$18 & & 1 193 & 0.07 & 185 & 0.01 & & 131 & 0.01 \\ & 9 & $-$1.085 478E$-$04 & & 2 451 & 0.13 & 291 & 0.02 & & 201 & 0.02 \\ \hline ($F_{10}$) & 3 & $-$3.906 250E$-$03 & 4 -- 12 & 192 & 0.01 & 57 & 0.01 & 5 -- 13 & 564 & 0.02 \\ & 5 & 1.050 303E$-$03 & & 1 082 & 0.05 & 137 & 0.01 & & -- & -- \\ & 7 & $-$5.835 262E$-$03 & & 3 672 & 0.17 & 252 & 0.02 & & -- & -- \\ & 9 & 2.466 559E$-$04 & & 9 402 & 0.44 & 402 & 0.03 & & -- & -- \\ \hline ($F_{11}$) & 3 & $-$7.812 500E$-$03 & 6 -- 15 & 930 & 0.08 & 180 & 0.01 & -- & -- & -- \\ & 5 & $-$2.118 683E$-$04 & & 13 962 & 0.57 & 942 & 0.05 & & -- & -- \\ & 7 & $-$7.073 045E$-$04 & & 92 673 & 4.05 & 3 123 & 0.14 & & -- & -- \\ & 9 & $-$5.603 515E$-$05 & & 392 973 & 17.50 & 7 923 & 0.37 & & -- & -- \\ \hline ($F_{12}$) & 3 & $-$1.171 875E$-$02 & 11 -- 22 & 49 585 & 2.60 & 1 735 & 0.11 & -- & -- & -- \\ & 5 & $-$1.495 112E$-$06 & & 8 080 900 & 380.01 & 38 080 & 1.72 & & -- & -- \\ \hline \end{tabular}} \end{center} \end{table} % % TABLE 3 % \begin{table}[p] \begin{center} \caption{Comparison of the parallel version {\tt NJSUMPAR} with the sequential version {\tt NJSUMMATION}. For each test case we list the results obtained for several sets of all equal $(2j_i+1)$ values, which are shown in column~2, while column~3 shows the corresponding value of the recoupling coefficient. We compare running times for the parallel (i.e.\ unoptimised) version, the sequential unoptimised version and the sequential optimised version (i.e.\ the version discussed in table~1). Running times are given in seconds and are obtained on a network of 4 transputers for the parallel version and on a single transputer for the sequential versions. We also list the number of 6-$j$ coefficients actually computed for both the optimised and unoptimised version.} \bigskip \bigskip {\small \begin{tabular}{|c|r|r|r|r|r|r|r|} \hline & & & \multicolumn{1}{c|}{parallel} & \multicolumn{4}{c|}{sequential} \\ \cline{4-8} & & & {\tt UNOPT} & \multicolumn{2}{c|}{{\tt UNOPT}} & \multicolumn{2}{c|}{{\tt OPT}} \\ \cline{4-8} & $j$'s & value & time & \#6-$j$'s & time & \#6-$j$'s & time \\ \hline ($F_4$) & 11 & $-$1.166 757E$-$2 & 0.074 & 456 & 0.205 & 326 & 0.155 \\ & 21 & $-$4.668 661E$-$3 & 0.247 & 1 656 & 0.773 & 1 146 & 0.544 \\ & 31 & 8.127 933E$-$3 & 0.561 & 3 606 & 1.833 & 2 466 & 1.246 \\ \hline ($F_5$) & 9 & $-$2.024 054E$-$2 & 0.055 & 357 & 0.155 & 189 & 0.093 \\ & 15 & $-$1.487 217E$-$2 & 0.138 & 987 & 0.436 & 483 & 0.230 \\ & 31 & 5.938 541E$-$3 & 0.553 & 4 207 & 2.124 & 1 927 & 0.996 \\ & 41 & $-$4.149 020E$-$3 & 0.953 & 7 357 & 4.077 & 3 317 & 1.828 \\ \hline ($F_7$) & 5 & 6.537 455E$-$3 & 0.291 & 1 953 & 0.781 & 868 & 0.397 \\ & 9 & 3.567 917E$-$3 & 2.535 & 19 719 & 7.711 & 7 701 & 3.412 \\ & 11 & 1.213 792E$-$2 & 5.458 & 43 749 & 17.300 & 16 594 & 7.386 \\ \hline ($F_8$) & 3 & $-$9.375 000E$-$2 & 0.721 & 4 277 & 1.812 & 1 585 & 0.786 \\ & 5 & $-$2.286 994E$-$4 & 16.176 & 124 020 & 47.502 & 37 212 &16.952\\ & 7 & 2.700 617E$-$3 & 145.639 & 1 233 596 & 462.954 & 340 368 & 151.947 \\ \hline ($F_9$) & 3 & 3.906 250E$-$3 & 3.662 & 20 672 & 9.070 & 7 482 & 3.765 \\ & 5 & $-$2.027 983E$-$3 & 207.980 & 1 634 346 & 628.240 & 480 752 & 214.418 \\ \hline \end{tabular} } \end{center} \end{table} \end{document}