Conference Program

Simplice Donfack (Università della Svizzera italiana ), Stanimire Tomov (University of Tennessee), Jack Dongarra (University of Tennessee), Olaf Schenk (Università della Svizzera italiana)

We present an efficient hybrid CPU/GPU approach that is portable, dynamically and efficiently balances the workload between the CPUs and the GPUs, and avoids data transfer bottlenecks that are frequently present in numerical algorithms. Our approach determines the amount of initial work to assign to the CPUs before the execution, and then dynamically balances workloads during the execution. Then, we present a theoretical model to guide the choice of the initial amount of work for the CPUs. The validation of our model allows our approach to self-adapt on any architecture using the manufacturer's characteristics of the underlying machine. We illustrate our method for the LU factorization. For this case, we show that the use of our approach combined with a communication avoiding LU algorithm is efficient. For example, our experiments on a 24 cores AMD opteron 6172 show that by adding one GPU (Tesla S2050) we accelerate LU up to 2.4 X compared to the corresponding routine in MKL using 24 cores. The comparisons with MAGMA also show significant improvements.

Nicholas Knight (University of California, Berkeley), James Demmel (University of California, Berkeley), Grey Ballard (University of California, Berkeley), Nicholas Knight (University of California, Berkeley), Edgar Solomonik (University of California, Berkeley)

Communication, i.e., data movement, typically dominates the algorithmic costs of computing the eigenvalues and eigenvectors of large, dense symmetric matrices. This observation motivated multistep approaches for tridiagonalization, which perform more arithmetic operations than one-step approaches, but can asymptotically reduce communication.

We study multistep tridiagonalization on a distributed-memory parallel machine and propose several algorithmic enhancements that can asymptotically reduce communication compared to previous approaches. In the first step, when reducing a full matrix to band form, we apply the ideas of the tall-skinny QR algorithm and Householder reconstruction to minimize the number of interprocessor messages (along the critical path). In subsequent steps, when reducing a band matrix to one with narrower bandwidth, we extend Lang's one-step band tridiagonalization algorithm to a multistep approach, enabling a reduction in both messages and on-node data movement. For thicker bands, we propose a pipelined approach using a 2D block cyclic layout (compatible with ScaLAPACK's), where we agglomerate orthogonal updates across bulge chases; for thinner bands, we apply the idea of chasing multiple (consecutive) bulges, previously studied in shared-memory. In addition to asymptotically decreasing algorithmic costs over a wide range of problem sizes, when the input data is so large that its storage consumes all but a constant fraction of the machine's memory capacity, our approach attains known lower bounds on both the number of words moved between processors and on the number of messages in which these words are sent (along the critical path). Moreover, our approach also benefits the sequential case, asymptotically reducing the number of words and messages moved between levels of the memory hierarchy, and, again, attaining lower bounds. We benchmark and compare our approach with previous works like ELPA to evaluate whether minimizing communication improves overall performance.

We also consider the limitations of our approach and discuss possible extensions. First, the tradeoffs of multistep approaches are much greater when eigenvectors are computed, due to the cost of back-transformation, i.e., applying the orthogonal transformations generated by the tridiagonalization steps. Our approach shows, perhaps surprisingly, that it can be beneficial to perform band reduction steps redundantly, on different subsets of processors, in order to reduce communication involved with the back-transformation. In contrast to previous approaches, which perform only one or two tridiagonalization steps in order to minimize the back-transformation costs, our results show that it can be worthwhile to perform more steps, even when eigenvectors are computed. In practice, we navigate the back-transformation tradeoffs with performance modeling and parameter tuning. Second, if more than a constant fraction of memory is unused, other approaches are needed to exploit this memory to further reduce communication and attain the (smaller) lower bounds that apply; we mention our ongoing work for this scenario. Third, we show how our approach extends to many other applications of multistep band reduction, like Hermitian tridiagonalization, bidiagonalization, and Hessenberg reduction.

Sheng Fang (University of Oxford), Raphael Hauser (University of Oxford)

A novel algorithm for leading part SVD computations of large scale dense matrices is proposed. The algorithm is well adapted to the novel computational architectures and satisfies the loosely coupled requirements(low synchronicity and low communication overheads). In contrast to randomized approaches that sample a subset of the data matrix, our method looks at all the data and works well even in cases where the former algorithms fail. We will start with a geometric motivation of our approach and a discussion of how it differs from other approaches. The convergence analysis is split into the derivation of bounds on the local error occurring at individual nodes, and bounds on the global error accumulation. Several variants of the algorithm will be compared, and numerical experiments and applications in matrix optimization problems will also be discussed.

Ren Xiaoguang (National University of Defense Technology)

With the development of the electronic technology, the processors count in a supercomputer reaches million scale. However, the processes scale of a application is limited to several thousands, and the scalability face a bottleneck from several aspects, including I/O, communication, cache access, etc. In this paper, we focus on the communication bottleneck to the scalability of linear algebraic equation solving. We take preconditioned conjugate gradient (PCG) as an example, and analysis the feathers of the communication operations in the process of PCG solver. We find that reduce communication is the most critical issue for the scalability of the parallel iterative method for linear algebraic equation solve. We propose a local residual error optimization scheme to eliminate part of the reduce communication operations in the parallel iterative method, and improve the scalability of the parallel iterative method. Experimental results on the Tianhe-2 supercomputer demonstrate that our optimization scheme can achieve a much signally effect for the scalability of the linear algebraic equation solve.

Maya Neytcheva (Uppsala University), Ali Dorostkar (Uppsala University), Dimitar Lukarski (Uppsala University), Björn Lund (Uppsala University), Yvan Notay (Université libre de Bruxelles), Peter Schmidt (Uppsala University)

We consider large scale computer simulations of the so-called Glacial Isostatic Adjustment model, describing the submerge or elevation of the Earth surface as a response to redistribution of mass due to alternating glaciation and deglaciation periods.

We consider a 2D model that is often used in the geophysics community, where self-gravitation is excluded. Due to the latter, the Earth is modelled as incompressible. The governing equations are where the coefficients , are the Lamé coefficients, is the strain tensor, is the displacement vector, is the so-called pre-stress, is a body force and is the so-called kinematic pressure, introduced in order to be able to simulate fully incompressible materials, for which .

The model is discretized using the finite element method. It gives raise to algebraic systems of equations that are large, sparse, nonsymmetric and indefinite, and have the following saddle point structure, where the pivot block is nonsymmetric.

These linear systems are solved using Krylov subspace iteration methods and a block preconditioner of block-triangular form, that requires inner-outer iterations. Algebraic multigrid is used as an inner solver. The efficiency of solving one system of the above type is crucial as it will be embedded in a time-evolution procedure, where systems with matrices with similar characteristics have to be solved repeatedly.

The simulations are performed using toolboxes from publicly available software packages - deal.ii, Trilinos, Paralution and AGMG. The numerical experiments are performed using OpenMP-type parallelism on multicore CPU systems as well as on GPU. We present performance results in terms of numerical and computational efficiency, number of iterations and execution time, and compare the timing results against a sparse direct solver from a commercial finite element package.

References:
deal.ii: http://www.dealii.org
Trilinos: http://trilinos.sandia.gov/
Paralution: http://www.paralution.com
AGMG: http://homepages.ulb.ac.be/ ynotay/AGMG/

Markus Blatt (HPC-Simulation-Software & Services)

We present a parallel algebraic (AMG) multigrid method based on on-smoothed aggregation for preconditioning the elliptic problem on bounded domains .

Our method is robust for highly variable or even discontinuous coefficients . It is developed for massively parallel supercomputers with a low amount of main memory per core. Due to a greedy heuristic used for the coarsening the method has a low memory footprint and allows for scalable subsurface flow simulations with up to 150 billion unknowns on nearly 300 thousand cores.

The implementation is available as open source within the "Distributed and Unified Numerics Environment" (DUNE).

Seiji Fujino (Kyushu University), Chiaki Itou

We want to solve a linear system of equations as , where is a sparse real matrix, and are solution and right-hand side vectors of order , respectively. We consider to split the matrix as , where , and denote lower, upper triangular matrices and diagonal matrix, respectively. Moreover, we treat with preconditioned linear systems as , where means a preconditioner matrix.

As well known, Eisenstat trick of Symmetric SOR is an efficient preconditioning. In this algorithm, and as splits of coefficient matrix are computed in place of incomplete factorization of matrix . Accordingly it does not need so much computational time in multiplication of preconditioner and matrix . However, it is not suited to parallel computation. In particular, preconditioning of Eisenstat SSOR is not executed efficiently on parallel computers with distributed memory.

In our talk, we will propose Cache-Cache balancing technique for Eisenstat type of preconditioning on parallel computers with distributed memory. “ Cache-Cache” in French means “ hide and seek” in English. The nonzero entries of once hidden for parallelism are found in computation of for balancing the original matrix . Through numerical experiments, we will make clear that the proposed Cache-Cache balancing technique for parallelism outperforms some conventional preconditionings.

Hartwig Anzt (University of Tennessee), Edmond Chow (Georgia Institute of Technology), Jack Dongarra (University of Tennessee)

Based on the premise that preconditioners needed for scientific computing are not only required to be robust in the numerical sense, but also scalable for up to thousands of light-weight cores, the recently proposed approach of computing an ILU factorization iteratively is a milestone in the concept of algorithms for Exascale computing. While ILU preconditioners, often providing appealing convergence improvement to iterative methods like Krylov subspace solvers, are among the most popular methods, their inherently sequential factorization is often daunting. Although improvement strategies like level-scheduling or multi-color ordering render some parallelism, the factorization often remains an expensive building block, especially when targeting large-scale systems. The iterative ILU computation offers an efficient workaround, as it allows for not only fine-grained parallelism, but also excessive asynchronism in the iterative component updates providing the algorithmic characteristics needed for future HPC systems: - The algorithm allows for scaling to any core count, while only the number of nonzero elements in the factorization sets an efficiency-bound. - Communication latencies are of minor impact, and hierarchical distribution of the data allows for reduction of the communication volume. - The asynchronism implies algorithmic error-resilience, as hardware-failures as well as soft errors neither result in the breakdown, nor in the need for restart. This intrinsic fault-tolerance also removes the need for checkpointing strategies or error-resilient communication, both detrimental to performance. At the same time, the operation characteristics of streaming processors like GPUs promote inherently asynchronous algorithms. While implementing linear algebra routines usually requires careful synchronization, the tolerance for asynchronism makes the iterative ILU factorization an ideal candidate for efficient usage of the computing power provided by GPUs. The increasing number of high-performance computing platforms featuring coprocessors like GPUs results in a strong demand for suitable algorithms. In this paper we review the iterative approach of computing the ILU factorization, and the potential of its GPU implementations. We evaluate different approaches that result in either totally asynchronous iteration, or algorithms in which subsets have a pre-defined update order, with iterating asynchronously between the subsets. For those block-asynchronous strategies, the convergence is obviously ensured as they arise as a special case of asynchronous iteration, their performance might benefit from appealing cache reusage on the GPUs. In this case, the granularity of the asynchronism is sallied to a block level, but in contrast to classical block-asynchronous methods, where blocks are usually formed out of adjacent vector components, this approach blocks matrix entries that are adjacent in memory, but not necessarily in the matrix. We compare the performance characteristics of the implementations featuring different asynchronism granularity for several matrix storage formats and matrix sparsity pattern. We also investigate the impact of the asynchronism. Using the shift representing the age of the oldest component used in an update as a control parameter we analyze the trade-off between performance and convergence rate. Careful inspection reveals that the performance-optimal choice of the shift depends primarily on the hardware characteristics, while the properties of the matrix have a minor impact. Although the potential of the asynchronous ILU computation primarily pays off when targeting large-scale problems, we also include a brief comparison against direct ILU factorizations for problems of moderate size.

Lukasz Szustak (Częstochowa University of Technology), Krzysztof Rojek (Częstochowa University of Technology), Roman Wyrzykowski (Czestochowa University of Technology), Pawel Gepner (Intel Corporation)

The multidimensional positive definite advection transport algorithm (MPDATA) belongs to the group of nonoscillatory forward in time algorithms, and performs a sequence of stencil computations. MPDATA is one of the major parts of the dynamic core of the EULAG geophysical model.

The Intel Xeon Phi coprocessor is the first product based on the Intel Many Integrated Core (Intel MIC) architecture. This architecture offers notable performance advantages over traditional processors, and supports practically the same traditional parallel programming model. Intel MIC architecture combines more than 50 cores (200 threads) within a single chip, it has at least 25 MB aggregated L2 cache and 6 GB of on-board memory. The Intel Phi coprocessor is delivered in form factor of a PCIe device.

In this work, we outline an approach to adaptation of the 3D MPDATA algorithm to the Intel MIC architecture. This approach is based on combination of loop tiling and loop fusion techniques, and allows us to ease memory and communication bounds and better exploit the theoretical floating point efficiency of target computing platforms. The proposed approach is than extended on clusters with Intel Xeon Phi coprocessors, using a mixture of MPI and OpenMP programming standards. For this aim, we use a hierarchical decomposition of MPDATA, including level of cluster as well as distribution of computation across cores of Intel Xeon Phi. We discuss preliminary performance results achieved for cluster with Intel Xeon Phi 5120D.

Yves Ineichen (IBM Research - Zurich), Costas Bekas (IBM Research - Zurich), Alessandro Curioni (IBM Research - Zurich)

Over the years graph analytics has become a cornerstone for a wide variety of research areas and applications, and a respectable number of algorithms have been developed to extract knowledge from big and complex graphs. Unfortunately, the high complexity of these algorithms poses a severe restriction on the problem sizes. Indeed, in modern applications such as social or bio inspired networks, the demand to analyze larger and larger graphs now routinely range in the tens of millions and out-reaching to the billions in notable cases.

We believe that in the era of the big data regime accuracy has to be traded for algorithmic complexity and scalability to drive big data analytics. With this goal in mind we developed novel near linear () methods for graph analytics. In this talk we will focus on two algorithms: estimating subgraph centralities (measure for the importance of a node), and spectrograms (density of eigenvalues) of the adjacency matrix of the graph. On one hand we discuss the underlaying mathematical models relaying on matrix trace estimation techniques. On the other hand we discuss the exhibited parallelism on multiple levels in order to attain good parallel efficiency on a large number of processes. To that end we explore the scalability of the proposed algorithms with benchmarks on huge datasets, e.g. the graph of the European street network containing 51 million nodes and 108 million non-zeros. Due to the iterative behavior of the methods, a very good estimate can be computed in a matter of seconds.

Both analytics have an immediate impact on visualizing and exploring large graphs arising in many fields.

Vaclav Hapla (Technical University of Ostrava), Martin Cermak (Technical University of Ostrava), David Horak (Technical University of Ostrava), Alexandros Markopoulos (Technical University of Ostrava), Lubomir Riha (Technical University of Ostrava)

Discretization of most engineering problems, described by partial differential equations, leads to large sparse linear systems of equations. However, problems that involve variational inequalities, such as those describing the equilibrium of elastic bodies in mutual contact, can be more naturally expressed as quadratic programming problems (QPs). A QP takes the form where is a quadratic function of the form and is a feasible set containg only satisfying zero or more prescribed constraints of any of the following forms \begin{align} B_Ex &= c_E,
B_Ix &\leq c_I,
l \leq x &\leq u, \end{align} where is the symmetric positive semidefinite Hessian matrix, is the right hand side vector, is the equality constraint matrix, is the equality constraint right hand side vector, is the inequality constraint matrix, is the inequality constraint right hand side vector, is the lower bound vector, and is the upper bound vector.

We shall introduce here our emerging software stack PERMON. It makes use of theoretical results in advanced quadratic programming algorithms, discretization and domain decomposition methods, and is built on top of the respected PETSc framework for scalable scientific computing, providing parallel matrix algorithms and linear system solvers. The core solver layer of PERMON has arised from splitting the former FLLOP (FETI Light Layer on Top of PETSc) package into several modules: PermonDDM (aka FLLOP) library module for (the algebraical part of) domain decomposition methods (primarily those from FETI family), PermonQP library module for quadratic programming, and PermonAIF which is a pure C array based wrapper of PermonDDM and PermonQP useful for users wanting to take advantage of Permon in their software but do not understand/want/need object-oriented PETSc API. Other layer include application-specific solver modules such as PermonPlasticity, discretization tools such as PermonCube, interfaces allowing calling FLLOP from external discretization software such as libMesh or Elmer. We will present a snapshot of the current state and an outlook for near future.

Alberto Garcia-Robledo (CINVESTAV), Arturo Diaz-Perez (CINVESTAV), Guillermo Morales-Luna (CINVESTAV)

There exists a growing interest in the reconstruction of a variety of technological, sociological and biological phenomena using complex networks. The emerging “science of networks” has motivated a renewed interest in classical graph problems such as the efficient computing of shortest paths and breadth-first searches (BFS). These algorithms are the building blocks for a new set of applications related to the analysis of the structure of massive real-world networks.

All-Sources BFS (AS-BFS) is a highly recurrent kernel in complex network applications. It consists of performing as many full BFS traversals as the number of vertices of the network, each BFS starting from a different source vertex. AS-BFS is the core of wide variety of global and centrality metrics that allow us to evaluate global aspects of the topology of huge networks and to determine the importance of a particular vertex.

The massive size of the phenomena that complex networks model introduce large amounts of processing times that might be alleviated by exploiting modern hardware accelerators, such as GPUs and multicores. However, existing results accelerating BFS-related problems reveal the existing gap between complex network algorithms and current parallel architectures.

There is an useful duality between fundamental operations of linear algebra and graph traversals. For example, in [Qian et al. 2012] it is proposed a GPU approach that overcome the irregular data access patterns on graphs from Electronic Design Automation (EDA) applications by moving BFS to the realm of linear algebra using sparse matrix-vector multiplications (SpMVs). Sparse matrix operations provide a rich set of data-level parallelism, high arithmetic loads and clear data-access patterns, suggesting that an algebraic formulation of BFS can take advantage of the GPU data-parallel architecture.

Unfortunately, there is neither a systematic study on the algebraic BFS approach on non-EDA networks nor a study that compares the algebraic and visitor strategies. Is the algebraic approach adequate for accelerating measurements on complex networks? How they compare on graphs with different structure?

In this paper, we develop an experimental study of the algebraic AS-BFS on GPUs and multicores, by methodologically comparing the visitation and the algebraic approaches on complex networks. We experimented on a variety of synthetic regular, random, and scale-free networks, in order to reveal the impact of the network topology on the GPU and multicore performance and show the suitability of the algebraic AS-BFS for these architectures.

We experimented with sparse and dense graphs with varying diameter, degree distribution and edge density. In order to create graphs with the desired structure, we made use of the well-known Gilbert and Barabási-Albert random graph models, which create uncorrelated random graphs with Poisson and scale-free degree distribution, respectively.

We found that structural properties present in complex networks influence the performance of AS-BFS strategies in the following (decreasing) order of importance: (1) whether the graph is regular or non-regular, (2) whether the graph is scale-free or not, (3) and the graph density. We also found that even though the SpMV approach is not well suited for multicores, the GPU SpMV strategy can outperform the GPU visitor strategy on low-diameter graphs.

In addition, we show that CPUs and GPUs are suitable for complementary kinds of graphs, being GPUs better suited for random sparse graphs. The presented empirical results might be useful for the design of an hybrid network processing strategy that considers how the different algorithmic BFS approaches perform differently on networks with different structure.

Hiroki Tanaka (Kyoto University), Hiroyuki Ishigami (Kyoto University), Kinji Kimura (Kyoto University), Yoshimasa Nakamura (Kyoto University)

Singular value decomposition facilitates to generate orthonormal basis of the row and null spaces. Singular vectors corresponding to non-zero and zero singular values form an orthonormal basis of the row space and the null space, respectively. If singular values computed do not have relative accuracy, we cannot distinguish easily whether singular values are zero or not. Thus, algorithms which can compute “relatively accurate” singular values and singular vectors with high accuracy are required. In order for us to compute orthonormal basis of the row and null spaces, the numerical algorithm must possess the following two properties. First, the algorithm must compute singular vectors with high accuracy. For this purpose, it is desirable that the algorithm is constituted by a sequence of orthogonal transformations. “Divide and conquer” algorithm (D&C algorithm) and Inverse iteration algorithm (II algorithm) sometimes fail to compute singular vectors. Hence, D&C algorithm and II algorithm are not suitable for the purpose. Secondly, the algorithm must be able to compute “relatively accurate” singular values.

The modified QR algorithm proposed by Demmel and Kahan, and the orthogonal qd-algorithm with shift (oqds algorithm) proposed by von Matt is performed by using sequences of orthogonal transformations. It is to be noted that there is no guarantee that the QR algorithm can compute tiny singular values accurately. On the other hand, the oqds algorithm is guaranteed to compute “relatively accurate” singular values theoretically. However, effective implementation of the oqds algorithm to achieve such performance are not known. Thus, in this talk, we propose two techniques for improving the performance of the oqds algorithm.

First, we improve to generate orthogonal transformations. The oqds algorithm employs the Givens rotation and the generalized Givens rotation. The Givens rotation is actualized by the BLAS1 subroutine "drotg", and the generalized Givens rotation by the original subroutine "drotg3". Instead of using "drotg", we can use the "dlartg" subroutine in LAPACK to actualize the Givens rotation. This subroutine is slower but it is more accurate than "drotg". Moreover, the generalized Givens rotation can be redefined by the product of the Givens rotation and an orthogonal transformation. We can apply "dlartg" to the generalized Givens rotation.

Secondly, we consider an improvement of the convergence speed of the oqds algorithm. We incorporate a new shift of origin into the oqd algorithm in order to accelerate its convergence. It is known that the convergence of the oqds algorithm is accelerated more effectively when the shift value is closer to the smallest singular value of the matrix on each iteration. Acceleration of convergence reduces the number of iterations. The smaller number of iterations provides singular value decomposition with high accuracy since accumulation of rounding errors can be suppressed. In the conventional method, Laguerre's and Newton's methods are used to compute the shift value. However, the shift value computed does not sufficiently accelerate the convergence. Thus, we propose a new shift strategy for the oqds algorithm. The new shift strategy compute a bound of the smallest singular value by using binary search. The shift values computed by the proposed shift strategy are much closer to the smallest singular value than that computed by the conventional shift strategy.

The results of numerical experiments indicate that the improved oqds algorithm performs singular value decomposition with higher accuracy than that of the modified QR algorithm. Both the more accurate Givens rotation and the new shift strategy contribute to compute “relatively accurate” singular values and singular vectors with high accuracy.

The improved oqds algorithm can be parallelized using the same technique as that used by the modified QR algorithm implemented in ScaLAPACK.

Vedran Novakovic (University of Zagreb)

We present a hierarchically blocked one-sided Jacobi algorithm for the singular value decomposition (SVD), targeting both a single and the multiple graphics processing units (GPUs). The blocking structure reflects the levels of GPU's memory hierarchy. The algorithm is faster than MAGMA's DGESVD, while retaining the high relative accuracy. To this end, we developed a family of parallel pivot strategies on GPU's shared address space, but applicable also to inter-GPU communication. Unlike common hybrid approaches, our algorithm in a single GPU setting needs a CPU for the controlling purposes only. When required by the problem size, the algorithm, in principle, scales to an arbitrary number of GPU nodes. The scalability is demonstrated by more than twofold speedup for sufficiently large matrices on a Tesla S2050 system with four GPUs vs. a single Fermi card.

Yuki Fujii (Kyoto University), Hiroyuki Ishigami (Kyoto University), Kinji Kimura (Kyoto University), Yoshimasa Nakamura (Kyoto University)

In this talk, we implement the bisection and inverse iteration (BI) algorithm for symmetric band matrices to singular value decomposition (SVD) for large-scale sparse matrices by using a Krylov subspace method.

We consider SVD of () matrices. In problems concerning scientific applications such as data search and rank- () approximation, () singular pairs, i.e. singular values () and the corresponding singular vectors, are often required. For a large-scale sparse matrix, it is difficult to perform SVD directly with respect to execution time. In this case, the target matrix is transformed into a small matrix whose singular values are approximated well to those of the original matrix as the preprocessing is performed by using Krylov subspace methods.

The Golub-Kahan-Lanczos (GKL) algorithm is a preprocessing algorithm based on the Krylov subspace, and generates an approximate bidiagonal matrix from the original matrix. A singular value problem for the bidiagonal matrix () is replaced with an eigenproblem for a symmetric tridiagonal matrix in terms of a extended matrix. From the extended matrix, () singular pairs of the original matrix are obtained. We see that the eigenpair of the extended matrix corresponding to the -th singular pair of the original matrix can be also used for stopping this algorithm, where the number of such an eigenpair is just one. For this purpose, it is suitable to employ the BI algorithm.

The GKL algorithm usually loses the orthogonality of the Krylov subspace owing to computational error. To improve the orthogonality, we need to incorporate reorthogonalization into the GKL algorithm. Such an algorithm is referred to as the reorthogonalized GKL (RGKL) algorithm. Moreover, the block GKL (BGKL) algorithm is a variant of the GKL algorithm, and the corresponding Krylov subspace is generated by using matrix-matrix multiplications. In general, we need to parallelize these algorithms in terms of the operations at each iteration. It is difficult to achieve high-performance parallel processing with the GKL and RGKL algorithm because it is mainly composed of matrix-vector multiplications. Hence, the BGKL algorithm is expected to achieve higher performance in parallel processing than those of the GKL and RGKL algorithm.

Along with the RGKL algorithm, the Reorthogonalized BGKL (RBGKL) algorithm improves the orthogonality of the Krylov subspace from that of the BGKL algorithm. The RBGKL algorithm generates a approximate block bidiagonal matrix (). As with the GKL algorithm, the -th singular pair of the approximate block bidiagonal matrix can also be used for stopping this algorithm. The subject of this talk is to propose efficient computation of the -th singular pair. Conventionally, we perform this computation by transforming the approximate block bidiaognal matrix into a bidiagonal matrix by the Murata method, and computing the -th singular pair with the BI algorithm for the extended matrix. In our approach, we compute the -th singular pair directly for block bidiagonal matrix by using the BI algorithm for a symmetric band matrix in terms of the extended matrix.

To evaluate the proposed implementation, numerical experiments are executed on shared-memory parallel processors. The results show the following properties. In the RBGKL algorithm, the proposed implementation is more suitable for computing the -th singular pair than the conventional implementation.　Then, in many cases, the RBGKL algorithm with the proposed implementation is faster than the RGKL algorithm with the BI algorithm for computing the -th singular pair.

A. Cristiano I. Malossi (IBM Research - Zurich), Yves Ineichen (IBM Research - Zurich), Costas Bekas (IBM Research - Zurich), Alessandro Curioni (IBM Research - Zurich), Enrique S. Quintana-Ortí (Universidad Jaime I)

While Moore's law is expected to remain valid for approximately one more decade, the end of Dennard scaling (i.e., the ability to shrink the feature size of integrated circuit technology while maintaining a constant power density), which partly motivated the transition to multicore architectures around 2004, is now pushing processor design towards “dark silicon” and, in consequence, heterogeneous systems. The scientific community, as responsible for a large fraction of the energy consumed in today's large-scale computing facilities, is not immune to these trends, and a number of research projects now place energy in par with performance among their goals, in a holistic effort for higher energy efficiency all the way from the hardware to the application software.

Along this line, optimizing the 7/13 “dwarfs”, introduced by UC Berkeley in 2006, represents a challenge with a potential tremendous impact on a vast number of scientific/engineering applications and libraries, which are indeed composed of these computational blocks. However, without a careful modeling and a subsequent understanding of the phenomena, the optimization process of these kernels is doomed to fail.

In this talk, we will analyze the interactions occurring in the performance-power-energy triangle for the sparse matrix-vector product (SpMV), a challenging kernel due to its indirect and irregular memory access pattern. SpMV constitutes the key ingredient of the sparse linear algebra dwarf, as well as for many iterative methods for the solution of linear systems arising in, e.g., boundary value problems, eigenvalue problems, finite element methods, economic modeling, web search, and information retrieval, among others.

In the first part of our analysis, we will identify a set of few critical parameters that drive the performance and instantaneous power consumption of a baseline implementation of the SpMV. We will show the link between these parameters and the algorithmic and architecture features, and on top of them we will build a general classification of sparse matrices, which represents a tool to obtain an instantaneous qualitative indication of the performance-power-energy features of SpMV (i.e., a characterizations along these three magnitudes), based only on a few statistical information about the sparse matrix structure, and independent from the original problem.

In a second step, we will use a small synthetic sparse benchmark collection (the training set) to build a general model which is able to provide accurate quantitative estimations of performance-power-energy for any matrix. To prove the robustness of our methodology, we will test the model against the entire University of Florida Matrix Collection, run on two different high-end multithreaded architectures, representative of state-of-the-art multicore technology.

The model and classification presented here can impact large clusters and cloud systems, where the prediction of performance and energy consumption yields important benefits in terms of, e.g., optimal jobs scheduling and costs prediction, among others.

George Delic (HiPERiSM Consulting, LLC)

This study reports on major performance enhancements for the community multi-scale air quality model (CMAQ) chemistry-transport model (CCTM) that add new levels of parallelism and replace the legacy algorithm in the Gear sparse matrix gas-phase chemistry solver. The role of air quality modeling is to study and predict the chemical interaction, atmospheric transport, and deposition of airborne pollutants on metropolitan, regional, continental and hemispheric scales. CMAQ is a major air quality model used in such studies and was developed by the U.S. Environmental Protection Agency (EPA). It has a world-wide user base of over 2,900 recorded downloads. By applying operator splitting, the CCTM implements multiple science modules that describe various physical processes such as advection, diffusion, photolysis, aqueous chemistry and cloud dynamics, gas-phase chemistry, etc. The gas phase chemistry solver in the CCTM is one of the most computationally intensive modules. It allows a choice of various gas-phase solvers that implement different algorithms for integration of a stiff system of ordinary differential equations (ODE), with sparse Jacobians, to describe production and loss of chemical species in reaction mechanisms. When the Gear algorithm is used for the system of ODEs it accounts for over 50% of the wall clock time of a simulation. To improve performance the new version makes two important changes in the standard EPA CMAQ distribution. The first change replaces the sparse matrix solver used for chemical species concentrations. The second modification integrates the new solver into the transit over grid cells so that separate blocks of cells are distributed to different threads. These modifications were successfully applied to CMAQ with the Rosenbrock solver [1] where the same legacy sparse Gaussian linear solver was replaced. For the Gear ODE algorithm this report describes how the sparse solver (FSparse) replaces the legacy JSparse solver method based on the work of Jacobson and Turco [2]. The FSparse solver is a Fortran implementation of Gaussian elimination procedures from the CSparse library of Davis [3]. The parallel processing limitations of the compressed column sparse storage scheme are overcome by transformations to the compressed row storage scheme after the decomposition phase, and before the solve phase. As a result, the revised Gear solver adds both instruction level parallelism and thread level parallelism to the existing distributed memory (message passing) level in the standard EPA release. Observed numerical differences between the two methods are related to the numerical precision achieved and this is discussed in detail. For a 24-hour simulation on a continental U.S.A. grid of 2.3 million cells, results show that with 8 threads the FSparse version of CMAQ is 1.8 times faster than the standard EPA release of CMAQ with the Gear ODE solver (and is also more accurate). This speedup represents 90% of the theoretical maximum of 2.

[1] G. Delic, 12th Annual CMAS conference, October 28-30, 2013. http://www.cmascenter.org/conference/2013/agenda.cfm. [2] M. Jacobson and R.P. Turco (1994), Atmos. Environ. 28, 273-284 [3] T.A. Davis, Direct Methods for Sparse Linear Systems, SIAM, Philadelphia, 2006.

Dan Gordon (University of Haifa), Rachel Gordon (Technion – Israel Institute of Technology), Eli Turkel (Tel Aviv University)

We consider the Helmholtz equation in 3D: , where is the displacement at , and is the wave number. This equation has applications in many fields, such as acoustics and electromagnetism. In geophysics, methods of full wave inversion reconstruct the structure of some volume of the earth from data received by many sensors recording the reflection of sound waves. The frequency domain approach to this reconstruction requires the repeated solution of the Helmholtz equation, in which the RHS function represents an impact source, so it has a small finite support.

Due to dispersion errors (also known as the pollution effect), it is well known that problems with wave propagation are solved more efficiently by using high order accurate methods. In addition to high order accuracy, it is also useful to consider only compact stencils ( in 3D). This has several advantages. Foremost, is that no non-physical boundary conditions need be considered. In addition this reduces the bandwidth of the matrix to be inverted. Another practical advantage is that huge 3D domains require parallel solvers, and compact schemes lend themselves more easily to parallel domain decomposition techniques, such as the CARP-CG algorithm [2] used in this work.

The main difficulty in applying a compact high order finite difference scheme to solve the Helmholtz equation on a geophysical problem is that large geophysical domains consist of different types of materials with different acoustic properties, i.e., the wave number varies within the domain. In this work we use the compact 6th order scheme of [4], which was developed for variable and shown to be far superior to lower order schemes.

However, high order finite difference schemes are not enough for correctly simulating the wave equation. The numerical solution is carried out in a finite domain, so in order to simulate an infinite domain, we need to impose so-called absorbing boundary conditions (ABCs) that will inhibit reflections from the boundaries. One of the first ABCs was formulated by Engquist and Majda. An equivalent BC was later formulated by Higdon [3]. For spherical domains, a series of absorbing boundary conditions was developed in [1]. We have developed several compact high order schemes based on these ABCs, including a novel approach for improved accuracy at the corners.

In order to evaluate the various ABCs, we ran experiments on a domain with a constant and a known analytic solution. The source was simulated by a function of the type for and zero elsewhere. is the amplitude, , and determines the size of the support. Experiments were ran with various frequencies.

It is well known that for constant , the analytic solution is obtained by a convolution of and , where is the Green's function for 3D. The convolution integral was calculated numerically and compared with the solutions obtained by various ABCs.

In order to determine the quality of the ABCs in absolute terms, the ABC results were compared with the solution obtained when Dirichlet boundary conditions were imposed on the exterior boundaries. The values on the boundaries were determined by the numerical convolution of and . Obviously, no ABC can do better than such a Dirichlet BC. We conclude with a geophysical application.

References

[1] Bayliss, Gunzburger // Turkel. SIAM J. Appl. Math., 1982.

[2] Gordon // Gordon. Paral. Comput., 2010.

[3] Higdon. Math. of Comput., 1987.

[4] Turkel, Gordon, Gordon // Tsynkov. J. Comput. Physics, 2013.

Francois Pacull (Claude Bernard University Lyon 1), Damien Tromeur-Dervout, Stephane Aubert

If parallel algebraic multigrid (AMG) methods are widely used for large sparse linear systems, they are not always the method of choice in the case of strongly non-symmetric matrices arising from systems of partial differential equations (PDE), with tightly coupled scalar unknowns. For example, convection-dominated compressible flows discretized with a high-order central scheme usually yield difficult cases for AMG, in part because of the negligible magnitude of diagonal blocks compared to extra-diagonal ones. Also, they may suffer from the strong anisotropy of flow and mesh.

Here is a brief description of the algorithm of interest: we use the GMRES Krylov subspace iterative solver, preconditioned by AMG. However, the operator for the preconditioning system is not the targeted matrix generated with a high-order scheme, but one discretized with a low-order scheme. In order to apply the preconditioner at each step of the outer GMRES loop, a simple aggregation-based AMG is used. A reduced primary matrix, which coefficients depends on the Frobenius norm of the matrix blocks, is constructed in the coarsening setup phase. The same interpolation strategy is used to interpolate each field. For the smoother, the point-block ILU(0) factorization combined with a line ordering is used locally within some Restricted Additive Schwarz (RAS) preconditioned Richardson iterations. The present line reordering technique consists of several steps: build a weighted quotient graph, each vertex corresponding to a diagonal block, then create sequences of strongly coupled vertices, and finally assemble the sequences in order to form a global permutation. This ordering is found to greatly improve the decomposition quality.

A parallel multilevel graph partitioning algorithm is used first along with a weighted graph, in order to distribute the matrix so as to minimize the coupling between the sub-domains and the amount of communication. Restrictions and prolongations steps can be performed in parallel without any communication. However, the smoothing steps and the coarse solver involve some significant communication. The GMRES method is also set as the coarse grid solver along with a RAS/ILU(0) preconditioner. Although the fine low-order matrix is significantly sparser that the higher-order one, the coarse-grid low-order matrix is quite denser that the fine one.

In order to test the robustness of the proposed method, the cases used in the following originate from compressible aerodynamics, more precisely from discrete direct numerical sensitivity analysis, where accurate flow derivatives are sought. However, the presented method apply to any system of coupled system of PDEs for which point-block algorithms apply. Indeed, the approach is purely algebraic, so that it automatically adjusts itself to any complex domains, discontinuous coefficients and anisotropies. The only requirement is that a linear system discretized with a low-order scheme, e.g. first-order upwind, must be generated alongside the targeted higher-order one. This does not represent a difficulty if the Jacobian matrix is generated by automatic differentiation.

Tom Goffrey (University of Exeter), Jane Pratt (University of Exeter), Maxime Viallet (Max Planck Institute for Astrophysics), Isabelle Baraffe (University of Exeter), Chris Geroux (University of Exeter), Mikhail Popov, Rolf Walder (ENS Lyon), Doris Folini (ETH Zürich)

The evolution of stellar interiors is governed by multi-dimensional non-linear phenomena which occur over a large number of dynamical time scales. To efficiently simulate these hydrodynamic effects a time implicit approach is desirable. Furthermore the physical regime is challenging, due to significant variations throughout a single computational domain. Characteristic parameters such as the Mach number, the density and the pressure scale height can vary by up to ten orders of magnitude over the volume of a star at a single point in its evolution. An outstanding challenge for the field of astrophysical hydrodynamics is to completely resolve flows where Mach numbers typically range from unity to . Without an efficient, scalable and parallel method one cannot expect to simulate stellar processes with the necessary resolution in time or space.

In order to effectively utilise the time implicit approach the resulting matrix system should be efficiently solved, both in terms of memory requirements and the number of iterations needed during the matrix solve. Algebraic preconditioners such as incomplete-LU (ILU) factorisation typically provide scalable, “black-box” tools to improve the solver efficiency. Combining an ILU factorisation with a preliminary block-Jacobi preconditioning step can optimize the preconditioning step and result in fewer iterations. We report on testing of this hybrid preconditioning technique in terms of a reduction in the fill level of the ILU step, and in terms of a simplification of diagonal perturbation methods when compared with pure ILU preconditioning. We then expand this study of the hybrid preconditioning technique by showing preliminary results for stellar processes, including accretion, shear mixing, and turbulent convection.

Cosmin Petra (Argonne National Laboratory), Mihai Anitescu (Argonne National Laboratory)

We present a scalable approach that computes in operationally compatible time the optimal energy dispatch under uncertainty for complex energy systems of realistic sizes. Complex energy systems, such as the US power grid, are affected by increased uncertainty of its target power sources, due for example to increasing penetration of wind power coupled with the physical impossibility of very precise wind forecast. The leading optimization-under-uncertainty paradigm for such problems, stochastic programming, requires thousands of simultaneous scenarios, giving problems with billions of variables that need to be solved within an operationally defined time interval.

To address this challenge we have developed PIPS-IPM, a hybrid MPI-SMP parallel solver for stochastic programming problems that uses an interior point method and implements specialized distributed-memory linear algebra. The linear algebra implementation uses a Schur complement technique to obtain a scenario-based parallelization of the computations and, despite the non-trivial decomposition pattern, achieves very high parallel efficiencies on HPC platforms .

To bridge the space between scalability and performance, needed for the real-time solution of stochastic power grid optimization problems, we have recently proposed several algorithmic and implementation advances inside PIPS-IPM. These include a novel incomplete augmented multicore sparse factorization implemented within PARDISO linear solver, coupled with BiCGStab Krylov subspace iterative method to ensure numerical stability and robustness, new multicore- and GPU-based dense matrix kernels, and an efficient implementation of the communication for Cray systems.

We will present large scale simulations of 24-hour horizon dispatch problems with more than 4 billion variables and more than 32 thousands of scenarios on "Titan" Cray XK7 (80% of the machine) of Oak Ridge National Laboratory, "Piz Daint" Cray XC30 (full machine) of Swiss National Computing Centre and "Intrepid" BG/P (80% of the machine) of Argonne National Laboratory. We will discuss several performance metrics of our new developments, such as time-to-solution, peak performance and parallel efficiency, and conclude that solution times compatible with the power grid industry operational practices are possible for this class of optimization problems.

David Horak (Technical University of Ostrava), Vaclav Hapla (Technical University of Ostrava), Lubomir Riha (Technical University of Ostrava), Alexandros Markopoulos (Technical University of Ostrava)

The presentation deals with the ingredients helping us to reach better scalability of the solvers for systems of linear equations and constrained quadratic programming problems. For their solution we develop at IT4Innovations at VSB-TU Ostrava novel software package FLLOP (FETI Light Layer on top of PETSc). It is an extension of PETSc, which is a suite of data structures and routines for the parallel solution of scientific applications modeled by PDE. These solvers can employ direct, iterative or hybrid (combining both previous) methods. Our team has a longterm experience with design and implementation of the scalable domain decomposition methods of FETI type, especially for contact problems, combining iterative a direct solvers. Each method has its pros and cons. We will present bottlenecks of these methods including the strategies for their elimination. A typical example is the coarse problem solution appearing in TFETI method or in deflated conjugate gradient method (DCG). We are able to reach significant improvements using parallel direct solvers in subcommunicators formed by CPUs only or by CPUs accelerated by GPUs/MICs, or using iterative pipelined CG solver applied to the coarse problem matrix better conditioned because of the orthonormalization of the rows of the constrained matrix or the columns of the null space of the stiffness matrix. Numerical experiments illustrate the efficiency of these approaches.

Nerijus Galiauskas (Vilnius University), Julius Žilinskas (Vilnius University)

Finding a step direction in an active-set algorithm for convex quadratic programming problems is the most expensive part of this algorithm. In order to find a step direction we construct and analyze KKT systems. If the KKT system is consistent and has exactly one solution, we solve it by using any solver for systems of linear equations. If the KKT system has more than one solution, we obtain any of them by using a solver for linear least squares problems. If the KKT system is inconsistent, we have to find any eigenvector of the KKT matrix corresponding to the zero eigenvalue. All of these operations (including the analyze of consistency) can be performed by using parallel solvers. Here we present results of numerical experiments related to various strategies how to perform these operations parallelly.

Lukas Pospisil (Technical University of Ostrava), Zdenek Dostal (Technical University of Ostrava)

Numerical solution of many engineering problems, including semicoercive problems of a particle dynamics, leads to the problem of minimizing a convex quadratic function subject to a simple linear inequalities. This problem is given by where is a symmetric positive semidefinite matrix, is a right-hand side vector, and is a vector of bound constraints. To include the possibility that not all the components of are constrained, we admit . In this paper, we propose a modification of recently developed Modified Proportioning with Reduced Gradient Projection (MPRGP) method. This algorithm combines the conjugate gradient steps with the reduced gradient projection steps. Optimality of MPRGP was recently demonstrated on the solution of a class of problems with spectrum of in a given positive interval. Our modification, MPRGP-SPS, converges for any well-posed problem.

The efficiency of our algorithm is illustrated on the solution of granular dynamics problem. Our benchmark consists of moving sphere particles and static box obstacles. In every time-step, we construct the object function from the contacts in the system. Because of the special structure of the matrix , we are also able to construct the projector onto null space. Then, this projector is used in our algorithm.

Michele Martone (Max Planck Institute of Plasma Physics),

Sparse Basic Linear Algebra Subroutines (or `Sparse BLAS') type computations (e.g. triangular solution and multiplication by a vector — respectively SpSV and SpMV) constitute the computational kernel of several scientific and technical computing applications, often being their bottleneck. The memory intensive and irregular features of these computations (foremost, their memory access patterns) make their optimization difficult. An approach that has recently proven to be efficient on shared memory parallel, cache based machines (that is, the `traditional' architectures) is that of a recursive memory layout of the matrix coefficients and coarse grained parallelism in the computations. This approach, that we name `Recursive Sparse Blocks' (RSB) is based on a `cache blocking' like heuristic. In recent work [PARCO_J_14], RSB proved to be competitive with respect to the industry standard Intel MKL CSR implementation and the state-of-the-art CSB prototype, especially when used in SpMV on large symmetric matrices, or with the transposed unsymmetric SpMV operation. Apart from the two operations above, the Sparse BLAS standard specifies also the SpSV operation and provides support for an arbitrary number of right hand sides (`NRHS') and array strides (respectively `INCX' and `INCY') for the two vector operands. An SpMV/SpSV approach that is optimal for the commonest case — one right hand side (NRHS=1) and unitary stride (INCX=INCY=1) — will likely not be optimal for other cases as well, for the memory access pattern will differ substantially. A goal of RSB is to be efficient on the different cases as well, and on the four canonical `BLAS' numerical types.

The present work addresses the problem of tuning our recursive layout for any of the variants of the SpMV and SpSV operations by means of an empirical auto-tuning approach.

That is, given an RSB matrix instance and a SpMV/SpSV operation specification (e.g. NRHS,INCX,INCY), we trade time of search of a `tuned' instance variant for time saved in subsequent, shorter iterations of the same operation. This is possible by re-blocking the matrix differently so it will exhibit a `better' (more efficient) memory access pattern in the given operation.

We present efficient multi-threaded algorithms to support auto-tuning in RSB. On a significant set of real-world matrices, we show that a tangible speedup (over what already exceeds the Intel MKL or the CSB) is still possible. Although results differ widely depending on the matrix structure, type and operation variant, in most of the cases the auto-tuning time can be amortized within a few hundred of iterations.

To confirm the Sparse BLAS oriented spirit of our work, our study extends to all the four the BLAS numerical types, transposition and symmetry and multiple right hand sides variants of the triangular solve and multiply operations.

Since the Sparse BLAS standard interface does not address the case of an `auto-tuning request' by the user, an appropriate extension would be appropriate; we comment this aspect as well.

[PARCO_J_14] http://dx.doi.org/10.1016/j.parco.2014.03.008

Weifeng Liu (University of Copenhagen), Brian Vinter (University of Copenhagen)

Sparse matrix-vector multiplication (SpMV) is perhaps the most widely used non-trivial sparse BLAS routine in computational science and modeling. Compared to CPUs, modern graphics processing units (GPUs) promise much higher peak floating-point performance and memory bandwidth. Thus a lot of research has focused on GPU accelerated SpMV. However, most of them were concentrating on developing new sparse matrix formats (e.g. HYB, Cocktail and BBCOO) constantly challenging well-established scientific software and legacy codes. Further, improving performance of SpMV based on the most widely supported compressed sparse row (CSR) format was largely ignored on GPUs.

Our work described in this paper particularly focuses on accelerating the CSR-format based SpMV on GPUs. Our main contribution is proposing an efficient and general segmented reduction algorithm that directly supports possible empty rows in the input matrix. Note that the empty-row problem actually prevents the classical parallel segmented reduction method from being a general building-block of fast GPU SpMV. Previous work on CSR-format based GPU SpMV either made a strong assumption (i.e. the input matrix is empty-row free) or used pre and post global permutation operations for squeezing the possible empty rows out and inserting them back. Our method, in contrast, deals with the empty rows at runtime and eliminates the costly global permutation operations. Therefore the performance of the CSR-format based SpMV can be significantly improved for general input matrices.

Our algorithm includes 5 steps: (1) The CSR column index array and value array are evenly decomposed to `t' tiles. These tiles are assigned to `b' thread-bunches (i.e. warps or wavefronts). (2) Each GPU thread-bunch executes binary search of `t/b+1' tile boundaries on the CSR row pointer array and obtains corresponding row indices. (3) Each GPU thread-bunch runs `t/b' iterations in a loop to accomplish the following work: (3a) generating localized bit-flag array, (3b) marking tiles to dirty while one or more bit-flag conflicts are detected, (3c) calculating and saving candidate result entries, (3d) executing segmented reduction on them and transmitting the sum of the last segment to next iteration, (3e) storing the segmented sums to known locations in the result vector, and (f) writing the first and the last sum results of all inner iterations to a calibration array of size `2b' in the global memory. (4) The calibration array writes its values to the corresponding entries in the result vector. Then the result vector is numerically correct, except some values generated by these dirty tiles are not in their correct locations. (5) Finally, only for the dirty tiles, our method checks empty rows in the CSR row pointer array and scatter the saved compact results to their correct positions in the result vector. Then the SpMV operation is done.

To evaluate our method, we choose 9 representative unstructured matrices (boyd2, dc2, ASIC_680k, ins2, rajat21, c-73b, transient, wiki-Talk and FullChip) from the University of Florida Sparse Matrix Collection. Compared with state-of-the-art CSR-format based SpMV methods in the nVidia CUSP library (as well as the CUSPARSE library with similar execution pattern and performance) on an nVidia GeForce GTX 680 GPU, our method obtains 5.1x on average and up to 7.9x speedup. Additionally, our method is implemented in OpenCL thus diverse platforms are supported. Compared to our SpMV performance on the nVidia GPU, we achieve 5% on average and up to 83% extra speedup on an AMD Radeon HD 7970 GPU.

Emanuel H. Rubensson (Uppsala University), Elias Rudberg (Uppsala University)

We present a library for efficient block-sparse matrix-matrix multiplication on distributed memory clusters with GPU:s. The library is based on the Chunks and Tasks programming model [Parallel Comput., in press (2013)], allowing it to be implemented without explicitly worrying about communication between the nodes in the distributed memory cluster. On each node, any available GPUs are used for matrix-matrix multiplication tasks. When all GPUs are busy, tasks are also running on the available CPUs, thus making use of the full computing capacity of each node.

Our algorithm is based on a sparse hierarchical data structure [J. Comput. Chem. 28, 2531 ( 2007)] represented by a hierarchy of chunks. The leaves in the hierarchy are block-sparse submatrices. Sparsity may occur at any level in the hierarchy and/or within the submatrix leaves. Both work and data is distributed dynamically over the compute resources, allowing the library to handle unknown and irregular sparsity patterns.

Overlap between computation and communication is achieved for both the data transfer between nodes in the distributed memory cluster and for the data transfer to/from GPUs.

We evaluate the performance of the matrix-matrix multiplication implementation for matrices from linear scaling electronic structure calculations; block-sparse matrices with small block size 16-32, and a priori unknown sparsity patterns. Matrix-matrix multiplication for such matrices is a major computational challenge for large-scale electronic structure calculations.

Pawan Kumar (Fraunhofer ITWM)

In a Schur complement based domain decomposition method, we need to solve for the Schur complement system. The solution for the domain is determined from the solution on the interface, i.e., from the solution of the Schur complement system. For elliptic PDEs under certain conditions, the condition number of the Schur complement system is which is already a much better conditioned system than condition number of the global problem. However, the Schur complement matrix itself is not available explicitly, consequently, some of the approaches such as incomplete factorization that are easily applicable when solving the global linear system is not directly applicable when solving the Shur complement system. Nevertheless, countless preconditioners have been proposed for the Schur complement system.

In this talk, we study a parallel deflation based preconditioning for the Schur complement system for problems with highly discontinuous coefficients. We use the existing FGMRES framework to obtain deflation vectors and update the preconditioner during the iterations. The deflation vectors are the approximate Ritz vectors corresponding to real part of eigenvalues close to zero. We show scalability of our approach using global address space programming framework (GPI) provided by GPI library for inter-node as well as inter-socket communication, and C++-11 threading on local cores residing within each sockets. We use the asynchronous and one sided communications provided by GPI to hide the latencies at various phases of the iterative solver. The parallel efficiency and robustness of the proposed preconditioner is validated for convection diffusion problems discretized using finite difference method. We also briefly discuss the parallel partitioning and matrix assembly of the test problem we consider.

Krzysztof Rojek (Częstochowa University of Technology), Lukasz Szustak (Częstochowa University of Technology), Roman Wyrzykowski (Częstochowa University of Technology)

The general purpose programming on graphics devices is currently supported by two main vendors of GPUs. AMD and NVIDIA deliver powerful platforms, which peak performance exceeds a few TFlop/s. The goal of our research is to develop an efficient and portable adaptation of stencil-based 3D MPDATA algorithm to both types of GPUs, including NVIDIA GeForce GTX Titan and AMD FirePro W8000. Differences between these architectures force users to develop distinct approaches and optimizations for each of them, including algorithm configuration, memory management and organization of computation. In this paper, we propose a performance model, which allows for the efficient distribution of computation across GPU resources. Since MPDATA is strongly memory-bounded, the main challenge of providing a high-performance implementation is to reduce GPU global memory transactions. For this purpose, our performance model ensures a comprehensive analysis of transactions based on local memory utilization, sizes of halo areas (ghost zones), data dependencies between and within stencils. The results of analysis performed using the proposed model are number of GPU kernels, distribution of stencils across kernels, as well as sizes of work-groups for each kernel. The best solution selection should be based on maximizing the GPU occupancy, and minimizing resource requirements simultaneously. However, such a solution does not exist, as the optimization criteria are conflicting. So to solve the problem of optimal mapping stencils onto kernels, we base on minimization of GPU memory transactions with constraints taking into account the GPU occupancy. The proposed model is one of the key components for developing an autotuning mechanism dedicated not only to the MPDATA computations, but also to the whole family of stencil algorithms. This mechanism aims at providing the performance portability for OpenCL codes across different GPU architectures.

Jens Saak (Max Planck Institute for Dynamics of Complex Technical Systems), Martin Köhler (Max Planck Institute for Dynamics of Complex Technical Systems)

The Basic Linear Algebra Subprogramms (BLAS) are the core pieces of all scientific computing codes that employ linear algebra operations. Their improvement has been the focus of many researchers and hardware manufacturers since their introduction in late 1980s. Today highly optimized versions for almost any plattform do exist. Still, especially for debugging or profiling purposes switching between the netlib.org reference implementation and optimized versions is a common task that can be rather time consuming when many parts of the project rely on the BLAS.

Reinventing basic ideas of the liftracc project (https://github.com/pc2/liftracc) with an entirely new focus, we present a wrapper library and tools for easy switching and basic profiling of BLAS libraries in the aforementioned code development process. FlexiBLAS allows switching between available BLAS implementations via user hosted configuration files or environment variables at runtime. Thus the need for time consuming recompilations is completely avoided and a smooth workflow is guaranteed. Moreover simple profiling, like counting BLAS calls and timing their execution, comes at basically no user effort and with minimal runtime overhead.

Inge Gutheil (Forschungszentrum Jülich), Johannes Grotendorst (Forschungszentrum Jülich)

For new exascale computers the degree of parallelism will increase leading to architectures with more and more cores per node and multithreading and SIMD on each core. Due to the limited memory per node the pure MPI parallelizing strategy may be no longer suited for these new architectures.

For applications that already use a hybrid parallelization due to memory limits, a pure MPI implementation will waste several cores for the eigensolver unless a hybrid version can be used.

Newly developed libraries for dense linear algebra computation such as the ELPA library for the computation of a large part of the eigenspectrum of dense real symmetric and complex Hermitan matrices thus offer versions with hybrid MPI and OpenMP parallelization. On Intel architectures it has been shown that this hybrid parallelization can perform better than the pure MPI parallel version. On BlueGene/Q the hybrid version is still experimental, but it shows promising results although the pure MPI version using all cores still performs better.

In this talk we will present performance evaluation results of the two eigensolvers of the library ELPA on BlueGene/Q and compare it to the results on the Intel Nehalem cluster JUROPA.

Thomas Auckenthaler (TU München), Thomas Huckle (TU München), Roland Wittmann (TU München)

In this talk, we present a new stable algorithm for the parallel QR-decomposition of "tall and skinny" matrices. The algorithm has been developed for the dense symmetric eigensolver ELPA, whereat the QR-decomposition of tall and skinny matrices represents an important substep. Our new approach is based on the fast but unstable CholeskyQR algorithm. We show the stability of our new algorithm and provide promising results of our MPI-based implementation on different clusters.

Edoardo Di Napoli (Forschungszentrum Jülich), Mario Berljafa (University of Manchester)

Sequences of eigenvalue problems consistently appear in a large class of applications based on the iterative solution of a non-linear eigenvalue problem. A typical example is given by the chemistry and materials science ab initio simulations relying on computational methods developed within the framework of Density Functional Theory (DFT). DFT provides the means to solve a high-dimensional quantum mechanical problem by representing it as a non-linear generalized eigenvalue problem which is solved self-consistently through a series of successive outer-iteration cycles. As a consequence each self-consistent simulation is made of several sequences of generalized eigenproblems . Each sequence, , groups together eigenproblems with increasing outer-iteration index .

In more general terms a set of eigenvalue problems is said to be a “sequence” if the solution of the -th eigenproblem determines, in an application-specific manner, the initialization of the -th eigenproblem. For instance at the beginning of each DFT cycle an initial function determines the initialization of the -th eigenproblem. A large fraction of eigenpairs are then use to compute a new which, in turn, leads to the initialization of a new eigenvalue problem . In addition to be part of a sequence, successive eigenproblems might possess a certain degree of correlation connecting their respective eigenpairs. In DFT sequences, correlation becomes manifest in the way eigenvectors of successive eigenproblems become progressively more collinear to each other as the -index increases.

We developed a subspace iteration method (ChFSI) specifically tailored for sequences of eigenproblems whose correlation appears in the form of increasingly collinear eigenvectors. Our strategy is to take the maximal advantage possible from the information that the solution of the eigenproblem is providing when solving for the successive problem. As a consequence the subspace iteration was augmented with a Chebyshev polynomial filter whose degree gets dynamically optimized so as to minimize the number of matrix-vector multiplications. The effectiveness of the Chebyshev filter is substantially increased when inputed the approximate eigenvectors , as well as very reliable estimates, namely , for the limits of the eigenspectrum interval to be filtered in. In addition the degree of the polynomial filter is adjusted so as to be minimal with respect to the required tolerance for the eigenpairs residual. This result is achieved by exploiting the dependence each eigenpair residual have with respect to its convergence ratio as determined by the rescaled Chebyshev polynomial and its degree. The solver is complemented with an efficient mechanism which locks and deflates the converged eigenpairs.

The resulting eigensolver was implemented in C language and parallelized for both shared and distributed architectures. Numerical tests show that, when the eigensolver is inputed approximate solutions instead of random vectors, it achieves up to a 5X speedup. Moreover ChFSI takes great advantage of computational resources by scaling over a large range of cores commensurate with the size of the eigenproblems. Specifically, by making better use of massively parallel architectures, the distributed memory version will allow DFT users to simulate physical systems quite larger than are currently accessible.

Andreas Stathopoulos (College of William and Mary), Lingfei Wu (College of William and Mary)

The Singular Value Decomposition (SVD) is a ubiquitous computational kernel in science and engineering. Many applications demand a few of the largest or smallest singular values of a large sparse matrix and the associated left and right singular vectors. The computation of the smallest singular triplets in particular presents challenges both to the speed of convergence and the accuracy that iterative methods can accomplish. We focus mainly on this difficult problem.

While traditional solutions were based on symmetric eigenvalue methods, over the last two decades the bidiagonalization Lanczos method has gained a much wider acceptance. Although in its original form it suffered from orthogonality loss, large memory demands, and problems with extracting the singular vectors through classical Rayleigh Ritz, recent techniques for restarting, refined and harmonic projections, and blocking have significantly increased its robustness. Despite this remarkable progress, current Lanczos bidiagonalization methods are still in development. Existing MATLAB implementations serve mainly as a testbed for mathematical research, and are not close to a general, high performance code. More importantly, preconditioning techniques cannot be used directly.

The JDSVD is the only method currently that can exploit preconditioning directly for the SVD problem. It employs several of the above techniques for restarting and projection, but this also is available only as a MATLAB code. Moreover, it relies on the solution of the correction equation to achieve good convergence, but the system it solves is maximally indefinite and thus difficult to solve. On the other hand, if no inner iterations are performed and no preconditioning is used, the outer iteration of JDSVD is much slower than bidiagonalization Lanczos.

We argue that a well designed, state-of-the-art eigensolver such as PRIMME, can offer faster speed of convergence than both of the above methods, even without preconditioning, while also exploiting the highly optimized implementation to achieve superior performance on sequential and parallel architectures. To accomplish this, we have engineered a hybrid, two-stage algorithm that addresses performance and accuracy at two successive stages.

One of the stumbling blocks of bidiagonalization Lanczos and JDSVD methods have been the irregular convergence to the smallest singular values. The problem is accentuated with restarting since the relevant singular vector directions may be discarded and have to be built again. A careful combination of harmonic and refined projections alleviate this problem, but their implementation can be very expensive. In contrast, a symmetric eigenvalue solver applied to the normal equations converges monotonically toward the smallest eigenvalues. Moreover, efficient restarting schemes exist that can achieve near-optimal convergence. PRIMME provides such methods that equal the convergence of the unrestarted bidiagonalization Lanczos, for a fraction of memory and time. To extract the eigenpairs from the normal equations, however, a squared projection must be used which limits the maximum achievable accuracy for the smallest eigenpairs. The idea is that after most of the digits of accuracy have been obtained through a very fast method on the normal equations, we resolve the rest of the accuracy by switching to an interior eigenvalue problem on the augmented matrix using a highly tuned, Jacobi-Davidson method. The availability of excellent shifts for the correction equation acts as an iterative refinement with the additional benefit of subspace acceleration.

We present theoretical and experimental evidence for the merit of this approach, and emphasize that a highly optimized, preconditioned, parallel code for the SVD problem now exists.

Martin Köhler (Max Planck Institute for Dynamics of Complex Technical Systems), Peter Benner (Max Planck Institute for Dynamics of Complex Technical Systems), Jens Saak (Max Planck Institute for Dynamics of Complex Technical Systems)

The solutions of (generalized) Lyapunov equations play a key role in many applications in systems and control theory. Their robust numerical computation, when the full solution is sought, is considered solved since the seminal work of Bartels and Stewart. The algorithm consists mainly of two phases: a (quasi-)triangular reduction and a forward/backward substitution. A number of variants of what is now called the Bartels-Stewart algorithm have been proposed. Here, we focus on the second phase which consists mainly of BLAS level-2 directives and which is used in standard implementations, as e.g. the solver lyap in the MATLAB^® Control Systems Toolbox (which is based upon the routine SG03AY from the systems and control library SLICOT). On modern computers, however, the formulation of BLAS level-3 type implementations is crucial to enable optimal usage of cache hierarchies and modern block scheduling methods, which are based on directed acyclic graphs describing the interdependence of single block computations.

Our contribution transforms the available level-2 variants for the solution of the (quasi-)triangular generalized Lyapunov equation to level-3 versions. Thus, by comparison on a standard multicore machine and selection of the most appropriate variant we can close the aforementioned gap in the implementations.

To this end, we analyze the implementation details of the current SLICOT version and identify the crucial operations which prevent a blocked implementation. We consider different algorithms to handle the arising generalized Sylvester equations and compare the flop counts and their performance on current computer architectures. Additionally, we focus on implementation details that avoid unnecessary computations and guarantee the symmetry of the solution of the Lyapunov equation. In a number of numerical experiments we point out which of the covered variants should replace the SLICOT implementation of the solver for (quasi-)triangular generalized Lyapunov equations with a comparable accuracy and much improved performance. Significant speed-ups over current standard solvers can be observed already for moderately sized Lyapunov equations.

Wilfried Gansterer (University of Vienna), Karl Prikopa (University of Vienna)

We present a novel truly distributed linear least squares solver PSDLS-IR for a rectangular matrix distributed row-wise. Upon termination, each node holds an approximation of the full solution vector. Since we are interested in a truly distributed algorithm, we do not employ a fusion centre or multi hop communication and limit the communication of each node to its immediate neighbourhood.

Mathematically, PSDLS-IR is based on the semi-normal equations (SNE). It is well known that SNE is in general not backward stable. In combination with iterative refinement the method can be stabilised. We consider mixed precision iterative refinement, which increases the performance of the algorithm because the computationally and communication intensive operations are reduced. This leads to a significant speed-up compared to existing truly distributed linear least squares solvers.

In PSDLS-IR, all communication between nodes takes place in the context of reduction operations. They can, for example, be performed using randomized gossip algorithms where each node only communicates with its immediate neighbourhood. Moreover, fault tolerant reduction algorithms become applicable. We focus on the push-sum algorithm and the fault tolerant push-flow algorithm. We analyse the communication cost of PSDLS-IR and compare it to existing distributed linear least squares solvers. Moreover, we summarize simulation experiments showing that our PSDLS-IR solver requires significantly fewer messages per node to reach a high accuracy solution than the existing methods.

Jonathan Hogg (Rutherford Appleton Laboratory), Evgueni Ovtchinnikov (Rutherford Appleton Laboratory), Jennifer Scott (Rutherford Appleton Laboratory)

The sparse LDL^T factorization of symmetric matrices is a key kernel in optimization methods. These methods often produce numerically challenging sparse indefinite matrices that require pivoting to obtain an accurate solution and estimate of the inertia.

We present a new GPU implementation of this kernel that includes full pivoting and performs all numerical operations on the GPU. This avoids the PCI-express bottleneck that has limited the effectiveness of such solvers in the past.

We also describe the use of some novel numerical techniques such as the use of a posteriori pivoting to minimize the overhead of pivot acceptance tests, and the use of a stable partial inverse of L to accelerate the solve phase.

Ioannis E. Venetis (University of Patras), Alexandros Kouris, Alexandros Sobczyk, Efstratios Gallopoulos (University of Patras), Ahmed Sameh (Purdue University)

Solving banded linear systems is an important computational kernel in many applications. Over the years, a variety of algorithms and modules in numerical libraries for this task have been proposed that are suitable for computer systems ranging from uniprocessors to multicores and high performance coprocessors. In this presentation we focus on tridiagonal systems as these are very common and are also encountered as kernels in several linear algebra computations. Interesting recent efforts focus on GPU implementations. Algorithms based on cyclic reduction, parallel cyclic reduction and recursive doubling and hybrids were studied in [4]; these depend on symmetric positive definiteness or diagonal dominance to proceed without breakdown.

For very large tridiagonal and narrow banded matrices, a well known approach for obtaining parallelism and high performance is to use partitioning. As shown recently by Chang et al., the Spike partitioning framework ([2, 3]) combined with a special pivoting algorithm for the resulting subsystems can lead to a scalable, highly performing tridiagonal solver [1].

In this presentation we revisit a Spike approach combined with Givens rotations to solve the subsystems, as first proposed in [3], present its CUDA implementation and evaluate its performance and stability.

[1] L.-W. Chang, J.A. Stratton, H.S. Kim, and W.-M.W. Hwu. A scalable, numerically stable, high-performance tridiagonal solver using GPUs. In Proc. Int'l. Conf. High Performance Computing, Networking Storage and Analysis, SC '12, pages 27:1-27:11, Los Alamitos, CA, USA, 2012. IEEE Computer Society Press. [2] E. Polizzi and A.H. Sameh. A parallel hybrid banded system solver: The SPIKE algorithm. Parallel Computing, 32(2):177-194, 2006. [3] A.H. Sameh and D.J. Kuck. On stable parallel linear system solvers. J. Assoc. Comput. Mach., 25(1):81-91, January 1978. [4] Y Zhang, Jonathan Cohen, and John D Owens. Fast tridiagonal solvers on the GPU. ACM Sigplan Notices, 45(5):127-136, 2010.

Bruno Carpentieri (University of Groningen), Jia Liao (University of Groningen), Masha Sosonkina (Old Dominion University ), Aldo Bonfiglioli (University of Basilicata)

The starting point for this project was the iterative solution of sparse and block structured linear systems that arise from the analysis of turbulent flows in computational fluid dynamics applications. In the last two years we have studied preconditioning techniques based on block multilevel incomplete LU factorization preconditioners for this problem class, and we have found them to be quite effective in reducing the number of GMRES iterations. These preconditioners exhibit better parallelism than standard ILU algorithms due to the recursive factorization and they may be noticeably more robust for comparable memory usage especially for solving large problems. Additionally, exploiting the available block structure in the matrix can maximize computational efficiency.

Sparse matrices arising from the solution of systems of partial differential equations often exhibit a perfect block structure, meaning that the nonzero blocks in the sparsity pattern are fully dense (and typically small), e.g., when several unknown quantities are associated with the same grid point. However, similar block orderings can be sometimes unravelled also on general unstructured matrices, by ordering consecutively rows and columns with a similar sparsity pattern, and treating some zero entries of the reordered matrix as nonzero elements, and the nonzero blocks as dense, with a little sacrifice of memory. The reordering results in linear systems with blocks of variable size in general.

Our recently developed parallel package pVBARMS (parallel variable block algebraic recursive multilevel solver) for distributed memory computers takes advantage of these frequently occurring structures in the design of the multilevel incomplete LU factorization preconditioner, and maximizes computational efficiency achieving increased throughput during the computation and improved reliability on realistic applications. The method detects automatically any existing block structure in the matrix without any user’s prior knowledge of the underlying problem, and exploits it to maximize computational efficiency.

We describe a novel parallel MPI-based implementation of pVBARMS for distributed memory computers based on the block Jacobi, the additive Schwarz and the Schur-complement methods. We address hybrid MPI and OpenMP implementations and gain further parallelism using Many-Integrated Codes (MIC) technology. Implementation details are always critical aspects to consider in the design of sparse matrix algorithms. Therefore, we revisit our original implementation of the partial (block) factorization step with a careful selection of the parameters, and we compare different algorithms for computing the block ordering based on either graph or matrix analysis.

Finally, we report on the numerical and parallel scalability of the pVBARMS package for solving turbulent Navier-Stokes equations on a suite of two- and three-dimensional test cases, among which the calculation of the flow past the DPW3-W1 wing configuration of the third AIAA Drag Prediction Workshop, which is the application that motivated this study. In our formulation the mean flow and turbulent transport equations are solved in coupled form using a Newton-Krylov algorithm. These analyses are carried with coarse- to medium-sized grids featuring up to 6.5 million nodes at Reynolds number equal to .

Ziad Sultan (University of Grenoble), Jean-Guillaume Dumas (University of Grenoble), Clément Pernet (University of Grenoble), Thierry Gautier (INRIA)

We propose efficient parallel algorithms and implementations of LU factorization over word size prime fields and target shared memory architectures. With respect to numerical linear algebra, we have identified three major differences, specific to exact linear algebra.

First, the arithmetic complexity could be dominated by modular reductions. Therefore, it is mandatory to delay as much as possible these reductions. Delaying the latter can be done by accumulating several multiplications before reducing while keeping the result exact. This induces splitting the matrices in blocks of large size to accumulate multiplication operations as much as possible, without overflow. Then one performs only one reduction operation for several of those accumulations. Simultaneously, to gain in parallelism, the splitting should mix fine-grain parallelizations of tiled iterative and recursive algorithms. A finer grain allows more flexibility in the scheduling when running numerous cores: different independent blocks can be treated at the same time. This trade-off also creates challenges on the efficiency of the underlying scheduler and on the minimization of distant memory accesses in order to reduce the dependency on the bus speed.

Second, in numerical linear algebra, the cost of arithmetic operations is more or less associative: with dimensions above a rather low threshold, sequential matrix multiplication of the BLAS reaches the peak efficiency of the processor. Hence the granularity has very little impact on the efficiency of a block algorithm run in sequential. Over a finite field numerical stability is not an issue, and asymptotically fast matrix multiplication algorithms, like Winograd's variant of Strassen's fast algorithm can be used on top of the BLAS. Their speed-up increases with matrix dimension. The cost of sequential matrix multiplication over finite field is therefore not associative: a larger granularity delivers better sequential efficiency. Thus, trade-offs are also to be made between size of blocks well suited to those fast variants or to load and communication balancing.

Third, many applications over finite fields require the rank profile of the matrix rather than the solution to a linear system. It is thus important to design parallel algorithms that preserve and compute this rank profile. This profile is a key invariant used in many applications using exact Gaussian elimination, such as Grobner basis computations and computational number theory. Now, in dense numerical linear algebra, pivoting is used to ensure good numerical stability and good data locality. In the context of exact linear algebra, only certain pivoting strategies will reveal the echelon form or, equivalently, the rank profile of the matrix. Moreover, the rank deficiency of some blocks also leads to unpredictable dimensions of the tiles or slabs as the rank profile is only discovered during the algorithm. Thus the computing load for each task is not known in advance (some panel blocks may have high rank deficiency), making the tasks very heterogeneous. The block splitting has then to be dynamically computed.

We thus propose and compare different parallelization strategies and block decomposition to tackle these issues: tile iterative with left-looking, right-looking and Crout variants with delayed modular reductions; slab and tile recursive variants; dynamic/static block cutting and adaptive coupling of the different variants. Experiments demonstrate that the tile recursive variant performs better and matches the performance of reference numerical software when no rank deficiency occur. Furthermore, even in the most heterogeneous case, namely when all pivot blocks are rank deficient, we show that it is possible to maintain a high efficiency.

Jiri Stary (Academy of Sciences of the Czech Republic), Radim Blaheta (Academy of Sciences of the Czech Republic), Ondrej Jakl (Technical University of Ostrava)

We deal with the solution of high resolution FEM systems arising from the analysis of the macroscale response of materials through investigation of processes in their microstructure. By other words, we test the smaller material samples with relatively complicated inner structure and try to determine effective (homogenized) material properties. In computational micromechanics, the testing of such samples means solution of boundary value problems on test domains involving the microstructure with loading provided by suitable boundary conditions.

We focus on X-ray CT image based micromechanics of geomaterials with the use of continuum mechanics and the finite element computation of the microscale strains and stresses. This means that basic information about the microstructure is provided by analysing (segmentation) of 3D images of real samples. This information should be completed by information on local material properties, i.e. material properties of the individual material constituents.

There is a strong need for high performance parallel computing at several stages of the computational micromechanics, namely at - analysis of CT scans, - high resolution finite element solution of boundary value problems, - solution of inverse problems for determination or calibration of local material properties.

In this contribution, we focus on the second point, i.e. solving the high resolution finite element systems with tens or hundreds degrees of freedom on available parallel computers at the Institute of Geonics and the IT4Innovations supercomputer centre in Ostrava. We describe efficiency of the in-house GEM solvers exploiting the Schwarz domain decomposition method with aggregation by performing computational experiments on the above parallel computers. The solution of these systems is also necessary for building efficient solution methods for inverse material identification problems.

Nam Zin Cho (KAIST), Han Jong Yoo (KAIST)

Nuclear reactor core analysis via neutron transport method in a three-dimensional whole-core domain is one of the challenges facing reactor physicists. Monte Carlo method is a stochastic method which is directly simulating particles with minimum assumptions. Method of characteristics (MOC) is a deterministic approach which is capable of analyzing heterogeneous geometry with energy groups. Both of these methods are computationally very expensive when applied directly to 3-D whole-core analysis. A typical nuclear reactor of current generation is radially very heterogeneous but is usually piecewise homogeneous in axial direction. As appropriate for this type of reactors, a 2D/1D fusion method was developed in the previous work for 3-D whole-core analysis, which employs two-dimensional method of characteristics in radial direction and one-dimensional SN-like method in axial direction. Both 2-D and 1-D problems are coupled by axial and radial leakage terms. With a small mesh size, however, these leakage terms may cause negative source term and ruin the convergence of the method. In this paper, a novel scheme of 3D/2D rotational plane slicing method is proposed for whole-core transport calculation. A 3-D problem is sliced along the characteristics line, forming multiple 2-D planes (characteristics planes). These 2-D planes usually consist of rectangular mesh cells in reactor core problems. Various 2-D SN solvers are available for these 2-D transport problems and the simplest diamond-difference (DD) scheme is chosen for this work. Parallel computing of this approach is efficient since each characteristic plane is independent from each other. The method is tested on several benchmark problems with encouraging results for accuracy and parallelization.

Remigijus Paulavičius (Vilnius University), Julius Žilinskas (Vilnius University)

The availability of computer power offers potential and challenges for solving difficult scientific optimization problems using global optimization algorithms. Many of these algorithms have large memory requirements and strong data dependency that degrade the program scalability and efficiency of parallel algorithms. Global optimization algorithms DIRECT (DIviding RECTangles) and DISIMPL (DIviding SIMPLices) are such kind of algorithms. This talk is concerned with a parallel implementation and the application of parallel algorithms to applied problems arising in civil engineering. Systems of linear equations often appear in computation of objective function and constraints of such problems. In our results we investigate the efficiency of parallel optimization algorithms for these problems.

Alvaro Frank (RWTH Aachen), Diego Fabregat-Traver (RWTH Aachen), Paolo Bientinesi (RWTH Aachen)

We present high-performance algorithms and implementations to carry out Genome Wide Association Studies (GWAS) of omics data on multi-threaded architectures. GWAS study the relation between phenotype traits and genetic variations of individuals, while the term omics typically implies large datasets and large-scale analyses. We consider a specific class of GWAS that rely on the solution of billions of ordinary least squares (OLS) problems , where determines the relation between a genetic marker and a trait . In GWAS, the above equation has to be solved for the cross product of matrices and vectors of size , resulting in the computation of the order of peta-floating point operations (petaflops) over terabytes of data. The proposed solvers deal with the high computational requirements of GWAS efficiently, reducing considerably the required time to solution when compared to existing solvers.

Current solvers, based on either traditional libraries (e.g., LAPACK) or domain-specific ones (e.g., ProbABEL from the GenABEL suite), suffer from a number of limitations. The main one comes from the fact that these libraries solve individual OLS problems in isolation, which hinders their efficiency and scalability. Additionally, these solvers stream data from/to disk synchronously which, given the size of the involved datasets, adds a noticeable overhead. Our solvers overcome these shortcomings and are thus able to treat large problems in reduced time frames. First, we exploit application-specific knowledge, such as the correlation among problems and the properties of the matrices, to remove redundant computations and lower the computational cost of the entire analysis.

Then, we model and rearrange the required data transfers, and make use of out-of-core techniques so that input/output (I/O) is, when possible, completely overlapped with computation.

Finally, we merge the solution of multiple OLS problems to cast the bulk of the computation in terms of efficient BLAS-3 kernels, and combine the use of multi-threaded BLAS with OpenMP to attain scalability.

As the experiments show, our solvers attain a high efficiency (from 70% to 85% of the peak performance), and when executed on a 40-core machine, they attain a speedup of about 32x. Also, I/O overhead is successfully overlapped with computations. As a result, our routines outperform current state-of-the-art tools by a factor of 100.

Gemma Sanjuan Gomez (Autonomous University of Barcelona), Tomàs Margalef (Autonomous University of Barcelona), Ana Cortés (Autonomous University of Barcelona)

Forests fires are a significant problem especially in countries of the Mediterranean basin. To fight against these disasters, the accurate prediction of forest fire propagation is a crucial issue. Propagation models try to describe the future evolution of the forest fire given an initial scenario. However, the data describing the real fire scenario are usually subject to high levels of uncertainty. Moreover, there are input parameters that present spatial and temporal variation that make the prediction less accurate.

Wind speed and direction are parameters that affect forest fire propagation dramatically. So, an accurate estimation of such parameters is crucial to predict the fire propagation precisely. WindNinja is a wind field simulator that can easily be coupled to a forest fire propagation simulator such as FARSITE. However, wind field simulators present to main drawbacks when the map describing the terrain is large ( Km): On the one hand the execution time is prohibitive and, on the other hand, the amount of memory required to tackle such problem is usually not available in a single node. Both problems are related to the number of cells describing the terrain.

To evaluate the wind field WindNinja apply mass conservation equations on the terrain surface and the mesh describing the terrain grows with the map size. However, the matrix describing the physical equations is a sparse matrix that is represented using CSR representation method. WindNinja solves the system applying a conjugate gradient method with LU preconditioner that converges quite fast, but scale very badly since it cannot be parallelized. To overcome this problem and reduce execution time and memory requirements a Domain Decomposition method with non-overlapping, commonly called Schur method has been applied.

The Schur method is a graph partitioning mathematical method that can be used to to divide the mesh representing the terrain map. Partitioning the mesh implies a partitioning of the equation system. So, the matrix describing the system of equations is arranged in such a way that there are submatrices representing the local interactions within any mesh partition and another one describing the interactions in the borders of the partitions. In this way the systems to be solved are smaller and can be solved in parallel.

The Schur method has been implemented using the METIS library to partition the CSR matrix describing the whole system. This library includes load balancing facilities that automatically distribute the non-zero elements of the matrix among the different partitions. The resulting system is solved applying the PARDISO library. This library different solvers combining direct and iterative methods and it presents very good scalability on SMPs and clusters of SMPs. Taking advantatge of PARDISO features different solvers have been applied including direct methods and iterative methods such as conjugate gradient with and without preconditioner.

A significant point to be considered is that partitioning the mesh in more parts to exploit more parallelism implies an increase in the points in the borders and, therefore, the submatrix representing the interactions in the borders includes more terms and is larger. Larger matrices require more complexity and more number of iterations to converge. So, it is necessary to determine the optimal number of partition that exploits parallelism, but does not penalizes execution time due to the increase in the matrix representing the borders among the different partitions of the mesh. Therefore, a partitioning methodology has been developed to determine the best partitioning.

The results show that the application of the Schur method with the proposed partitioning reduces execution time significantly and offers a good scalability.