IntroductionSparse matrices: A user's viewData structures and storageOperations on distributed sparse matricesConstructorsElement-wise matrix arithmeticMatrix multiplicationSparse matrix dense vector multiplicationSparse matrix sparse matrix multiplicationSparse matrix indexing, assignment, and concatenationSparse matrix transposeDirect solvers for sparse linear systemsIterative solvers for sparse linear systemsEigenvalues and singular valuesVisualization of sparse matricesSSCA #2 graph analysis benchmarkScalable data generatorKernel 1Kernel 2Kernel 3Kernel 4Visualization of large graphsExperimental ResultsLooking forward: A next generation parallel sparse libraryConclusionDistributed Sparse Matrices for Very High LevelLanguagesJohn R. GilbertDepartment of Computer ScienceUC Santa BarbaraSteve ReinhardtInteractive SupercomputingViral B. ShahDepartment of Computer ScienceUC Santa BarbaraandInteractive SupercomputingNovember 27, 2007AbstractSparse matrices are first class objects in many VHLLs (very high level languages)used for scientific computing. They are a basic building block for various numerical andcombinatorial algorithms. Parallel computing is becoming ubiquitous, specifically dueto the advent of multi-core architectures. As existing VHLLs are adapted to emergingarchitectures, and new ones are conceived, one must rethink tradeoffs in languagedesign. We describe the design and implementation of a sparse matrix infrastructurefor Star-P, a parallel implementation of the MatlabRprogramming language. Wedemonstrate the versatility of our infrastructure by using it to implement a benchmarkthat creates and manipulates large graphs. Our design is by no means specific toStar-P— we hope it will influence the design of sparse matrix infrastructures in otherlanguages.1Contents1 Introduction 32 Sparse matrices: A user’s view 33 Data structures and storage 44 Operations on distributed sparse matrices 64.1 Constructors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64.2 Element-wise matrix arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . 64.3 Matrix multiplication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74.3.1 Sparse matrix dense vector multiplication . . . . . . . . . . . . . . . 74.3.2 Sparse matrix sparse matrix multiplication . . . . . . . . . . . . . . . 74.4 Sparse matrix indexing, assignment, and concatenation . . . . . . . . . . . . 104.5 Sparse matrix transpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114.6 Direct solvers for sparse linear systems . . . . . . . . . . . . . . . . . . . . . 114.7 Iterative solvers for sparse linear systems . . . . . . . . . . . . . . . . . . . . 124.8 Eigenvalues and singular values . . . . . . . . . . . . . . . . . . . . . . . . . 134.9 Visualization of sparse matrices . . . . . . . . . . . . . . . . . . . . . . . . . 135 SSCA #2 graph analysis benchmark 145.1 Scalable data generator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155.2 Kernel 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155.3 Kernel 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165.4 Kernel 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165.5 Kernel 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175.6 Visualization of large graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . 195.7 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206 Looking forward: A next generation parallel sparse library 227 Conclusion 2421 IntroductionTwo trends have emerged of late in scientific computing. The first one is the adoption of highlevel interactive programming environments such as MatlabR[27], R [22] and Python [34].This is largely due to diverse communities in physical sciences, engineering and social sciencesusing simulations to supplement results from theory and experiments.Computations on graphs combined with numerical simulation is the other trend in sci-entific computing. High performance applications in data mining, computational biology,and multi-scale modeling, among others, combine numerics and combinatorics in a variety ofways. Relationships between individual elements in complex systems are typically modeledas graphs. Links between webpages, chemical bonds in complex molecules, and connectivityin social networks are some examples of such relationships.Scientific programmers want to combine numerical and combinatorial techniques in in-teractive VHLLs, while keeping up with the increasing ubiquity of parallel computing. Adistributed sparse matrix infrastructure is one way to address these challenges. We de-scribe the design and implementation of distributed sparse matrices in Star-P, a parallelimplementation of the MatlabRprogramming language.Sparse matrix computations allow structured representation of irregular data structuresand access patterns in parallel applications. Sparse matrices are also a convenient way torepresent graphs. Since sparse matrices are first class citizens in modern programming lan-guages for scientific computing, it is natural to take advantage of the duality between sparsematrices and graphs to develop a unified infrastructure for numerical and combinatorialcomputing.The distributed sparse matrix implementation in Star-P provides a set of well-testedprimitives with which graph algorithms can be built. Parallelism is derived from opera-tions on parallel sparse matrices. The efficiency of our graph algorithms depends upon theefficiency of the underlying sparse matrix infrastructure.We restrict our discussion to the design and implementation of the sparse matrix in-frastructure in Star-P, trade-offs made, and lessons learnt. We also describe our imple-mentation of a graph analysis benchmark, using Gilbert, Reinhardt and Shah’s “Graph andPattern Discovery Toolbox (GAPDT)”. The graph toolbox is built on top of the sparsematrix infrastructure in Star-P [16, 29, 32, 33].2 Sparse matrices: A user’s viewThe basic design of Star-P and operations on dense matrices have been discussed in earlierwork [8, 20, 21]. In addition to MatlabR’s sparse and dense matrices, Star-P providessupport for distributed sparse (dsparse) and distributed dense (ddense) matrices.The p operator provides for parallelism in Star-P. For example, a random parallel densematrix (ddense) distributed by rows across