pcxx_ug (1158314), страница 11
Текст из файла (страница 11)
It turns out that our method is as fast as CMF when all axes are declared \NEWS"and can be a few times faster than CMF when some axes are declared \SERIAL." For example, the speed for accessing a 64 64 double precision array is 5 68 105 bytes/second inpC++ and 3 81 105 bytes/second in CMF when the rst dimension is declared \NEWS"and the second \SERIAL."Besides memory allocation, we also provide methods for inter-element communication forthe Fortran arrays. We design a special pC++ collection for this purpose:::48Collection Fortran : public SuperKernel {public:Fortran(Distribution *D,Align *A);void GetFArray(int index,int offset,FArrayDouble &buffer);void GetFArray(int index1,int index2,int offset,FArrayDouble &buffer);};The Fortran collection is derived from the special root collection SuperKernel. The communication function GetFArray fetches an FArray array from a collection element denotedby the given index(es).A Fortran subroutines is declared as a method of an element class.
An invocation of aFortran subroutine declared this way results in separate invocations of the Fortran subroutineon each processor. The Fortran subroutine is thus similar to a message-passing CMF \nodeprogram" though explicit message-passing is strongly discouraged in the Fortran subroutine.Communication between collection elements is accomplished in the pC++ control program.Since Fortran compilers usually add under scores to subroutine names, the pC++ compilerwill need to generate a \wrapper function" for each Fortran subroutine.By way of illustration, the following shows how a Fortran subroutine is called from pC++utilizing the Fortran interface.extern "C" {//Fortran subroutine is declared externalvoid integrate_seg_(double*,double*,int&,double&);}class Segment {public:double seg_sum;int length;FArrayDouble x, y;Segment();void integrate_seg();};Segment::Segment():x(length),y(length) {//FArrays are allocated in constructor}void Segment::integrate_seg() {// wrapper invokes Fortran subroutine49integrate_seg_(x.d,y.d,length,seg_sum);}void main() {Processors P;Distribution D(64,&P,CYCLIC);Align A(64,"[ALIGN(T[i],D[i])]");Fortran<Segment> F(&D,&A);F.integrate_seg();}The Fortran subroutine is given below:cmfsubroutine integrate_seg(x,y,&length,seg_sum)double precision seg_suminteger lengthdouble precision,array(length)::x,ylayout x(:news), y(:news)seg_sum = SUM(x*y)returnendIn this example, the Fortran subroutine computes the dot-product between two one-dimensionalarrays x and y.
The x and y arrays are passed to the Fortran subroutine as \explicit-shape"arrays. However, the interface does not prevent us from declaring the arrays as \assumedshape" arrays.One shortcoming of the current Fortran interface implementation is that the pC++ control program cannot access Fortran COMMON blocks. This requires that all Fortran globalvariables be declared in the pC++ control program and passed to Fortran subroutines asarguments. Nevertheless, it does allow communication between Fortran subroutines throughCOMMON blocks. It is important to note, however, that pC++ programming model allows more than one pC++ collection element to be allocated on one processor, meaningFortran subroutines called by dierent collection elements will access the same COMMONblock.
Data in the COMMON block are shared by elements on that processor. Potentialrace-condition may arise if a programmer is not careful. The best way to solve this problemis to allocate only one collection element on one processor.At the time of this writing, the pC++ compiler cannot generate the wrapper for theFortran subroutine as well as the extern statement, so a user has to hand code them. In thefuture, a user would write something likeextern "HPF" void Segment::integrate_seg(FArrayDouble &x,FArrayDouble &y,int &length,double &seg_sum);to indicate that integrate seg is a Fortran subroutine and the compiler will generate thenecessary code.505.6 Working with Intel node Fortran on the ParagonWe have also ported the Fortran interface to the Intel Paragon.
Unfortunately, Fortran90type compiler is not available on the Paragon, so we can only experiment with the IntelFortran77 compiler. Since the Fortran interface is designed to be portable, most of whatwe have discussed for the CM-5 still applies for the Paragon. The pC++ program in theprevious section can be directly ported to the Paragon without any change except that thelayout specications are ignored on the Paragon. FArrays become normal C++ arrays sincethere is no vector unit attached to the i860 processor.
Nevertheless, the overloaded operatorsmake array element access conceptually easier. The Fortran90 style subroutine has to bechanged with vector expressions being converted to do-loops.516 The Demo Examples. Ver 1.0+The standard distribution of pC++ comes with a series of Demo \Benchmark" examplesin the subdirectory TestSuite that have been run on a variety of machines. In the sectionsbelow we discuss a number of these codes.6.1 BM-1: Block Grid CGThe rst benchmark illustrate a \typical toy problem", grid computation. The computationconsists of solving the Poisson equation on a two dimensional grid using nite dierenceoperators and a simple conjugate gradient method without preconditioning. Though thisproblem is very simple, it does illustrate important properties of the runtime system associated with an important class of computations.In the program we have used an early prototype of our DistributedGrid collection class.In addition we have also used a common blocking technique that is often used to increasethe computational granularity.
This involves placing a square sub-grid on each processor.Given a problem of size 2 and 2 processors with P<N, there are two ways to do this.The most natural is to dene the grid as being of size N by N and mapping it to a P byP processor array by a BLOCK BLOCK distribution through the template. Unfortunatelyour compiler is not yet powerful enough to do the optimization required to eliminate theelement-to-element access functions on the local sub-grids that will yield vectorizable code.An alternative is to make the grid size P by P and set the grid elements to be sub-grids ofsize M by M; = NP .The collection element then takes the formNPMclass GridBlock{double f[M][M], u[M][M], p[M][M],q[M][M], r[M][M], Ax[M][M];double up_edge[M], lo_edge[M], le_edge[M], ri_edge[M];....};where f is the initial data, u is the solution, and the remaining variables are temporaries usedin the conjugate gradient algorithm or are, as in the case of the edge vectors, buers used tostore boundary values from neighboring grid blocks.
The heart of the algorithm is given bythe subroutine below. (We have overloaded operators -, +, and * for arrays here to simplifythe text presented here. Our next version of the compiler will do this automatically.)void conjugategradient(int num_iters, DistributedGrid<GridBlock> &G){int i;double gamma, alfa, tau, gamnew, beta;G.Ap();// q = A*p52G.r = G.f-G.q;gamma = G.dotprod(R,R);G.p = G.r;for(i = 0; i < num_iters && gamma > 1.0e-12 ; i++){G.Ap();// q = A*ptau = G.dotprod(P,Q);// tau = <p,q>alfa = gamma/tau;G.r = G.r - alfa*G.q;gamnew = G.dotprod(R,R); // gamnew = <r,r>beta = gamnew/gamma;G.u = G.u + alfa*G.p;G.p = G.r + beta*G.p;gamma = gamnew;}}Communication occurs only in the function Ap which applies the nite dierence operatorand in the dotprod function.
In Ap the communication is all based on nearest neighbors andin the dotprod function the communication is all based on tree reductions.6.2 BM-2: A Fast Poisson SolverThe Block Grid CG algorithm represents one type of PDE algorithm. Another commonmethod is to use line oriented communications. For example, a standard fast Poisson solverinvolves applying a Fast Forier Transform based sine transform to each column of the array ofdata representing the right hand side of the equation. This is followed by solving tridiagonalsystems of equations associated with each row.
When this is complete, another sine transformof each column will yield the solution.In this case it is best to view the grid as a distributed vector of vectors, where the vectorclass takes the formclass Vector{double data[N];public:fft();sineTransform();};The distributed vector collection must have a special function for solving tridiagonalsystems of equations. In our case we use a standard cyclic reduction scheme.
This is accomplished by building a collection class DistributedTridiagonal which is a subclass of DistributedVector. This class has a public function cyclicReduction() which takes two parameterswhich correspond to the diagonal and o-diagonal elements of the matrix which are storedin each element.53Collection DistributedTridiagonal: public DistributedArray{public:DistributedTridiagonal(Template *T, Align *A, int n);void cyclicReduction(ElementType *diagonal,ElementType *offdiagonal);...MethodOfElement:ElementType *diag;// the matrix diagonalElementType *offd;// the subdiagonalvoid backSolve(int s);void forwardElimination(int n, int s);...};The cyclic reduction function, shown below, is organized as two phases of log(n) paralleloperations.
The rst phases is the forward elimination, the second phase is the back solvestep. At the rst stage of the elimination, only the even elements participate. At the nextstage only the elements with index a multiple of 4 participate, etc. We use the pC++ vectornotation to select the subsets.void DistributedTridiagonal::cyclicReduction(ElementType *diagonal,ElementType *offdiagonal){int s;int n = this->dim1size-1;this->diag = new ElementType(diagonal);this->offd = new ElementType(offdiagonal);for(s = 1; s < n/2; s = s*2){this[2*s:n-1:2*s].forwardElimination(n, s);}for(s = n/2; s >= 1; s = s/2){this[s:n-1:2*s].backSolve(s);}}Communication only happens within the element function forwardElimination and backSolve. In the forward elimination case, the active elements copy the element data s unitsaway and algebraically eliminate them from the system.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
