CSc-387 Parallel Processing 
	
CHAPTER 10 : NUMERICAL ALGORITHMS
=================================

MATRICES
========
                (* T Figure 10.1 *)

Matrix-Vector Multiplication:  Given NxN matrix A and Nx1 vector b. Compute

               C = A*b         (* T Figure 10.3 *)

Matrix Multiplication: Given (n x n) square matrices A and B. Compute

               C = A x B       (* T Figure 10.2 *)

Sequential code for Matrix-Multiplication:
------------------------------------------
   (**    C = A*B  **)
for (i=0; i<n; i++)
  for (j=0; j<n; j++){
    /* Compute C[i,j] as the vector product of A[i,*]B[*,j]  */
    C[i][j] = 0;
    for (k=0; k<n; k++)
         C[i][j] = C[i][j] + A[i][k]*B[k][j] ;
  }

       Tseq = O(n^3)

Parallel :
-----------
With N processors, we can expect a time complexity of   O(N^2)
With N^2 processors, each vector product (C[i,j]) can be computed
                                        in parallel yielding a Tpar = O(N)

- If I/O problem is solved efficiently, these two parallel implementations 
  are both cost-optimal (resulting in linear speedup)

With N^3 processors, each multiplication in vector product (C[i,j]) 
    can be computed in parallel and then reduced  ====>        Tpar = O(logN)

If there are fewer than N processors, then we assign submatrices to 
processors, and use BLOCK MATRIX MULTIPLICATION:

            (* T Figure 10.4  *)          (* T Figure 10.5  *)

s^2 submatrices each having a size of (n/s x n/s).

Sequential Code:
----------------
for (p=0; p<s; p++)
  for (q=0; q<s; q++)
    C[p][q] = 0;               /* clear elements of submatrix */
    for (r=0; r<s; r++)   /* s submatrix mult - incorrect in the textbook */
         C[p][q] = C[p][q] + A[p][r]*B[r][q] ;
  }

       Tseq = s * s * s * (n^3 / s^3) =  O(n^3)


Parallel Matrix Multiplication:  P = N^2
=========================================
With N^2 processors, each vector product (C[i,j]) can be computed
in parallel by a processor yielding a  Tpar = O(N). 
However, each PE must get a row of A and a column of B. 

                   (* T Figure 10.6 *)

Unlike what is given in the textbook, rows and columns
can be distributed to processors of a hypercube in  O(N^2) time:

 1) Rows of A is scattered to N submasters (SM-1) 
    each submaster receiving 1 row              ---  O(N^2)
 2) Columns of B is scattered to N other submasters (SM-2)
    each receiving 1 column                     ---  O(N^2)
 3) Submasters in (SM-1) broadcast their rows 
    in their respective subcubes in parallel    ---  O(NlogN)
 4) Submasters in (SM-2) broadcast their columns
    in their respective subcubes in parallel    ---  O(NlogN)

If I/O time is ignored, then   Tpar = O(N)         (linear speedup)
If I/O time is not ignored, then   Tpar = O(N^2)   (poor performance)

We can also obtain a computational time complexity of O(logN) using
P=N^3 processors:

                    (* T Figure 10.7 *)


Matrix Multiplication Using Submatrices
---------------------------------------
s^2 submatrices each having a size of (n/s x n/s).  
P=s^2 and each PE is responsible computing on submatrix of C.

  1) Each PE gets s submatrices corresponding to one row of A
     and s submatrices corresponding to one column of B 
  2) Each PE computes one submatrix of C --  O(s*(n^3/s^3)) = O(n^3/s^2)

 If I/O time is ignored,   Tpar = O(n^3/s^2)   (* linear speedup *)


Recursive Implementation:     (* T Figure 10.8 *) 
------------------------



MESH (2D-TORUS) IMPLEMENTATION
===============================

Cannon's Algorithm:    
-------------------

- Partition the matrices into n^2  square blocks.
- Initially assign A(i,j) and B(i,j) to P(i,j)
- Let P(i,j) compute C(i,j). 
  To do this, elements of A and B must be aligned in such a way that
  P(i,j) should initially have a pair A(i,*) and B(*,j) whose product
  contributes to the final value of C(i,j). To do this ith row is
  shifted i places to left and jth column is shifted j places up
                         
                       (* T Figure 10.10 *)

- After the initial alignment, if we shift A's to the left, and B's
  upwards in a circular fashion, P(i,j) will always get a unique pair 
  A(i, k) and B(k,j) the product of which contributes to C(i,j).

                       (* T Figure 10.11 *)
 
- after n such circular shifts, processor P(i,j) will compute C(i,j).

- If we assume that we have a 2D torus, each shift can be implemented
  in O(1) routing steps and the total time it takes:

      Tpar = O(n)    on   P=n^2  processors. 

      Speedup = n^3/n = n^2    (linear speedup)

- Since  COSTseq = COSTpar; Cannon's Algorithm is COST-OPTIMAL
- Furthermore, memory requirements for each PE is not more than 3 elements.

Q-1: Could you perform the same algorithm on a hypercube of P=n^2
     processors with the same time complexity?

Q-2: Could you modify Cannon's algorithm such that it can run on  
     P < n^2 processors? For simplicity, you may assume that P = n^2/k^2.


Two dimensional Systolic Array Implementation:    (* T Figure 10.12 *)
----------------------------------------------

Systolic Array Implementation of Matrix-Vector product: (* T Figure 10.13 *)
------------------------------------------------------


10.3. SOLVING A SYSTEM OF LINEAR EQUATIONS
==========================================

     Given A, and b, solve     Ax = b   for n unknowns, x.

Gaussian Elimination:     (* T Figure 10.14 *)
---------------------
 for(i=0; i < n-1; i++)      /* for each row except last */
   for(j=i+1; j < n; j++){   /* step through subsequent rows */
      m = a[j][i]/a[i][i]    /* compute multiplier */
      for (k=i; k<n; k++)    /* modify last n-i-1 elements of row j */
          a[j][k] = a[j][k] - a[i][k] * m;
      b[j] = b[j] - b[i]*m   /* modify righthand side */
   }
                           Tseq = O(n^3)

-----------------------
Parallel Implementation
-----------------------
If  p=n, each PE holds one row of the matrix and operates on that row.

1) P(i) broadcasts Row-i to everyone in the order i=0, 1, 2, ... n-1
2) a processor, P(j) can eliminate its row based on Row-i making its 
   first element a zero, a[j][i]=0, and modifying the remaining
   elements accordingly.

The procedure is repeated with each of P1, P2, P3, ..., P(n-2)
broadcasting elements of their rows. 
 
                  (* T Figure 10.15 *)

Time Analysis:
--------------
Communication: To broadcast (n-i-1) data items takes  O((n-i-1)logp)

    Total broadcast time: SUM(i=0, n-2) (n-i-1)logp = O(n^2*logp)

Computation: SUM(j=1, n-1) (n-j+2) = O(n^2)

          Tpar = O(n^2*logp)

- Notice that efficiency will be low, because all the processors
  before the processor holding row i do not participate in the
  computation again. At each step, less work need to be done. 

Partitioning: 
-------------
If p < n, give each PE  n/p rows. However, if we assign the rows in a
cyclic manner, the workload would be more evenly distributed.
     P0 gets 0, n/p, 2n/p, ... 
     P1 gets 1, n/p+1, 2n/p+1, ... 
     P2 gets 2, n/p+2, 2n/p+2, ... 
       ..............                     (* T Figure 10.18 *)

 Tpar = Tcomm + Tcomp = n^2*logp + n^3/p

because  Tcomp = SUM(j=1, n-1) (n-j+2)*(n/p)  =  n^3/p


ITERATIVE METHODS
=================

Jacobi Iteration - covered earlier. Convergence is guaranteed if the
                   coefficient matrix is Diagonally Dominant. In this
                   case, it is proved to converge in logN steps.

Time for iterative methods depends on two things:
  1) How dense the coefficient matrix is
  2) How many steps does it take to converge? Does it converge at all?

For example, coefficient matrix for solving Laplace's equation has
only 5 non-zeros per row:
                           (* T Figure 10.21 *)

FASTER Convergence Methods
--------------------------
GAUSS-SEIDEL RELAXATION:
 Similar to Jacobi with one important difference:
 To compute Xi, use values of X0, X1, ..., X(i-1) from this iteration
 and values of X(i+1), X(i+2), ..., X(n) from previous iteration.

- has higher convergence rate, however, it is strictly sequential (Why?)

However, in solving Laplace's equation, it can be easily made
parallel by changing the ordering of the equations.

Sequential order of computation:               (* T Figure 10.22 *)

Red-Black ordering for Parallel computation:   (* T Figure 10.23 *)

 forall(i=1; i<n; i++)
   forall(j=1; j<n; j++)
      if((i+j) % 2 == 0)            /* compute red points */
        f[i][j]=.25*(f[i-1][j] + f[i][j-1] + f[i+1][j] + f[i][j+1]);
 forall(i=1; i<n; i++)
   forall(j=1; j<n; j++)
      if((i+j) % 2 != 0)            /* compute black points */
        f[i][j]=.25*(f[i-1][j] + f[i][j-1] + f[i+1][j] + f[i][j+1]);

Comments/Questions:
-------------------
- Would anything go wrong if "forall" construct were not used in this code?
- To be more computationally efficient, avoid using the mod operation
- If we were to assign one PE per mesh point, processors would be idle
  half of the time. Could you think of a solution to prevent that?
  (A: assign one Red and one Black point per PE. Use half as many PEs) 

HIGHER ORDER DIFFERENCE METHODS        (* T Figure 10.24 *)
    use more distant points (i.e. 9-point stencil) in the computation
    in order to increase the accuracy and the speed of convergence.   

MULTIGRID METHOD                (* T Figure 10.25 *)
----------------
First, a coarse grid of points is used with which the iterations converge
pretty quickly. Then, a finer grid is used which includes the points of
the coarser grid. The initial values of extra points is obtained through
interpolation. Computation can alternate between fine and coarse grids.