












Parallel Solution Methods using Hierarchical Dense Approximation





Ahmed Sameh, Ananth Grama, and Vivek Sarin, 
Purdue University,
Department of Computer Sciences,
West Lafayette, IN 47907

http://www.cs.purdue.edu/people/{sameh, ayg, sarin}


Overview


o Problem Formulation

o Preconditioning

o Dense Methods for Preconditioning

o Preconditioning Results

o Error Control in Dense Methods

o Parallelism in Dense Methods
Problem Formulation




o Arises in simulation of time-dependent Navier-Stokes equations for incompressible fluid flow

o One of the most time consuming steps

o Large scale 3D domains

o Multiprocessing is indispensable

o Robust, parallel preconditioners required

Nature of the System and Preconditioning

	 Linear System
     where
	 	 ,  in which  and  
	 M : Mass matrix; T : Laplace operator  
	 B : Gradient operator (n x k)


Properties of the Linear System
	 
	 Symmetric indefinite (n +ve and k -ve eigenvalues)
	 Typically    (2D)  or    (3D)
	 


Previous Results on Preconditioning

Multilevel Methods [Sarin and Sameh]

	 Obtains a well-conditioned basis for zero-divergence space 
	 Recursive multilevel strategy to coarsen mesh
	 Preconditioning : robust, effective for a wide range of 
	 Efficient parallel implementation
	 Divergence residual satisfied to machine precision



Reference: "An efficient iterative method for the generalized Stokes problem'', 
                   SIAM Journal on Scientific Computing,  Vol. 19, No. 1, pp. 206--226, Jan. 1998
Hierarchical Methods for Solving Dense Linear Systems


Simple Illustrative Case: A is the Laplacian.

Derive solution from potential theory.


 is the potential and G represents the Green's function of the Laplacian
(1/r in three dimensions, -log r in two dimensions).
Hierarchical Methods: Issues


The underlying kernel is dense.

Solution strategy: Use approximations techniques such as 
Fast Multipole Method, Barnes-Hut Method.

Related Issues:

o Computational complexity for unstructured domains.

o Error control for approximation techniques.

o Parallel processing issues.
Adapting Dense Methods to the Preconditioning Problem


Use a dense solver to compute the preconditioner for the matrix A.

The dominant behavior of matrix A is . 

The Green's function of this operator is .





Issue: Implementing boundary conditions?
Implementing Boundary Conditions for Dense Preconditioner


Analogous problem in potential theory: Compute the potential over a domain resulting from a set of given charges provided the boundary potential is prespecified.

Solution strategy:

Assume (unknown) charges residing on the boundary. 

The result of the boundary and internal charges result in the boundary conditions. Use this to compute the values of the unknown boundary charges.

Finally, use these boundary charges along with given internal charges to compute the potential in the interior.


Computational steps: a dense boundary element solve of an n x n system (for n boundary nodes) followed by a dense mat-vec with an n x n system.
Computational Complexity:

o The dense solve corresponds to a dense system of n1/2 in 2-D and n2/3 in 3-D.

o This follows from a simple boundary-interior argument.

o Therefore, a mat-vec with the dense matrix itself would require O(n) computation in 2D and 
O(n4/3) in 3D.

o Using hierarchical approximations, we can reduce this complexity to O(n1/2 log n) in 2D and 
O(n2/3 log n) in 3D.

o Furthermore, we can derive rigorous bounds on the error associated with these solutions.

o We will also demonstrate excellent parallelization properties of these schemes.



Consequently, this forms an ideal preconditioner.
Comparison of Results:
Preconditioning Effect of Dense Solvers:



  Preconditioning of Hierarchical Approximate Techniques
	Incomplete Factorization 	Hierarchical Approximation
297	14	9
653	20	14
1213	25	14
2573	35	16
4953	45	18


Improved Methods for Solving Dense Linear Systems:

Conventional Hierarchical Methods for Particle Simulations

 o Represent the domain hierarchically using a spatial tree






	set of particles	Can they be approximated	If not, divide domain and
	impacting x.	by their center of mass?	recursively apply the same
			criteria to sub-domains.

o Accurate formulation requires O(n2) force computations

o Hierarchical methods reduce this complexity to O(n) or O(n log n)

o Since particle distributions can be very irregular, the tree
	can be highly unbalanced

Commonly Used Hierarchical Methods

Barnes-Hut method

		construct hierarchical representation of domain

		for each particle {
				apply multipole acceptance criteria (MAC) to the root

				if MAC is satisfied then
				compute direct interaction
					else 
				recursively apply to each of its sons
		}
	

		o Multipole acceptance criteria: 

				is (l / d < a)

		d 	: distance of particle from the center of the box
		l 	: length of the box
		a 	: preset constant (normally between 0.7 and 1.0)




Theorem 2. In Barnes-Hut method, a particle interacts with a bounded number of boxes of fixed size. 

Proof: Since    , the centers of all boxes of size   lie within annulus defined by  
and the boxes lie completely within the annulus defined by  
The ratio of the volumes of the annulus and a box gives the following upper bound on the number of boxes of size   : 


Commonly Used Hierarchical Methods

Fast Multipole Method (FMM)


		construct hierarchical representation of domain
		/* top down pass */
		for each box in tree {
			construct well separated box list

			for each box in well separated list 
				translate multipole expansion of box and add
	 	}
		for each leaf node
				translate series to particle position and apply
	

		o Well separatedness criteria: 



Performance and Errors in Hierarchical Methods:


Theorem 1. In Barnes-Hut method, the ratio   for particle-cluster interactions is bounded as follows: 
where   and   are constants.



Proof:   Since a particle p interacts with a box b, but not with its parent box B,
where   and  .  Use triangle inequality    to show that 



Interactions in Hierarchical Methods:




Using Hierarchical Methods for Computing Dense Matrix-Vector Products:




Interactions (direct force computations) are replaced by Gaussian Quadratures.

Far-field interactions are replaced by multipole series representing remote Gauss points.


Each matrix-vector product now becomes a single tree traversal taking O(n log n) time.

This GMRES can be encapsulated in GMRES (or other iterative solvers) to solve dense systems.




Observation: The systems arising in dense methods tend to be well conditioned because of the diagonal dominance of the operator.



Theorem 5. Error in the piecewise approximate B-H method   

This follows naturally from the following results:

- The number of interactions with subdomains at any level are constant.

- The number of subdomain interactions is logarithmic in number of particles.

- The error associated with a single particle-subdomain interaction is constant.




Theorem 3. Suppose that k charges of strengths { , j=1,...,k} are located within a sphere of radius .  Then, for Barnes-Hut method with -criterion for well-separatedness, the error in potential outside the sphere at a distance  from the center of the sphere due to these charges is bounded by
where p is the degree of the truncated multipole expansion such that ,  and 


Theorem 4. Controlling error in B-H     


The theorem defines a growth rate for the polynomial degree with the net charge inside a subdomain such that the total error associated with the subdomain remains constant.

For a uniform distribution, this growth rate can also be expressed in terms of the domain sizes.


Theorem 6. Computational complexity of the piecewise approximate B-H method  
where l is the levels of the hierarchical decomposition. 

For a uniform distribution, l grows as log8n.

For l < p, we can show that the operation count of the new method is within a fraction 7/3 of the fixed-degree multipole method. For smaller values of l this difference is smaller.

For typical values of p (6 - 7 degree approximations), this corresponds to between 256K - 2M node points. Thus, even for very large scale simulations, the improved method is within a small constant off the fixed-degree method while yielding significant improvements in error.



Experimental Results

Number of term evaluations and associated errors:

(alpha = 0.71)	

	n	Error(orig)	Error(new)	Terms(orig)	Terms(new)

	10000	0.012027	0.012027	12293460	12293460
	15000	0.017326	0.010399	24665616	25274289
	20000	0.025982	0.016820	59631516	60443782
	30000	0.036880	0.017386	94953456	96171621
	80000	0.098395	0.019327	253741860	297318482

	ip45K	2.479027	0.334691	70094952 	102045637
	ip82K	2.307508	0.268452	179039880	224057225


Comparison of Errors:



Comparison of Runtimes:



Concurrency in Hierarchical Methods


o Each node evaluation can potentially be executed as a thread.

o Since there may be a large number of nodes, this may lead to too many threads. In general, we build some granularity into the thread by gathering k nodes into a single thread. Here, k is a user supplied parameter.

o The Barnes-Hut tree is a read-only data structure. Therefore, good performance can be achieved if we can build spatial and temporal locality (and the working set size does not exceed local memory).

Cache Performance of Hierarchical Methods:

Spatially proximate particles are likely to interact with nearly identical parts of the tree. Therefore particles must be traversed in a spatial-proximity preserving order.

Cache Performance of Hierarchical Methods


o If the nodes are traversed in a proximity preserving order, significant amounts of computation 
can be performed on data within the cache.

o When additional data is brought into the cache, it is in the form of incremental branches in the 
tree. Since there are also likely to be close to the branches that have already been fetched, the tree 
itself must be stored in a proximity-preserving ordering.

o The same principle can be used to develop out-of-core treecodes and parallel treecodes for 
multi-processing platforms.


o Since the profile of tree nodes changes slowly across particles, significant amounts of computation can be performed before cache misses.

o The tree itself is ordered in a proximity-ordering to ensure spatial and temporal locality.


Minimizing False Sharing

o Use thread-local data as much as possible (especially read-write data).

o Note that the hierarchical tree itself is a read-only data structure in the tree traversal stage.

o For read-write data, make sure that thread-specific global data is appropriately spaced (with 
padding where needed).

Performance of the SGI Origin Implementation:


Problem Instance : Prop; Number of triangles:    140800
Number of nodes:        70439; Number of Gauss points: 6

	Algorithm	Degree	Time	Error

	Original	4	31.83	0.000406
	Improved	4*	33.60	0.000026
	Reference	9	100.81	-
	Exact		> 900.0

Problem Instance : Gripper2; Number of triangles:    185856
Number of nodes:        92918; Number of Gauss points: 6

	Algorithm	Degree	Time	Error

	Original	4	46.40	0.000516
	Improved	4*	50.13	0.000028
	Reference	9	151.21	-
	Exact		> 900.0

Concluding Remarks:


o Hierarchical dense methods can provide very effective preconditioners.

o Approximation techniques for dense methods have excellent error and computational characteristics.

o Hierarchical dense methods have excellent parallel processing characteristics.



... work is under way to extend to more complex kernel functions and to fully estimate the performance of these preconditioned solvers.


Uzawa-type Methods
	 Solve   using iterative method
	 
	 Accelerate by CG
	 
	 Assumption : div-stability   is well-conditioned (steady state)
	 
	 Single-level iterative schemes with suitable preconditioning
	 

Key Issues
	 
	 Two-level nested iteratitve schemes
	 
	 Expensive iterations due to inner iterative solver


Ongoing Work: Simulation of Particulate Flows
