IBM SP2 Research Articles

A parallel Monte Carlo code for planar and SPECT imaging: implementation, verification and applications in 131I SPECT

This paper reports the implementation of the SIMIND Monte Carlo code on an IBM SP2 distributed memory parallel computer. Basic aspects of running Monte Carlo particle transport calculations on parallel architectures are described. Our parallelization is based on equally partitioning photons among the processors and uses the Message Passing Interface (MPI) library for interprocessor communication and the Scalable Parallel Random Number Generator (SPRNG) to generate uncorrelated random number streams. These parallelization techniques are also applicable to other distributed memory architectures. A linear increase in computing speed with the number of processors is demonstrated for up to 32 processors. This speed-up is especially significant in Single Photon Emission Computed Tomography (SPECT) simulations involving higher energy photon emitters, where explicit modeling of the phantom and collimator is required. For 131I, the accuracy of the parallel code is demonstrated by comparing simulated and experimental SPECT images from a heart/thorax phantom. Clinically realistic SPECT simulations using the voxel-man phantom are carried out to assess scatter and attenuation correction.

Scheduling GEN_BLOCK Array Redistribution

This article is devoted to the redistribution of arrays that are distributed in a GEN_BLOCK fashion over a processor grid. While GEN_BLOCK redistribution is essential for load balancing, prior research about redistribution has been focused on block-cyclic redistribution. The proposed scheduling algorithm exploits a spatial locality in message passing from a seemingly irregular array redistribution. The algorithm attempts to obtain near optimal scheduling by trying to minimize communication step size and the number of steps. According to experiments on CRAY T3E and IBM SP2, the algorithm shows good performance in typical distributed memory machines.

Architecture independent parallel algorithm design: theory vs practice

We propose architecture independent parallel algorithm design as a framework for writing parallel code that is scalable, portable and reusable. Towards this end we study the performance of some dense matrix computations such as matrix multiplication, LU decomposition and matrix inversion. Although optimized algorithms for these problems have been extensively examined before, a systematic study of an architecture independent design and analysis of parallel algorithms and their performance (including matrix computations) has not been undertaken. Even though more refined algorithms and implementations (sequential or parallel) for the stated problems exist, the complexity and performance of the introduced algorithms is sufficient to raise the issues that are important in architecture independent parallel algorithm design. Two established distributions of an input matrix among the processors of a parallel machine are examined and the particular theoretical and practical merits of each one are also discussed. The algorithms we propose have been implemented and tested on a variety of parallel systems that include the SGI Power Challenge, the IBM SP2 and the Cray T3D. Our experimental results support our claims of efficiency, portability and reusability of the presented algorithms.

Parallel approaches for singular value decomposition as applied to robotic manipulator Jacobians

The system of equations that govern kinematically redundant robotic manipulators is commonly solved by finding the singular value decomposition (SVD) of the corresponding Jacobian matrix. This can require a considerable amount of time to compute, thus a parallel SVD algorithm reducing execution time is sought. The approach employed here lends itself to parallelization by using Givens rotations and information from previous decompositions. The key contribution of this research is the presentation and implementation of parallel SVD algorithms to compute the SVD for a set of Jacobians that represent various different joint failure scenarios. Results from implementation of the algorithm on a MasPar MP-1, an IBM SP2, and the PASM prototype parallel computers are compared. Specific issues considered for each implementation include: how data is mapped to the processing elements, the effect that increasing the number of processing elements has on execution time, the type of parallel architecture used, and trade-offs between modes of parallelism.

Development and Application of a Parallel Implicit Solver for Unsteady Viscous Flows

This work investigates the performance and application of a parallel version of a three-dimensional second-order time accurate Navier-Stokes solver based on an implicit approximate-factorization Beam-Warming algorithm. A systematic incremental approach for parallelizing the serial code was developed which ensures that the parallel version of the code produces identical results to the original serial code. The current parallel scheme decomposes the grid using two-dimensional multipartitioning to evenly distribute the work across multiple processors with parallel communication via Message-Passing Interface (MPI) library. The code's performance has been assessed on three supercomputers: the IBM SP2, IBM SP3 and the Silicon Graphics Origin 2000. The solver is validated for Couette flow, and both steady and unsteady flow over a circular cylinder. Additional applications include both two- and three-dimensional flow over a stationary and a rotationally oscillating circular cylinder. This new solver enables the efficient simulation of large-scale unsteady viscous flows employing grids containing on the order of 10 7 points using available parallel supercomputers.

A Parallel Implementation of the Nonsymmetric QR Algorithm for Distributed Memory Architectures

One approach to solving the nonsymmetric eigenvalue problem in parallel is to parallelize the QR algorithm. Not long ago, this was widely considered to be a hopeless task. Recent efforts have led to significant advances, although the methods proposed up to now have suffered from scalability problems. This paper discusses an approach to parallelizing the QR algorithm that greatly improves scalability. A theoretical analysis indicates that the algorithm is ultimately not scalable, but the nonscalability does not become evident until the matrix dimension is enormous. Experiments on the Intel Paragon system, the IBM SP2 supercomputer, the SGI Origin 2000, and the Intel ASCI Option Red supercomputer are reported.

Fast pseudorandom-number generators with modulus 2kor 2k-1 using fused multiply-add

Many numerically intensive computations done in a scientific computing environment require uniformly distributed pseudorandom numbers in the range (0, 1) and (−1, 1). For multiplicative congruential generators with modulus 2 <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">k</sup> , k ≤ 52, and period 2 <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">k-2</sup> , we show that the cost per random number for these two distributions is 3 and 3.125 multiply–adds on RS/6000® processors. Our code, on the IBM POWER2 Model 590, produces more than 40 million uniformly distributed pseudorandom numbers per second for both ranges (0, 1) and (−1, 1). Additionally, our code sustains the 40 million per second rate for data out of cache. The Numerical Aerodynamic Simulation (NAS) parallel benchmarks use a linear congruential generator with modulus 2 <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">46</sup> . Our result is about 50 times faster than the generic implementation given in the benchmarks. The extra-accuracy fused multiply-add instruction of RS/6000 machines combined with a few algorithmic innovations gives rise to the 50-fold increase. If IEEE 64-bit arithmetic is used with our Fortran code on POWER and PowerPC® architectures, the results we obtain are bit-wise identical to the generic algorithms. The paper gives several illustrations of a general technique called the Algorithm and Architecture approach. We demonstrate herein that programmer-controlled unrolling of loops is equivalent to “customized vectorization of RISC-type code.” Customized vectorization is more powerful than ordinary vectorization, and it is only possible on RISC-type machines. We illustrate its use to show that RS/6000 processors can compute the distribution (−1, 1) at the rate of 3.125 multiply–adds. We also specify a linear congruential generator that is related to the multiplicative congruential generator referred to above. It has a full period of 2 <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">k</sup> , where 2 <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">k</sup> is the modulus. The cost per random number [in the range (0, 1)] for this generator is four multiply–adds on RS/6000 processors. Our code, on the IBM POWER2 Model 590, for this generator produces more than 30 million uniformly distributed pseudorandom numbers per second for the range (0, 1). We show that this generator is “embarrassingly parallel,” or EP. Using the Algorithm and Architecture approach, we describe a new concept called “generalized unrolling.” Finally, we present a multiplicative congruential generator for which the modulus is not a power of 2. Such a generator, as well as one with modulus 2 <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">k</sup> , is selectable as the generator used in the RANDOM_NUMBER intrinsic function of IBM XL Fortran and XL High Performance Fortran. All of the generators reported here are EP. Using an IBM SP2 machine with 250 wide nodes, it is possible to compute more than ten billion uniform random numbers in a second.

Parallel Algorithms for Particles–Turbulence Two-Way Interaction Direct Numerical Simulation

Understanding the demixing effect on the dispersion of particles by large-scale turbulence is very important in practical applications. Using pseudo-spectral and Lagrangian approaches, we have simulated a three-dimensional particle-laden mixing layer under one-way coupling effect. However, the computer resource required to simulate such a two-phase flow with high Reynolds number and two-way momentum coupling effect exceeds the limit of the current single processor. In this paper, the computation of particles and the two-way momentum coupling terms are partitioned in the span-wise direction because particles are distributed most evenly in this direction. The computation of the tree-dimensional flow field is first partitioned into three groups of processors because of the most independence of the computation among the three spatial dimensions. In each group, the domain is then partitioned using two different schemes based on the property of the fast Fourier transformation. The first one, the master–slave scheme, is employed for Algorithm MS due to its simplicity and overlapping of communication and computation. The second one, the transpose approach, is used for Algorithm TP order to partition all of the flow field computation. An analysis shows that compared to Algorithm MS, Algorithm TP can also reduce nearly a half of the amount of communication work. Experiments show that Algorithm MS has obtained a speedup of 4.3 using 9HP workstations and a speedup of 6.4 using 15 nodes of IBM SP-2 for a problem size on the order of 643, and Algorithm TP has achieved speedups 44% higher than Algorithm MS.

An I/O-Conscious Tiling Strategy for Disk-Resident Data Sets

This paper describes a tiling technique that can be used by application programmers and optimizing compilers to obtain I/O-efficient versions of regular scientific loop nests. Due to the particular characteristics of I/O operations, a straightforward extension of the traditional tiling method to I/O-intensive programs may result in poor I/O performance. Therefore, the technique presented in this paper adapts iteration space tiling for I/O-performing loop nests to deliver high I/O performance. The generated code results in huge savings in the number of I/O calls as well as the volume of data transferred between the disk subsystem and main memory. Our experimental results on the IBM SP-2 distributed-memory message-passing multiprocessor demonstrate that the reduction in these two parameters, namely, the number of I/O calls and the transferred data volume, can lead to a marked decrease in overall execution times of I/O-intensive loop nests. In a number of loop nests extracted from several benchmarks and math libraries, we were able to improve the execution times by an average 42.5% for one data set and by an average 47.4% for another.

Parallel Shear-Warp Factorization Volume Rendering Using Efficient 1-D and 2-D Partitioning Schemes for Distributed Memory Multicomputers

3-D data visualization is very useful for medical imaging and computational fluid dynamics. Volume rendering can be used to exhibit the shape and volumetric properties of 3-D objects. However, volume rendering requires a considerable amount of time to process the large volume of data. To deliver the necessary rendering rates, parallel hardware architectures such as distributed memory multicomputers offer viable solutions. The challenge is to design efficient parallel algorithms that utilize the hardware parallelism effectively. In this paper, we present two efficient parallel volume rendering algorithms, the 1D-partition and 2D-partition methods, based on the shear-warp factorization for distributed memory multicomputers. The 1D-partition method has a performance bound on the size of the volume data. If the number of processors is less than a threshold, the 1D-partition method can deliver a good rendering rate. If the number of processors is over a threshold, the 2D-partition method can be used. To evaluate the performance of these two algorithms, we implemented the proposed methods along with the slice data partitioning, volume data partitioning, and sheared volume data partitioning methods on an IBM SP2 parallel machine. Six volume data sets were used as the test samples. The experimental results show that the proposed methods outperform other compatible algorithms for all test samples. When the number of processors is over a threshold, the experimental results also demonstrate that the 2D-partition method is better than the 1D-partition method.

Achieving Extreme Resolution in Numerical Cosmology Using Adaptive Mesh Refinement: Resolving Primordial Star Formation

As an entry for the 2001 Gordon Bell Award in the "special" category, we describe our 3-d, hybrid, adaptive mesh refinement (AMR) codeEnzodesigned for high-resolution, multiphysics, cosmological structure formation simulations. Our parallel implementation places no limit on the depth or complexity of the adaptive grid hierarchy, allowing us to achieve unprecedented spatial and temporal dynamic range. We report on a simulation of primordial star formation which develops over 8000 subgrids at 34 levels of refinement to achieve a local refinement of a factor of 1012in space and time. This allows us to resolve the properties of the first stars which form in the universe assuming standard physics and a standard cosmological model. Achievingextreme resolutionrequires the use of 128-bit extended precision arithmetic (EPA) to accurately specify the subgrid positions. We describe our EPA AMR implementation on the IBM SP2 Blue Horizon system at the San Diego Supercomputer Center.

Scalable Atomistic Simulation Algorithms for Materials Research

A suite of scalable atomistic simulation programs has been developed for materials research based on space‐time multiresolution algorithms. Design and analysis of parallel algorithms are presented for molecular dynamics (MD) simulations and quantum‐mechanical (QM) calculations based on the density functional theory. Performance tests have been carried out on 1,088‐processor Cray T3E and 1,280‐processor IBM SP3 computers. The linear‐scaling algorithms have enabled 6.44‐billion‐atom MD and 111,000‐atom QM calculations on 1,024 SP3 processors with parallel efficiency well over 90%. production‐quality programs also feature wavelet‐based computational‐space decomposition for adaptive load balancing, spacefilling‐curve‐based adaptive data compression with user‐defined error bound for scalable I/O, and octree‐based fast visibility culling for immersive and interactive visualization of massive simulation data.

Large-scale structural analysis by parallel multifrontal solver through Internet-based personal computers

Internet supercomputing methodology is introduced, and the concept is realized for large-scale finite element analysis. This is enabled by an efficient out-of-core parallel solver, which is based on domainwise or elementwise multifrontal approach. The primary resources of Internet supercomputing are numerous idling personal computers connected by the Internet irrespective of their locations. The computing ability of hundreds or thousands personal computers connected by the Internet can be as powerful as that of supercomputer such as the IBM SP system if these personal computers are utilized for solving a problem simultaneously through an efficient parallel computing algorithm. Under the described concept, a virtual supercomputing system InterSup I is constructed and tested. To establish the InterSup I system, 64 personal computer nodes, which are located in several places and connected by the Internet, are conscripted. When the established InterSup I system is used, linear static analyses of a finite element model having around 2 × 10 6 degrees of freedom are done through the developed parallel multifrontal solver. Some numerical tests have been carried out to investigate the affordability and effectiveness of Internet supercomputing. The possibility of Internet supercomputing through Internet personal computers as an alternative supercomputing power for high-performance computing is provided in this research.

A 3D incompressible Navier–Stokes velocity–vorticity weak form finite element algorithm

AbstractThe velocity–vorticity formulation is selected to develop a time‐accurate CFD finite element algorithm for the incompressible Navier–Stokes equations in three dimensions.The finite element implementation uses equal order trilinear finite elements on a non‐staggered hexahedral mesh. A second order vorticity kinematic boundary condition is derived for the no slip wall boundary condition which also enforces the incompressibility constraint. A biconjugate gradient stabilized (BiCGSTAB) sparse iterative solver is utilized to solve the fully coupled system of equations as a Newton algorithm. The solver yields an efficient parallel solution algorithm on distributed‐memory machines, such as the IBM SP2. Three dimensional laminar flow solutions for a square channel, a lid‐driven cavity, and a thermal cavity are established and compared with available benchmark solutions. Copyright © 2002 John Wiley &amp; Sons, Ltd.

New computer promises to accelerate climate change research

The acquisition of a new IBM SP supercomputer and advanced technologies, announced on November 9 by a Boulder, Colorado‐based research organization, could accelerate research into global and regional climate change, short‐ and long‐range weather prediction, atmospheric chemistry, space weather, and other areas.The equipment, announced by the National Center for Atmospheric Research (NCAR) and its parent organization, the University Corporation for Atmospheric Research (UCAR), will provide atmospheric scientists with access to new computational, storage, and communications capabilities. The U.S. National Science Foundation (NSF) purchased the machine for use at NCAR to advance geoscience research.

Optimizing noncontiguous accesses in MPI–IO

The I/O access patterns of many parallel applications consist of accesses to a large number of small, noncontiguous pieces of data. If an application's I/O needs are met by making many small, distinct I/O requests, however, the I/O performance degrades drastically. To avoid this problem, MPI–IO allows users to access noncontiguous data with a single I/O function call, unlike in Unix I/O. In this paper, we explain how critical this feature of MPI–IO is for high performance and how it enables implementations to perform optimizations. We first provide a classification of the different ways of expressing an application's I/O needs in MPI–IO – we classify them into four levels, called levels 0–3. We demonstrate that, for applications with noncontiguous access patterns, the I/O performance improves dramatically if users write their applications to make level-3 requests (noncontiguous, collective) rather than level-0 requests (Unix style). We then describe how our MPI–IO implementation, ROMIO, delivers high performance for noncontiguous requests. We explain in detail the two key optimizations ROMIO performs: data sieving for noncontiguous requests from one process and collective I/O for noncontiguous requests from multiple processes. We describe how we have implemented these optimizations portably on multiple machines and file systems, controlled their memory requirements, and also achieved high performance. We demonstrate the performance and portability with performance results for three applications – an astrophysics-application template (DIST3D), the NAS BTIO benchmark, and an unstructured code (UNSTRUC) – on five different parallel machines: HP Exemplar, IBM SP, Intel Paragon, NEC SX-4, and SGI Origin2000.

Design Alternatives for Virtual Interface Architecture and an Implementation on IBM Netfinity NT Cluster

The Virtual Interface Architecture (VIA) specification has been developed to standardize user-level network interfaces that provide low-latency, high-bandwidth communications. Few hardware and software implementations of VIA exist. Since the VIA specification is flexible, different choices exist for implementing various components of VIA such as doorbells, address translation methods, and completion queues. Although previous studies have evaluated the overall performance of different VIA implementations, there has not been a comparative study on the performance of VIA components. In this paper, we evaluate and compare the performance of different implementations of essential VIA components. We discuss the pros and cons of each design approach and describe the required support for implementing each of them. Then, we discuss an experimental implementation of the Virtual Interface Architecture for the IBM SP Switch-Connected NT cluster, one of the newest clustering platforms available. We discuss different design issues involved in this implementation. In particular, we explain how the virtual-to-physical address translation is implemented efficiently with a minimum Network Interface Card (NIC) memory requirement. We show how caching the VIA descriptors on the NIC can reduce the communication latency. We also present an efficient scheme for implementing the VIA doorbells without any hardware support. We provide a comprehensive performance evaluation study and discuss the impact of several hardware improvements on the performance of our implementation. The performance of the implemented VIA surpasses that of other existing software implementations of the VIA and is comparable to that of a hardware VIA implementation. The peak measured bandwidth for our system is 101.4MBytes/s and the one-way latency for short messages is 18.2μs.

Generalized least-squares polynomial preconditioners for symmetric indefinite linear equations

Polynomial preconditioners are frequently used in a parallel environment for the computation of the solution of large-scale sparse linear equations ( Ax= b) because of their easy implementation and trivial parallelization. With respect to symmetrical indefinite (SID) linear systems, the use of generalized least-squares (GLS) polynomial preconditioning is preferable to other polynomial preconditioning methods because of the ability to use a three-term recurrence relationship and the low implementation costs. The GLS preconditioning polynomial and its influence on the flexible generalized minimized residual (FGMRES) solver are discussed in this paper. The orthogonal polynomials required in the solution of the least-squares approximation problem are constructed using the Stieltjes procedure in multiple disjoint intervals which exclude the origin. The time-consuming numerical integration associated with this procedure is computed efficiently using Chebyshev polynomials of the first kind and the GLS polynomial reconditioned FGMRES algorithm is implemented using MPI in a highly parallel IBM SP2 environment. Experimental results using classical benchmark systems are presented and compared with those obtained using the recently developed SPAI preconditioned Bi-CGSTAB iterative method. The performance of the GLS preconditioned FGMRES solver is critically accessed.

Linear-scaling density-functional-theory calculations of electronic structure based on real-space grids: design, analysis, and scalability test of parallel algorithms

We have implemented parallel algorithms for density-functional-theory (DFT) based electronic-structure calculations. These include a plane-wave based algorithm, a real-space-grid algorithm based on a high-order finite difference method, and a linear-scaling real-space algorithm using localized orbitals. Parallelization schemes are described for these algorithms, and the computational complexity and the communications involved in the resulting parallel algorithms are analyzed. Scalability tests of these algorithms on massively parallel computers show that the linear-scaling DFT algorithm is highly scalable. For a 110,592-atom gallium arsenide system on 1024 IBM SP3 processors, the parallel efficiency is as high as 93%.

Compariing parallel multidimensional integration algorithms

We consider the design and implementation of new parallel globally adaptive multidimensional numerical integration algorithms. We describe a pair of algorithms (implemented in MPI on an IBM SP2); one uses a master-slave paradigm and the other is genuinely distributed. Then, we describe modifications of these algorithms to take advantage of a singularity detection scheme. The most obvious competitors for globally adaptive schemes in high dimensions are algorithms based on lattice rules. We compare our globally adaptive strategy with the straightforward implementation of a lattice rule algorithm paying attention to the quality of the respective error estimates. The lattice rule algorithm implementations involve computing new seeds for larger numbers of points than used previously, a significant though straightforward computational task.