Next: BM-3: The NAS Up: The Demo Examples. Previous: BM-1: Block Grid

BM-2: A Fast Poisson Solver

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.


    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 units away and algebraically eliminate them from the system. In the back solve, the values of the elements that are s units away are copied and are used to compute the local value.

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;
    }



Next: BM-3: The NAS Up: The Demo Examples. Previous: BM-1: Block Grid


beckman@cica.indiana.edu
Mon Nov 21 09:49:54 EST 1994