The Block Grid CG algorithm represents one type of PDE algorithm. Another common method is to use line oriented communications. For example, a standard fast Poisson solver involves applying a Fast Forier Transform based sine transform to each column of the array of data representing the right hand side of the equation. This is followed by solving tridiagonal systems of equations associated with each row. When this is complete, another sine transform of each column will yield the solution.

In this case it is best to view the grid as a distributed vector of vectors, where the vector class takes the form

class Vector{ double data[N]; public: fft(); sineTransform(); };

The distributed vector collection must have a special function for
solving tridiagonal systems 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
parameters which correspond to the diagonal and off-diagonal elements
of the matrix which are stored in each element.

Collection DistributedTridiagonal: public DistributedArray{ public: DistributedTridiagonal(Template *T, Align *A, int n); void cyclicReduction(ElementType *diagonal, ElementType *offdiagonal); ... MethodOfElement: ElementType *diag; // the matrix diagonal ElementType *offd; // the subdiagonal void backSolve(int s); void forwardElimination(int n, int s); ... };

The cyclic reduction function, shown below, is organized as two
phases of *log(n)* parallel operations. The first phases is
the forward elimination, the second phase is the back solve step.
At the first stage of the elimination, only the even elements
participate. At the next stage only the elements with index a multiple
of 4 participate, etc. We use the pC++ vector notation to select
the subsets.

Communication only happens within the element function forwardElimination and backSolve. In the forward elimination case, the active elements copy the element datavoid 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); } }

The final poisson solver uses two collections *F* and *U* which
represent the initial data and the final solution. Both are of
type distributed tridiagonal of vector which is easily mapped
to a one dimensional template. The main part of the program
is shown below.

void FastPoissonSolve( DistributedTridiagonal<Vector> &U, DistributedTridiagonal<Vector> &F){ Vector diag, offdiag; .... F.sineTransform(U); // U <- sineTform(F); U.cyclicReduction(diag, offdiag); U.sineTransform(F); // F <- sineTform(U); U = F; }

Mon Nov 21 09:49:54 EST 1994