Abstract A Flexible Class of Parallel Matrix Multiplication Algorithms

a rXiv:cs /0727v1[cs.PF]14J ul2Efficient cache use for stencil operations on structured discretization grids Michael A.Frumkin and Rob F.Van der Wijngaart ∗Numerical Aerospace Simulation Systems Division NASA Ames Research Center Abstract We derive tight bounds on the cache misses for evaluation of explicit stencil oper-ators on structured grids.Our lower bound is based on the isoperimetrical property of the discrete octahedron.Our upper bound is based on a good surface to volume ratio of a parallelepiped spanned by a reduced basis of the interference lattice of a grid.Measurements show that our algorithm typically reduces the number of cache misses by a factor of three,relative to a compiler optimized code.We show that stencil calculations on grids whose interference lattice have a short vector feature abnormally high numbers of cache misses.We call such grids unfavorable and suggest to avoid these in computations by appropriate padding.By direct measurements on a MIPS R10000processor we show a good correlation between abnormally high numbers of cache misses and unfavorable three-dimensional grids.1Introduction On modern computers the gap between access times to cache and to global memory amounts to several orders of magnitude,and is growing.As a result,improvement in usage of the memory hierarchy has become a significant source of enhancing application performance.Well-organized data traffic may improve performance of a program,without changing the actual amount of computation,by reducing the time the processor stalls waiting for data.Both data location and access patterns affect the amount of data movement in the program,and the effectiveness of the cache.A number of techniques for improvements in usage of data caches have been developed in recent years.The techniques include improvements in data reuse (i.e.temporal locality)[3,4,5,13],improvements in data locality (i.e.spatial locality)[13],and reductions in conflicts in data accesses [4,5,9,10].In practice,these techniques are implemented through code and data transformations such as array padding and loop unrolling,tiling,and fusing.Tight lower and upper bounds on memory hierarchy access complexity for FFT and matrixmultiplication algorithms are given in[8].However,questions concerning bounds on the number of cache misses and how closely current optimization techniques approach those bounds for stencil operators remain open.In this paper we consider improvement of cache usgae through maximizing temporal locality in evaluations of explicit stencil operators on structured discretization grids.Our contribution if twofold.First,we prove lower and upper bounds on the number of cache misses for local operators on structured grids.Our lower bound(i.e.the number of un-avoidable misses)is based on the discrete isoperimetric theorem.Our upper bound(i.e.the achievable number of misses)is based on a cachefitting algorithm which utilizes a special basis of the grid interference lattice.As shown by example,the lower bound can be achieved in some cases.The second contribution is the identification of grids unfavorable dimensions which cause significant increases in cache misses.We provide two characterizations of these unfavorable grids.Thefirst one,derived experimentally,states that the product of all rele-vant grid dimensions is close to a multiple of half the cache size.The second characterization is that the grid interference lattice has a short vector.2Cache model and definitionsWe consider a single-level,virtual-address-mapped,set-associative data cache memory,see [7].The memory is organized in a sets of z lines of w words each.Hence,it can be characterized by the parameter triplet(a,z,w),and its size S equals a∗z∗w words.A cache with parameters(a,1,w)is called fully associative,and with parameters(1,z,w)it is called direct-mapped.The cache memory is used as a temporary fast storage of words used for processing.A word at virtual address A is fetched into a(a(A),z(A),w(A))cache location,where w(A)= A mod w,z(A)=(a/w)mod z,and a(A)is determined according to a replacement policy (usually a variation of least recently used).The replacement policy is not important within the scope of this paper.If a word is fetched,then w−1neighboring words are fetched as well tofill the cache line completely.In practice,a,z,and w are often powers of2in order to simplify computation of the location in cache.For example,the MIPS R10000processor for which we report some measurements in Section6,has a cache with parameters(2,512,4),which makes S equal to 4K double precision words,or32KB.Our lower bound for the minimum number of cache misses that must be suffered during a stencil computation holds for any cache,including fully associative caches.The upper bound shows that a particular number of cache misses can be achieved by choosing a special sequence of computations.A cache miss is defined as a request for a word of data that is not present in the cache at the time of the request.A cache load is defined as an explicit request for a word of data for which no explicit request has been made previously(a cold load),or whose residence in the cache has expired because of a cache load of another word of data into the exact same location in the cache(a replacement load).The definitions of cold and replacement loads match those of cold and replacement cache misses,respectively [4],and if w equals1they completely coincide.If a piece of code featuresφcache misses andµcache loads,it can easily be shown thatµ≤wφ.For a code with good spatial locality we typically haveµ≈wφ.As can be shown by a simple example,no bound of the formφ≤cµ(c constant)can be derived for arbitrary code segments,but if the code implements a non-redundant stencil operation,we haveφ≤|K|µ, where|K|is the total number of points within the stencil.This is shown as follows.Let the stencil operation be written as q(x)=Ku(x),with x∈Ω.HereΩis the(not necessarily contiguous)point set on which array q is evaluated.LetΩ|.The number of cache missesφdoes not exceed the total number of accesses to array u(may included repeated accesses to the same element),which equals|K||Ω|,so|φ≤w,which can be used to bound the number of cache misses in termsof the number of cache loads.3A lower bound for cache loads for local operatorsIn this section we consider the following problem:for a given d-dimensional structured grid and a local stencil operator K,how many cache loads have to be incurred in order to compute q=Ku,where q and u are two arrays defined on the grid.We will provide a lower boundµwhich asserts that,regardless of the order the grid points are visited for the computation of q, at leastµcache loads have to occur.In the next section we provide a cachefitting algorithm for the computation of q whose number of cache loads closely approaches the lower bound.We use the following terminology to describe the operator K.The vectors k1,...,k s defined such that q(x),the value of q at the grid point identified by the vector x,is a function of the values u(x+k1),...,u(x+k s),are called stencil vectors.Locality of K means that the stencil vectors are contained in a cube{x||x i|≤r,i=1,...,d}(r is called the radius of K,and2r+1its diameter).In this section we assume that K contains only the star stencil (i.e.the{0,e1,...,e d,−e1,...,−e d}stencil).A lower bound for cache loads for the star stencil will give us a lower bound for any stencil containing it.Let q be computed in the K-interior R of a rectangular region(a grid)G.We assume that computation of q is performed in a pointwise fashion,that is,at any grid point the value of q is computed completely before computation of the value of q at another point is started.In order to compute the value of q at a grid point x,the values of u in neighbor points of x must be loaded into the cache(a point y is a neighbor of x if y−x is a stencil vector of K).If x is a neighbor of y and u(y)has been loaded in cache to compute q(z)but is dropped from the cache before q(x)is computed,then u(y)must be reloaded,resulting in a replacement load associated with x.To estimate the number of elements,ρ,of array u that must be replaced,we choose a partition of R into a disjoint union of grid regions R i,with R=∪k i=1R i,in such a way that q is computed in all points of R i before it is computed at any point of R i+1,see Figure1.Let B ij be the set of points in R j which are neighbors of R i.Since the star stencil is symmetrical, the B ij are neighbor points of B ji.Because any point of B ij can have at most2d neighbors in B ji,we have the following inequalities:1GFigure1:The boundaries B ij of already computed values of q in a sequence of regions R i. Reloading of some values of u on the boundary of R3results in at least max(|B31|+|B32|−S,0) cache loads.For computation of q in R i we have to replace at leastρi values of u,whereρi equals max i−1j=1|B ij|−S,0 .The total number of replaced values in the course of computing q on the entire grid will be at leastρ,whereρequals B−kS,and B equals k i=1 i−1j=1|B ij|. Summing all the terms|B ij|,taking into account Equation(1)and the fact that|B ii|=0, we get:B≥14dki=1(|δR i|−|δi|)−kS=k σ4d k i=1|δi|≥V4d|D|.(3)We subsequently chooseσin such a way thatσ=|δO(d,t)|=dk=12k d k t k−1 ≥8dS(4)for some t,where O(d,t)is the standard d-dimensional octahedron of radius t(see Appendix A).It follows from Equation21,Appendix A,that t can be chosen in such a way thatσisless than8d(2d+1)S.Now the value ofνcan be estimated using the isoperimetric property of the octahedron(see again Appendix A),namely:ν≤|O(d,t)|.Hence,wefindS8d(2d+1)ν≥|δO(d,t)|d−1(5)where c d equals1/(d(2d+1)2d+2).This gives the following lower bound:ρ≥V|O(d,t)|−1d−1−1d−1−1l+(1−2dd−1 .(7)In general,assuming that the cache associativity a is larger than the diameter of the operator K,the order of this lower bound can not be improved,as shows the following example(remember that our lower bound is valid for a cache with any associativity,including a fully associative cache).Let the spatial extents of a two-dimensional grid be n1and n2, respectively,with n1equal to kS and n2arbitrary,and perform calculations of the star stencil (i.e.r=1)in the following order:do i=0,k*a-1do j=2,n2-1do i1=max(2,1+i*(S/a)),min(n1-1,(i+1)*(S/a))q(i1,j)=u(i1,j)+···end doend doend doSince n1equals kS,all values of q and u having the same value of the second index are mapped into the same cache location within a set.Since a exceeds2r+1,none of the values required for the computation of q will be replaced in the cache,except those at a distance r around the line defined by i1=i*S/a.The total number of elements of u read into the cache for execution of this loop nest will therefore be n1n2+(n2−2)2r(ka−1)−4,which equals n1n2(1−2/n1+2a(1−2/n2)/S).Similar examples in higher dimensions show that the order of our lower bound(Equation7)can not be improved.4An upper bound for cache loads for local operators.Cachefitting algorithmIn order to obtain an upper bound we present a cachefitting algorithm which has a small number of replacements.Wefind a set P of cache conflict-free indices of u and calculate Ku at the points of P.Then we tile the index space of u with P to minimize the total number ofreplacements.For the analysis we assume an cache associativity of one,which is the worst case for replacement loads.Let L be a set in the index space of u having the same image in cache as the index (0,...,0),Figure2.L is a lattice in the sense that there is a generating set{b i},i=1,...,d, such that L is the set of grid points{(0,...,0)+ d i=1x i b i|x i∈Z}.We call this the interference lattice of u.It can be defined as the set of all vectors(i1,...,i d)such thati1+n1i2+n1n2i3+···+n1···n d−1i d≡0mod S.(8)In[4]this lattice is defined as the set of solutions to the cache miss equation.Let P be a fundamental parallelepiped of L∗.For future reference we note that vol(P)= det L=S.The second equality follows form the fact that L has a basis{v i}of the form:v1=S e1,v i=−m i e1+e i,2≤i≤d,m i+1=ij=1n j.(9)Obviously,the vectorsνi satisfy Equation8.Conversely,any vector satisfying Equation8can be represented as a linear combination of v1,...,v d,with coefficients x k=i k for k=2,...,d, and x1=(i1+m2i2+···+m d i d)S−1.x1is an integer number,since i1+m2i2+···+m d i d is divisible by S according to Equation8.Since v1,...,v d are linearly independent vectors, they form a basis of the lattice.Let F be a face of P(see Figure2),and let v be a basis vector of L such that P= {f+x v|f∈F,0≤x<1}.Then shifts F+(k/g)v,k=...,−1,0,1,...contain all integer points of a pencil Q,with Q={f+x v|f∈F,x is any number}for an appropriate value of g†.The values of q at the points of Q can be computed without replacing reusable values of u except at a distance of r or less from the boundary of Q.Let h1,...,h s be the signed projections along F of the stencil vectors of K onto v,and let h+and h−be the maximum and the minimum of the projections,respectively.We assume also that|h+−h−|/g<|v|a, meaning that the extent of P in the direction of v is big enough to allow to compute q on F without replacements.It may be impossible to satisfy this condition when the shortest vector in L is shorter than the diameter of K divided by the cache ttices with short vectors are discussed in Section6.The associated grids are called unfavorable grids.The Cache Fitting Algorithm for computing q is as follows(see Figure2);here K(R)is the set of points where u must be available in order to compute q in all points of R(i.e.the K-extension of R):set w=(1/g)vdo Q=Qmin,Qmaxdetermine face F inside pencil Qdo k=kmin,kmaxload in cache all values of u inside K(F+k∗w)compute q at F+k∗wend doend doIn this algorithm the parameters Qmin,Qmax,kmin and kmax are determined such that the scanning face F sweeps out the entire grid.Whenever a point is not contained in the grid,it is simply skipped in the nest.Since we defined the algorithm in such a way that the scanning face in the direction of v with step g passes through all integer points of Q,the values of q at all points inside Q will be computed.Figure2:The Interference Lattice.Cachefitting set F+k w,k∈Z,sweeps across pencil Q in the direction of v.Only values of u at points at a distance r or less from the pencil boundariesβ1andβ2will be replaced in the cache when K is evaluated inside of Q.Replacements misses can occur only at points at a distance of r or less from the boundaries of the pencils.For each of these points at most s replacement need to take place,where s, the size of the stencil,is defined by s=|K|≤(2r+1)d.So the number of replacements will not exceed r(2r+1)d A,where A is the total surface area of all pencils.To minimize A we choose P so that Q has a good surface to volume ratio.Let P be the fundamental parallelepiped of a reduced basis of L.A basis b1,...,b d of a d-dimensional lattice L is called reduced ifi||b i||≤c d det L,(10)where c d is a constant which depends only on d‡.Let b1be the shortest vector of the basis, and let the eccentricity e of the basis be defined by e=max(||b i||/||b1||).If we define∂P as the surface of P,we can derive an estimation for the surface-to-volume ratio of P: |∂P|det L≤2c d i1||b1||≤ec′d i||b i|| −1d,(11)where we twice used the Hadamard inequality: i||b i||≥det L,and the abovementioned fact that L=S.The constant c′d is defined by c′d=2dc d.Since A does not exceed the surface area of all fundamental parallelepipeds covering the grid,the total number of these parallelepipeds(which equals|G|/det L)gives us:A≤|∂P||G|/det L,so that the total number of replacementsρcan be bounded by:ρ≤r(2r+ 1)d|∂P||G|/det L.This,combined with Equation11,gives an upper bound for the total number of elements to be loaded into the cache in the cachefitting algorithm:µ≤|G|+ρ≤|G| 1+ec′′d S−1f d−1||b1||≤fc d.In Appendix B we show that there are grids whose interference lattices feature f’s that are independent of S(provided that S is a prime power,which is true in most practical cases).For these lattices the relative gap between the upper bound(Equation12)and the lower bound(Equation7)of the previous section goes to zero as S increases.When the cache associativity exceeds the diameter of K,this gap can be closed.In that case a parallelepiped, built on a reduced basis of the interference lattice of the array indices with x d=0,can be swept in the d th coordinate direction,similar to the example at the end of Section3.In general,the cachefitting algorithm gives full cache utilization,in contrast to the algorithm forfinding grid-aligned parallelepipeds devoid of interference lattice points,as proposed in [4].See Table2,[4],where the sizes of blocks without self interference are approximately 20%smaller than S.5Lower and upper bounds for multiple RHS arrays In this section we consider the case where there multiple arrays involved in the computation of q.Let p be the number of arrays(we call these the RHS arrays),all having the same sizes, and let the stencil of each RHS array include the star stencil.This means,in particular, that for each boundary point of any region R i(see Figure1)values of p RHS arrays are necessary§for computation of q in R i.Hence,we have to replace at leastρi values,withρi=max(p( j<i|B ij|)−S,0)values of RHS arrays.Now we can repeat the arguments of Section3,with|V|and|G|replaced by p|V|and p|G|,respectively,and S replaced by⌈S/p⌉, to obtain the following lower bound for the number of cache loads for stencil computations with p RHS arrays:µ≥p|V|+ρ≥p|V| 1+c d S d−1−1l+ 1−2d p −1S ⌉,i=2,...,p.Sweeping through the pencil by units of tile P1in thedirection of v we can compute Ku without any cache conflicts,except on the boundary of the pencils.The number of replacement loads of this algorithm can be estimated similarly to the number of replacement loads of a single-array algorithm,taking into account that for calculation of a value u at any point values of all p RHS arrays in the neighbor points may have to be in cache,thus reducing the effective cache size to[S/p]:µ≤p|G|+ρ≤p|G| 1+ec′′d S d (14) where c′′d is a constant which depends only on d,and e is the eccentricity of L.Figure3:Tiling of a fundamental parallelepiped of a reduced basis of the lattice L.We assume that|h+−h−|≤a|v/p|(a is the cache associativity).The tiling effectively reduces the size of the parallelepiped by a factor of at most2p(since x/p≥⌊x/p⌋≥x/(2p)),and increases the cost of a replacement in the cache per point of the boundary of the pencil by at most a factor of p,since elements of all p RHS arrays will be replaced at the same time. 6Unfavorable array sizesWe have implemented our cachefitting algorithm and compared its actually measured num-ber of cache misses with those of the compiler-optimized code for the corresponding naturally ordered loop nest on a MIPS R10000processor(SGI Origin2000).For comparison we chose a second order difference operator(the common13-point star stencil)an a test set includ-ing three-dimensional grids of sizes40≤n1<100,n2=91,and n3=100(the value of the second dimension was chosen to show a typical picture;that of the third dimension is irrelevant).A plot of measured cache misses for both codes is shown in Figure4.The pro-gram was compiled with options“-O3-LNO:prefetch=0,”using the MIPSpro f77compiler, version7.3.1.1m.The prefetchflag disables the prefetching compiler optimization.Without this option the number of cache misses increases significantly,because the compiler does aggressive prefetching to try to reduce execution time.The upper bounds for the cache misses from the previous sections would suggest that the number of replacement cache misses will increase in the cases where the interference lattice has a very short vector.Very short means that the length is smaller than the diameter of the operator divided by the cache associativity.In this case the self interference would increase significantly.This result suggests how to pad arrays to improve cache performance:the padding should be organized in such a way that the shortest vector in the lattice is not too short,though short enough to minimize the number of pencils(large index of scanning face F).The sweeping is organized such that pencils are as wide as possible(i.e.the smallest total number of pencils),while avoiding—in the case of multiple RHS arrays—tiles that are406080100second array dimension (n )0500100015002000c a c h e m i s s e s (n o r m a l i z ed )natural loop order cache fitting algorithm1Figure 4:Plot of measured cache misses for 40≤n 1<100,n 2=91for 13-point star stencil.The top line corresponds to the naturally ordered nest,optimized by the SGI Fortran com-piler.The bottom line corresponds to our cache fitting algorithm.A typical ratio between the two is 3.5.The large fluctuations correspond to grids with short lattice vectors (n 1=45and n 1=90yield shortest vectors (1,0,1)and (2,0,1),respectively).The fluctuations of cache misses of the cache fitting algorithm for such grids can be so big that their cache misses become more numerous than for the compiler-optimized nest.thinner than the diameter of the stencil operator divided by the cache associativity.To demonstrate these unfavorable grids we again choose the second order stencil and force computations in the nest to follow the natural order ¶.Figure 5a shows the correlation between spikes in the number of cache misses and the presence of a very short vector in the lattice.We call these lattices unfavorable for cache utilization.Arrays having such lattices should be avoided on the target machine.When the shortest vector of the interference lattice is shorter than the diameter of the operator,the number of cache misses sharply increases.The application developer should avoid such unfavorable array sizes,and compilers should avoid the sizes using appropriate padding of array dimensions.Note that similar unfavorable cache effects have been mentioned in [1].nn A)B)n 1n1Figure 5:Plot A shows measured fluctuations of cache misses (above 15%of the upper bound).Plot B shows the interference lattices with short (less than 8in the L 1norm)vectors.Array sizes are 40≤n 1,n 2<100.The plots can be fitted well by hyperbolae defined by n 1n 2=1of the scanning face direction within each pencil(positive or negative sweep direction),and of the visit order of subsequent pencils.This is always possible for a one-dimensional data dependency.Our results can also be extended to arrays that store more than one word per grid point (tensor arrays).The lower bound of Section3for operations with multiple right hand sides immediately applies to tensor arrays.The upper bound of that section also applies,provided the tensor components can be stored as independent subarrays.In a future study we plan to extend the results of this paper to more general implicit operators,to operators on unstructured grids,and to tensor arrays with restricted storage models.We intend to study more closely the dependence of cache misses on the size of the operator’s stencil.We also plan to enhance the presented results by taking into account a secondary cache and TLB,and to formulate bounds for cache misses more directly than through the determination of cache loads.Appendix A:The simplex and the octahedronIn this section we list some basic facts on the number of integer points in the octahedron and simplex.The standard octahedron is defined as:O(d,t)= x∈Z d|d i=1|x i|≤t (15) and the standard simplex as:S(d,t)= x∈Z d|0≤x1,...,x d,d i=1|x i|≤t .(16)If we consider sections of the octahedron by planes x1=k,k=−t,...,t,then for the number of integer points in the octahedron we get the following recurrence relation:|O(d,t)|=|O(d−1,t)|+2t−1k=0|O(d−1,k)|.(17)This relation can be used to prove that|O(d,t)|=dk=02k d k t k (18)and that|δO(d,t−1)|=|O(d,t)−O(d,t−1)|=dk=12k d k t−1k−1 .(19)Also,the relation|δO(d,t)|=|δO(d,t−1)|+|δO(d−1,t)|+|δO(d−1,t−1)|(20)shows that|δO (d,t )|≤(2d +1)|δO (d,t −1)|.(21)For the number of integer points in the simplex we have the following recurrence relation:|S (d,t )|=|S (d −1,t )|+|S (d,t −1)|.(22)This can be used to prove that cf.[6],Table 169,seealso[12],Section5:|S (d,t )|=d k =1 dk t k = d +t d .(23)From Equations 18and 21it follows that |O (d,t )|≤2d |S (d,t )|.Also,since δO (d,t −1)contains at least two nonoverlapping simplices S (d −1,t )and can be covered by 2d such simplices,we see that2|S (d −1,t )|≤|δO (d,t −1)|≤2d |S (d −1,t )|,d ≥2.(24)Hence,if |S (d −1,t )|equals S ,we have for d ≥2:|δO (d,t )||O (d,t )|≥2|S (d −1,t )|d −1≥2−d +1S −1The standard sphere,defined in [12],is the integer point set of minimal surface area for any given number of interior points.Corollary:Since grids with dimensions n i+k i S,i=1,...,d have the same interference lattice for any non-negative integers k i,any grid can be embedded in a favorable larger grid. Proof:Step a:Let a lattice L have a basis of the form of Equation9.Any lattice vector x,whichincludes all basis vectors of L,with L∞norm at most l must be a solution of the following system of inequalities:|x i|≤l,i=2,...,d|Sx1+m2x2+···+m d x d|≤l .(26) Existence of a solution to this system is equivalent to that of the system|x i|≤l,i=2,...,dm2S x d ≤ll d−1(28)for all nonzero x satisfying|x i|≤l,i=2,...,d.If we choose the nonzero integers m i in such a way that|m i−Sµi|≤2for i=2,...,d, then we getm2S x d ≥c′′′d S= S S(29)which shows that Equation26has no integer solutions if l<(c′′′d/(Sd))1/d.Hence,f in Section4can be chosen as:f=d/c′′′d,and the lattice with the basis given by Equation 9has a reduced basis with eccentricity depending only on d.Step b:In order tofind a grid whose interference lattice is L,wefirst sort the m i in order of increasing gcd(m i,S).Since we assume that S=p n,where p is prime,we know that gcd(m i,S)divides gcd(m i+1,S),and the appropriate grid dimensions n i can be found directly by solving the congruencies(n i m i−m i+1)mod S=0.Acknowledgements.We are grateful to professor Leonid Khachijan for some help with the theoretical part of the paper,and to Jerry Yan and Henry Jin for discussions on the practical aspects of the paper. This work is partially supported by the HPCC/CAS NASA program,and was executed under NAS task order A61812D.References[1]D.H.Bailey.Unfavorable Strides in Cache Memory Systems.Scientific Programming,Vol.4,pp.53–58,1995[2]J.W.S.Cassels.Introduction to Diophantine Approximations.Cambridge Univ.Press,1957(Ch.I,Theorem VIII)[3]S.Coleman,K.S.McKinley.Tile Size Selection Using Cache Organization and DataLayout.In Proc.SIGPLAN’95,Conf.on Programming Language Design and Imple-mentation,June1995,pp.279–289[4]S.Gosh,M.Martonosi,S.Malik.Cache Miss Equations:An Analytical Representationof Cache Misses.ACM ICS1997,pp.317–324[5]S.Gosh,M.Martonosi,S.Malik.Automated Cache Optimization using CME DrivenDiagnostics.ACM ICS2000,11p.[6]R.Graham,D.Knuth,O.Patashnik.Concrete Mathematics.Addison-Wesley,1989[7]J.L.Hennessy,puter Organization and Design.Morgan KaufmannPublishers,San Mateo,CA,1994[8]J.W.Hong,H.T.Kung.I/O Complexity:The Red-Blue Pebble Game.IEEE Symposiumon Theoretical Computer Science,1981,pp.326-333[9]G.Rivera,C.W.Tseng.Data Transformations for Eliminating Conflict Misses.PLDI1998,pp.38–49[10]G.Rivera,C.W.Tseng.Eliminating Conflict Misses for High Performance Architec-tures.ACM ICS1998,pp.353–360[11]A.Schrijver.Theory of Linear and Integer Programming.John Wiley and Sons,1986[12]D.L.Wang,P.Wang.Discrete Isoperimetric Problems.SIAM J.Appl.Math.,Vol.32,pp.860–870,1977[13]M.E.Wolf,m.A Data Locality Optimizing Algorithm.In Proc.SIGPLAN’91,Conf.on Programming Language Design and Implementation,June1991,pp.30–44。

Parallel Computing Toolbox 4.3Perform parallel computations on multicore computers and computer clustersParallel Computing Toolbox™ lets you solve computationally and data-intensive problems using MATLAB® andSimulink® on multicore and multiprocessor computers. Parallel processing constructs such as parallel for-loops,distributed arrays, parallel numerical algorithms, and message-passing functions let you implement task- anddata-parallel algorithms in MATLAB at a high level without programming for specific hardware and networkarchitectures. As a result,converting serial MATLAB applications to parallel MATLAB applications requires fewcode modifications and no programming in a low-level language. You can run your applications interactively oroffline, in batch environments.Key Features▪Support for data-parallel and task-parallel application development▪Ability to annotate code segments with parfor(parallel for-loops) and spmd(single program multiple data) forimplementing task- and data-parallel algorithms▪High-level constructs such as distributed arrays, parallel algorithms, and message-passing functions for processinglarge data sets on multiple processors▪Ability to run eight workers locally on a multicore desktop▪Integration with MATLAB Distributed Computing Server for cluster- and grid-based applications that use anyscheduler or any number of workers▪Interactive and batch execution modesDeveloping parallel applications with Parallel Computing Toolbox. The toolbox enables application prototyping on thedesktop with up to eight local workers (left), and, with MATLAB Distributed Computing Server (right), applications canbe scaled to multiple computers on a cluster or a grid or cloud computing service.You can use the toolbox to execute applications on a single multicore or multiprocessor desktop. Without changing the code, you can run the same application on a computer cluster or a grid or cloud computing service (using MATLAB Distributed Computing Server™). Parallel MATLAB applications can be distributed as executables or shared libraries (built using MATLAB Compiler™) that can access MATLAB Distributed Computing Server.Programming Parallel Applications in MATLABParallel Computing Toolbox provides several high-level programming constructs that let you convert your serial MATLAB code to run in parallel on several workers(MATLAB computational engines that run independently of MATLAB clients).Constructs such as parallel for-loops (parfor) and distributed arrays simplify parallel code development by abstracting away the complexity of managing coordination and distribution of computations and data between a MATLAB client and workers, as well as between workers. These constructs function even in the absence of workers, letting you maintain a single version of your code for both serial and parallel execution.You can run an application that uses up to eight workers on your multicore desktop. Using MATLAB Distributed Computing Server, you can run the same application without modification on large resources such as computer clusters or grid and cloud computing services.Taking Advantage of Built-in Parallel Computing Support in MathWorks ProductsKey functions in a growing number of products in the MATLAB and Simulink product families, such as Statistics Toolbox™ and Simulink Design Optimization™, provide built-in support for parallel computing. In the presence of Parallel Computing Toolbox, these functions exploit parallel computing resources by distributing computations across available workers, enabling you to take advantage of parallel computing resources without additional programming effort.Implementing Task-Parallel AlgorithmsYou can parallelize Monte Carlo simulations and other coarse-grained or embarrassingly parallel problems by organizing them into independent tasks(units of work). Parallel for-loops in the toolbox offer one way to distribute tasks across multiple MATLAB workers. Using these loops, you can automatically distribute independent loop iterations to multiple MATLAB workers. The parfor construct manages data and code transfer between the MATLAB client session and the workers. It automatically detects the presence of workers and reverts to serial behavior if none are present.Using parallel for-loops for a task-parallel application. You can use parallel for-loops in MATLAB scripts and functions and execute them both interactively and offline.You can also explicitly program tasks either as MATLAB functions or as MATLAB scripts. You execute the tasks, when specified as functions, by manipulating task and job objects in the toolbox or, when specified as scripts, by using the batch function.Implementing Data-Parallel AlgorithmsFor MATLAB algorithms that require large data set processing, Parallel Computing Toolbox provides distributed arrays and parallel functions that can operate on these arrays. Distributed arrays let you process significantly larger data sets than you could in a single MATLAB session by distributing data across multiple workers and using parallel functions to perform operations on the distributed data. The toolbox also provides functions for moving distributed data to and from MAT-files.You can work with distributed arrays interactively from your MATLAB session as you would with MATLAB arrays. Working with distributed arrays and their associated parallel functions does not require knowledge oflow-level programming tools, such as MPI. The toolbox provides more than 150 parallel functions for operating on distributed arrays, including linear algebra routines based on ScaLAPACK. You can use the familiar operatorsfor arrays in MATLAB to perform indexing, matrix multiplication, and transforms directly on distributed arrays.Programming with distributed arrays. Distributed arrays and parallel algorithms let you create data-parallel MATLAB programs with minimal changes to your code and without programming in MPI.The toolbox also provides the spmd(single program multiple data) construct, which you can use to designate sections of your code to run concurrently across workers participating in a parallel computation. During program execution,spmd automatically transfers data and code used within its body to the workers and, once the execution is complete, brings results back to the MATLAB client session.For explicit, fine-grained control over your parallelization scheme, which also requires explicitly managing interworker coordination, Parallel Computing Toolbox functions provide access to message-passing routines based on the MPI standard (MPICH2), including functions for send, receive, broadcast, barrier, and probe operations.Working in an Interactive Parallel EnvironmentParallel Computing Toolbox extends the MATLAB interactive environment, letting you prototype and develop task- and data- parallel applications in your familiar environment. The matlabpool command allocates a set of dedicated computational resources for you, connecting your MATLAB session to a pool of MATLAB workers that can run either locally on your desktop or on a computer cluster. The command sets up an interactive parallel execution environment in which you can execute your parallel MATLAB code from the command prompt and retrieve results immediately as computations finish.In this environment, parallel for-loops and distributed arrays automatically detect the presence of workers and distribute computations and data between your MATLAB session and the workers. Commands can be collected into MATLAB functions or scripts and submitted for offline execution.Working in Batch EnvironmentsParallel Computing Toolbox lets you use batch environments for offline execution, enabling you to free your MATLAB client session for other activities while you execute large MATLAB and Simulink applications. You can also shut down your MATLAB client session while your application executes and retrieve results later.The batch function lets you execute MATLAB scripts offline. It provides a mechanism for data transfer betweenMATLAB client and worker workspaces, which frees you from explicitly managing data in multiple workspaces.The toolbox also provides job and task objects. These objects provide a lower-level but more general mechanism to execute parallel MATLAB applications in batch environments. Both the batch function and the job and task objects can be used to offload the execution of serial MATLAB programs from a desktop MATLAB session to a worker.Executing MATLAB scripts offline. With the Configurations Manager, you can direct program execution to your workstation using eight local workers or to a cluster using any scheduler and any number of workers.Scaling Up to Clusters, Grids, and Clouds Using MATLAB Distributed Computing ServerParallel Computing Toolbox provides the ability to use up to eight local workers on a multicore or multiprocessor computer using a single toolbox license. When used with MATLAB Distributed Computing Server, you can run an application using any number of workers on large-scale resources such as computer clusters or grid and cloud computing services. The server also supports both interactive and batch workflows.Using MATLAB Compiler, you can build MATLAB applications that use toolbox functions into standalone executables or shared software components and distribute them royalty-free. These executables and libraries can connect to MATLAB Distributed Computing Server workers and perform MATLAB computations on large computing resources.Product Details, Demos, and System Requirements/products/parallel-computingTrial Software/trialrequestSales/contactsalesTechnical Support/support Running a gene regulation model on a cluster using MATLAB Distributed Computing Server. The server enablesapplications developed using Parallel Computing Toolbox to harness computer clusters for large problems. Compiled parallel MATLAB applications can also use the server.Integrating Products into Your Computing EnvironmentMATLAB Distributed Computing Server provides a job manager and directly supports several third-partyschedulers such as Platform LSF®, Microsoft® Windows® HPC Server, Altair PBS Pro®, and TORQUE. Using the Configurations Manager in the toolbox, you can maintain named settings such as scheduler type, path settings,and cluster usage policies. Switching between clusters or schedulers typically requires changing the configuration name only.MATLAB Distributed Computing Server dynamically enables the required product licenses when an application runs on the cluster, based on the user’s profile. As a result, administrators need to manage only the server license on the cluster rather than separate toolbox and blockset licenses for every cluster user.ResourcesOnline User Community /matlabcentral Training Services /training Third-Party Products and Services /connections Worldwide Contacts /contact© 2010 The MathWorks, Inc. MATLAB and Simulink are registered trademarks of The MathWorks, Inc. See /trademarks for a list of additional trademarks. Other product or brand names may be trademarks or registered trademarks of their respective holders.。

f matrices and through a transformation to a product of two symmetric matrices if ∙ ¸ ∙ ¸ 1 −1 0 1 = and = 2 3 −1 2 10. a. Can one use this section’s formulas that indicate the complexity equivalence of multiplication and squaring of integers to show the complexity equivalence of multiplication and squaring of square matrices? b. Show that multiplication of two matrices of order can be reduced to squaring a matrix of order 2 11. Find a tight lower bound class for the problem of finding two closest numbers among real numbers 1 2 ..., . 12. Find a tight lower-bound class for the number placement problem (Problem 9 in Exercises 6.1).
This file contains the exercises, hints, and solutions for Chapter 11 of the book ”Introduction to the Design and Analysis of Algorithms,” 3rd edition, by A. Levitin. The problems that might be challenging for at least some students are marked by B; those that might be difficult for a majority of students are marked by I

MULTIPLE SOURCE DIRECTION FINDING WITH AN ARRAY OF M SENSORS USINGTWO RECEIVERSE.Fishler and H.MesserDepartment of Electrical Engineering–Systems,Tel Aviv University,Tel Aviv69978,Israel,E-mail::{eranf,messer}@eng.tau.ac.il.ABSTRACTMultiple source directionfinding algorithms(e.g.,MUSIC)are applied on simultaneous measurements collected by Msensors.However,practical considerations may dictate us-ing less receivers than sensors,such that the measurementscannot be collected simultaneously.In such cases,data is collected sequentially from the different array elements in a process which is referred to as”time varying preprocessing”[3],or”switching”.In this paper we study multiple source directionfinding(DF)with an array of M>2elements, where only two receivers are available.1.INTRODUCTION AND PROBLEMFORMULATIONDirectionfinding with fewer receivers then sensors via time varying processing is a very important issue[3].In may practical scenarios the number of receivers is considerably less then the number of sensors.Moreover,the tendency is to use the minimum number of receivers possible,that is only two receivers.Reducing the number of receivers results in a cheaper and much more simple designs.We in-vestigate the question of who many sources can be localized with only two receivers and how to perform this task with a computationally efficient algorithm.Assume q far-field narrow band sources impinge on the array,with p>q sensors,from directions{θ1,...,θq}. Using complex signal representation,the vector of received signals can be written as:x(t)=A(θ)s(t)+n(t)(1) where s(t)is the signals slowly varying complex envelope,n(t)is the additive noise,θis the vector of DOA’s,andA(θ)=[a(θ1),...,a(θq]where a(θ)is the array steering vector in directionθ.We denote by(x(t))i the i-th element of vector x(t).We assume that the source signals,s(t),are stationary and ergodic complex-valued,white,Gaussian random pro-cess with zero mean,uncorrelated with the noise process and with correlation matrix R S.The noise is assumed to be stationary and ergodic complex-valued,white,Gaussian random process with varianceσ2.It is well known that the correlation matrix of x(t),denoted by R x(θ),is given by:R x(θ)=A(θ)R S A(θ)H+σ2I(2)where(·)H denote the complex conjugate transpose opera-tion.The problem of estimatingθfrom a set of N snapshots of the array,x(t1),...,x(t N),is usually refereed to as the localization problem.The case of samples which are time dependent linear transformation of the array output is dis-cussed in[3].The resulted model for the measurements is y(t i)=G(t i)x(t i),where G(t i)is the time dependent linear transformation.Note that G(t i)is actually a matrix where the number of rows is the number of receivers in the system at time t i.We are interested in the special case where G(t i)is a 2×p matrix such that each row is a vector with all elements equal to zero and one element equal to1.For the sake of simplicity,we assume that we sample each combination of two elements N times.The total number of samples isp2·N samples,we denote this number by L.This assumption will be eased later.Note,y(t i)=[x(t i)k x(t i)l]T for some k=l∈{1,...,p}2.KNOWN RESULTS AND RELATION TO OTHERPROBLEMSIn[3]the ML estimate for general G(t i)was presented. This procedure involves maximization over all unknown pa-rameters,that isθ,σ2,R S.The maximization problem be-come extremely difficult even for as little as 2sources.The authors presented an ad-hoc approach,the GLS ,which re-duce the complexity of the estimator to a search over q pa-rameters.This approach is still very complicated and the computation load is huge.By noting that the problem can be modeled as the prob-lem of direction finding with time-varying array,one can apply all the results from [4]which includes,among others,the Cramer Rao lower bound (CRLB ).In this paper con-jectures about the complexity of the ML estimator were also presented.This conjectures state that in the general case the ML estimator is not separable.In [2]it was shown that in cases where the number of sensors is smaller then the number or sources the CRLB for the DOA ’s,θ1,...,θq ,does not approach to zero as the SNR approach to infinity.It is well known that when the number of sources is smaller than the number of sensors the CRLB for the DOA ’s approach to zero when the SNR approach to infinity or the number of snapshots approach to infinity.Eigenvector based method for the case of time-varying arrays had been proposed in [1].In this paper two possible Eigenvector based method had been proposed.One is based on interpolating matrix and the other is based on focus ma-trix.Both this methods can not be used in our case due to the big difference in the steering vectors between different time instances.3.THE IDENTAFABILITY PROBLEMIt is well known that under the model (1)and some very weak conditions on the array,one can localize up to p −1sources.The following theorem establish the maximal number of sources that can be localized when using only two receivers.Theorem 1The maximal number of possible sources which can be localized is p −1,that is the problem is identifiable up to p −1parameters.Proof 1Proving that up to p −1sources could be local-ized is equivalent to proving that f (y (t 1),...,y (t L )|θ)=f y (t 1),...,y (t L )|θfor every θ=θ ,wheref (y (t 1),...,y (t L )|θ)is the joint p.d.f.of y (t 1),...y (t L )given the vector of parameters θ.Let ˜y (t i )be a column vector with 2 p2elements,given by:˜y(t i )=[y (t i )Ty (t i +N )T···y (t i +L −N )T ]T(3)It is clear that {˜y (t i )}N i =1contain all the available samplesand thus it is an equivalent representation of the data.It can be easily verified that {˜y (t i )}N i =1are -plex gaussian vectors with block diagonal correlation ma-trix given by the following:(R ˜y )j j (R ˜y )j j +1(R ˜y )j j +1(R ˜y )j +1j +1 = (R x )kk (R x )kl (R x )lk (R x )ll(4)were k,l are the sensors sampled at the j -th switching.It is clear from the structure of R ˜y that a simple one to one mapping between R x and R ˜y exists.Since both y and ˜yare gaussian random vectors it is suffice to prove that for every θ=θR ˜y |θ=R ˜y |θ.Assumethat exsist a pair θ=θ such that R ˜y |θ=R ˜y |θ.Since,as can easily seen from (4),there exists a one to one mappingbetween R ˜y and R x also R x |θ=R x |θ.We found a pair θ=θ such that R x |θ=R x |θ which contradict the factthat we can localize up to p −1sources.Thus R ˜y |θ=R ˜y |θThis theorem provides a very important result,at each time instant we are sampling a sub array of size two which in turn enable us to localize only one source.but coher-ently combining all the results from the sub arrays enables to localize p −1sources,the same number as if we were sampling the all array with p receivers.4.EIGENVECTORS BASED METHODS The ML estimator for θrequires q dimensional search at lest.The eigenvector based method like the MUSIC offers a way to reduce this complexity to one dimensional search.This reduction in complexity might be crucial since still to-day with the most advanced DSP,searching in more then two dimensional can not be performed in real time.We next describe a new eigenvector based procedure which can be used in our problem.We start with the follow-ing equivalent description of the data.Let y (t i )be a column vector with p elements.Let all the elements be equal to zero except,say the k,l elements,which are equal to x (t i )k and x (t i )l respectively,that is k,l are the two array elements which are sampled at time t i .It is clear that defining y (t i )as a p element vector is an equivalent way to describe theavailable data.Denote by ˆR Y =1NL i =1y (t i )y (t i )H the empirical correlation matrix.The expected value of ˆRY is given by the following:R Y =A (θ)R S A (θ)H +σ2I +Λ(5)where Λis a diagonal matrix with diagonal equal to (p −1)·diag (A (θ)R S A (θ)H ).Note that if Λwas equal to zeroone could use the MUSIC or any other eigenvector basedmethod directly on ˆRY for obtaining an estimate of θ.It is quite easy to see that if all the elements of Λare equal.An eigenvector based method of obtaining an esti-mate for θcan still be used.This is since R Y can be divided to a signal subspace and a noise subspace.Two conditions insure that all the elements of Λare equal.1.The sources are uncorrelated,that is R S is a diagonal matrix.2.All the array elements has a constant gain in every direction θ,that is |a (θ)i |are equal for every i .Although this scenario is very important it can be rarely found in practice.In all other cases an eigenvector based method can not be used directly with R Y .We suggest to use the following ad-hoc solution for the problem.By not-ing the R Y (θ)differ from R X (θ)only by multiplication of the diagonal elements by p −1we suggest to divide the di-agonal elements of ˆRY by p −1.Denote by ˆˆR Y the resulted matrix,it can be easily verified that E {ˆˆRY }=R X (θ)and thus ˆˆRY can serve as approximation to R X and thus we can use ˆˆRY with all the eigenvector based methods,e.g MUSIC .We refer to the MUSIC used with our prepro-cessing as MMUSICExample:Consider a uniform linear array with 4omni-directional elements.Assume two equi-power,partially cor-related (ρ=0.25)sources at bearings 0◦,15◦and N =100.In Figure 1a typical spectrum of the MMUSIC is shown,together with a typical spectrum of the MUSIC (when ap-plied on ˆRY without preprocessing).Direction [Deg]C o s t F u n c t i o nFigure 1:Typical MUSIC and MMUSIC cost functionsIt is obvious that this method can be extended to cases when the number of samples taken from each sensor is not equal.Let n i be number of samples taken at the i -th switch-ing.Let ˆRY = y (t i )y (t i )H .It is easily verified that the expectation of ˆRY is given by the following formula:E {ˆR Y }=(A (θ)R S A (θ)H +σ2I )Ψ(6)where (Ψ)ij is the total number of snapshots taken from thei,j sensors simultaneously,also (Ψ)i,i is the total numberof snapshots taken from i -th sensor,anddenote elementby element ing similar preprocessing,ˆRY can be divided element by element with Ψand the resultedmatrix,denoted again by ˆˆR Y can be used with any eigen-value based method.5.SIMULATION STUDYAssume a uniform linear array with 4omni-directional ele-ments.Assume two equi-power sources at bearings 0◦,15◦,and N =100.The following two graphs depicts the MSE of the first source and the probability of detecting two sources for various correlation coefficients as a function of the SNR.SNRP c Figure 2:Probability of detecting two sources as a function of the SNRThe following two graphs depicts the MSE of the first source and the probability of detecting two sources as the a function of the number of snapshots taken and for various correlation coefficient.The SNR was kept fix at 10dB .6.REFERENCES[1]B.Friedlander and A.Zeira,“Eigenstructure-based al-gorithms for direction finding with time-varying ar-SNRM S E [d e g ]Figure 3:The MSE of the First source as a function of theSNRNP cFigure 4:Probability of detecting two sources as a function of the number ofsnapshotsSNRM S E [d e g ]Figure 5:The MSE of the First source as a function of the number of snapshotsrays,”IEEE Trans.on AES.V ol.32,pp.689-701,April 1996.[2]J.Sheinvald,“On Detection and Localization of Mul-tiple Signals by Sensor Arrays,”Ph.D.Dissertaion,Tel Aviv University,Israel.[3]J.Sheinvald and M.Wax,“Direction Finding withFewer Receivers via Time-Varying Preprocessing”,IEEE Trans.on SP .V ol.47,pp.2-10,January 1999.[4]A.Zeira and B.Friedlander,“Direction Finding withTime Varying Arrays”,IEEE Trans.on SP .V ol.43,pp.927-938,1995.[5]M.A.Doron,A.J.Weiss and H.Messer,”MaximumLikelihood direction finding of wide band sources”,IEEE Trans.on SP .V ol.41,pp.411-414,1993.。