We develop efficient algorithms for reliable and accurate evaluatins of the complex arcsine and a... more We develop efficient algorithms for reliable and accurate evaluatins of the complex arcsine and arccosine functions. A tight error bound is derived for each algorithm; the results are valid for all machine-representable points in the complex plane. The algorithms are presented in a pseudocode that has a convenient exception-handling facility. Corresponding Fortran 77 programs for an IEEE environment have also been developed to illustrate the practicality of the algorithms, and these programs have been tested very carefully to help confirm the correctness of the algorithms and their error bounds. The results of these tests are included in the article, but the Fortran 77 programs are not (these programs are available from Fairgrieve). Tests of other widely available programs fail at many points in the complex plane, and otherwise are slower and produce much less accurate results.
Proceedings of the 2003 workshop on Computer architecture education Held in conjunction with the 30th International Symposium on Computer Architecture - WCAE '03, 2003
The Intel® Itanium® architecture is increasingly becoming one of the major processor architecture... more The Intel® Itanium® architecture is increasingly becoming one of the major processor architectures present in the market today. Launched in 2001, the Intel Itanium processor was followed in 2002 by the Itanium 2 processor, with increased integer and floating-point performance. Measured by the SPEC CINT2000 benchmarks, the Itanium 2 processor still trails by about 25% the Intel P4 processor in integer performance, albeit P4 runs at more than three times Itanium's clock frequency. However, its floating-point performance clearly leads in the SPEC CFP2000 charts, and its rating is about 25% higher than that of the P4 processor. While the general features of the Itanium architecture such as large register sets, predication, speculation, and support for explicit parallelism [1] have been presented in several papers, books, and mainstream college textbooks [2], its floating-point architecture has been less publicized. Two books, [3] and [4], cover well this area. The present paper focuses on the floating-point architecture of the Itanium processor family, and points out a few remarkable features suitable to be the focus of a lecture, lab session, or project in a computer architecture class.
[1991] Proceedings 10th IEEE Symposium on Computer Arithmetic
Table-lookup algorithms for calculating elementary functions offer superior speed and accuracy wh... more Table-lookup algorithms for calculating elementary functions offer superior speed and accuracy when compared with more traditional algorithms. \Vith careful design, we show that it is feasible to implement table-lookup algorithms in hardware. Furthermore, we present a uniform approach to carry out tight error analysis for such implementations. I _ I or reproduce the publir, hedothers (ormto do Ofso, thiS[for[ /contribution, or allow | ' _'" / II c; Gnvernrnenl purposes. ..__]
This paper presents a generalization of incremental condition estimation, a technique for trackin... more This paper presents a generalization of incremental condition estimation, a technique for tracking the extremal singular values of a triangular matrix. While the original approach allowed for the estimation of the largest or smallest singular value, the generalized scheme allows for the estimation of any number of extremal singular values. For example, we can derive estimates for the three smallest singular values and the corresponding singular vectors at the same time. When estimating k singular values at the same time, the cost of one step of our generalized scheme on an n n matrix is Onk 2. Experimental results show that the resulting estimator does a good job of estimating the extremal singular values of triangular matrices and that, in particular, it leads to an inexpensive, y et very accurate and robust condition estimator.
We consider the solutions of Hermitian Toeplitz systems where the Toeplitz matrices are generated... more We consider the solutions of Hermitian Toeplitz systems where the Toeplitz matrices are generated by nonnegative functions f. The preconditioned conjugate gradient method with wellknown circulant preconditioners fails in the case when f has zeros. In this paper, we employ Toeplitz matrices of xed band-width as preconditioners. Their generating functions g are trigonometric polynomials of xed degree and are determined by minimizing the maximum relative error jj(f ?g)=fjj1. We show that the condition number of systems preconditioned by the band-Toeplitz matrices are O(1) for f with or without zeros. When f is positive, our preconditioned systems converge at the same rate as other well-known circulant preconditioned systems. We also give an a priori bound of the number of iterations required for convergence.
The 64-bit Intel® Itanium® architecture is designed for high-performance scientific and enterpris... more The 64-bit Intel® Itanium® architecture is designed for high-performance scientific and enterprise computing, and the Itanium processor is its first silicon implementation. Features such as extensive arithmetic support, predication, speculation, and explicit parallelism can be used to provide a sound infrastructure for supercomputing. A large number of high-performance computer companies are offering Itanium® -based systems, some capable of peak performance exceeding 50 GFLOPS. In this paper we give an overview of the most relevant architectural features and provide illustrations of how these features are used in both low-level and high-level support for scientific and engineering computing, including transcendental functions and linear algebra kernels.
18th IEEE Symposium on Computer Arithmetic (ARITH '07), 2007
Most implementations of the modular exponentiation, M E mod N , computation in cryptographic algo... more Most implementations of the modular exponentiation, M E mod N , computation in cryptographic algorithms employ Montgomery multiplication, ABR −1 mod N , instead of modular multiplication, AB mod N , even the former requires some transformational overheads. This is so because a state-of-the-art Montgomery multiplication implementation has a performance advantage over direct modular multiplication based on the Barrett algorithm that more than compensates for the overhead. In this paper, we present a direct modular multiplication method that is comparable in speed to Montgomery multiplication. One consequence is that when the exponent in small, direct computation (which does not incur the transformational overhead) using the modular multiplication algorithm presented here results in practical performance gain. For the exponent 17, for instance, which requires five modular multiplication, a saving of up to 40% can be achieved.
The Fast Fourier Transform (FFT) algorithm that calculates the Oiscrete Fourier Transform (OFT) i... more The Fast Fourier Transform (FFT) algorithm that calculates the Oiscrete Fourier Transform (OFT) is one of the major breakthrough in scientific computing and is now an indispensable tool in a vast number of fields. Unfortunately, software that provide fast computation of OFT via FFT differ vastly in functionality as well as uniformity. A widely accepted Applications Programmer Interface (API) for OFT would advance the field of scientific computing significantly. In this paper, we formulate an API for OFT computation that encompasses all the functionality that are offered by a number of popular packages combined, allows easy porting from existing codes, and exhibits a systematic naming convention with relatively short calling sequences.
SIAM Journal on Matrix Analysis and Applications, 1990
This paper presents an improved version of incremental condition estimation, a technique for trac... more This paper presents an improved version of incremental condition estimation, a technique for tracking the extremal singular values of a triangular matrix as it is being constructed one column at a time. We present a new motivation for this estimation technique using orthogonal projections. The paper focuses on an implementation of this estimation scheme in an accurate and consistent fashion. In particular, we address the subtle numerical issues arising in the computation of the eigensystem of a symmetric rank-one perturbed diagonal 22 matrix. Experimental results show that the resulting scheme does a good job in estimating the extremal singular values of triangular matrices, independent of matrix size and matrix condition number, and that it performs qualitatively in the same fashion as some of the commonly used nonincremental condition estimation schemes.
We propose a new algorithm for finding best minimax polynomial approximations in the complex plan... more We propose a new algorithm for finding best minimax polynomial approximations in the complex plane. The algorithm is the first satisfactory generalization of the well-known Remez algorithm for real approximations. Among all available algorithms, ours is the only quadratically convergent one. Numerical examples are presented to illustrate rapid convergence.
The IEEE Standard 754-1985 for Binary Floating-Point Arithmetic [1] was revised [2], and an impor... more The IEEE Standard 754-1985 for Binary Floating-Point Arithmetic [1] was revised [2], and an important addition is the definition of decimal floating-point arithmetic. This is intended mainly to provide a robust, reliable framework for financial applications that are often subject to legal requirements concerning rounding and precision of the results, because the binary floating-point arithmetic may introduce small but unacceptable errors. Using binary floating-point calculations to emulate decimal calculations in order to correct this issue has led to the existence of numerous proprietary software packages, each with its own characteristics and capabilities. IEEE 754R decimal arithmetic should unify the ways decimal floating-point calculations are carried out on various platforms. New algorithms and properties are presented in this paper which are used in a software implementation of the IEEE 754R decimal floatingpoint arithmetic, with emphasis on using binary operations efficiently. The focus is on rounding techniques for decimal values stored in binary format, but algorithms for the more important or interesting operations of addition, multiplication, division, and conversions between binary and decimal floating-point formats are also outlined. Performance results are included for a wider range of operations, showing promise that our approach is viable for applications that require decimal floating-point calculations.
high performance computing symposium, Apr 13, 2014
Calculating portions of eigenvalues and eigenvectors of matrices or matrix pencils has many appli... more Calculating portions of eigenvalues and eigenvectors of matrices or matrix pencils has many applications. An approach to this calculation for Hermitian problems based on a density matrix has been proposed in 2009 and a software package called FEAST has been developed. The density-matrix approach allows FEAST's implementation to exploit a key strength of modern computer architectures, namely, multiple levels of parallelism. Consequently, the software package has been well received and subsequently commercialized. A detailed theoretical analysis of Hermitian FEAST has also been established very recently. This paper generalizes the FEAST algorithm and theory, for the first time, to tackle non-Hermitian problems. Fundamentally, the new algorithm is basic subspace iteration or Bauer bi-iteration, except applied with a novel accelerator based on Cauchy integrals. The resulting algorithm retains the multi-level parallelism of Hermitian FEAST, making it a valuable new tool for large-scale computational science and engineering problems on leading-edge computing platforms.
A detailed new upgrade of the FEAST eigensolver targeting non-Hermitian eigenvalue problems is pr... more A detailed new upgrade of the FEAST eigensolver targeting non-Hermitian eigenvalue problems is presented and thoroughly discussed. It aims at broadening the class of eigenproblems that can be addressed within the framework of the FEAST algorithm. The algorithm is ideally suited for computing selected interior eigenvalues and their associated right/left bi-orthogonal eigenvectors, located within a subset of the complex plane. It combines subspace iteration with efficient contour integration techniques that approximate the left and right spectral projectors. We discuss the various algorithmic choices that have been made to improve the stability and usability of the new non-Hermitian eigensolver. The latter retains the convergence property and multi-level parallelism of Hermitian FEAST, making it a valuable new software tool for the scientific community.
New microprocessor architectures often require software support for basic arithmetic operations s... more New microprocessor architectures often require software support for basic arithmetic operations such as divide, or square root. The Intel R XScale processor, designed for low power mobile devices, provides no hardware support for floating-point. We show that an efficient software imple- mentation of the basic operations and math library routines can achieve competitive performance, and effectively hide the lack of
Table-driven techniques can be used to test highly accurate implementation of EXP LOG. The larges... more Table-driven techniques can be used to test highly accurate implementation of EXP LOG. The largest error observeed in EXP and LOG accurately to within 1/500 unit in the last place are reported in our tests. Methods to verify the tests' reliability are discussed. Results of applying the tests to our own as well as to a number of other implementations of EXP and LOG are presented.
The IA-64 architecture provides new opportunities and challenges for implementing an improved set... more The IA-64 architecture provides new opportunities and challenges for implementing an improved set of transcendental functions. Using several novel polynomial-based table-driven techniques, we are able to provide new algorithms for the transcendental functions. Major improvements include an accuracy level of about 0.6 ulps (units in the last place) and forward trigonometric functions that have a period of 2 π. The accuracy enhancements are achieved at improved speed, yet without an increase in the table size. In this paper, we highlight the key IA-64 architectural features that influenced our designs and explain the main ideas used in our new algorithms.
We develop efficient algorithms for reliable and accurate evaluatins of the complex arcsine and a... more We develop efficient algorithms for reliable and accurate evaluatins of the complex arcsine and arccosine functions. A tight error bound is derived for each algorithm; the results are valid for all machine-representable points in the complex plane. The algorithms are presented in a pseudocode that has a convenient exception-handling facility. Corresponding Fortran 77 programs for an IEEE environment have also been developed to illustrate the practicality of the algorithms, and these programs have been tested very carefully to help confirm the correctness of the algorithms and their error bounds. The results of these tests are included in the article, but the Fortran 77 programs are not (these programs are available from Fairgrieve). Tests of other widely available programs fail at many points in the complex plane, and otherwise are slower and produce much less accurate results.
Proceedings of the 2003 workshop on Computer architecture education Held in conjunction with the 30th International Symposium on Computer Architecture - WCAE '03, 2003
The Intel® Itanium® architecture is increasingly becoming one of the major processor architecture... more The Intel® Itanium® architecture is increasingly becoming one of the major processor architectures present in the market today. Launched in 2001, the Intel Itanium processor was followed in 2002 by the Itanium 2 processor, with increased integer and floating-point performance. Measured by the SPEC CINT2000 benchmarks, the Itanium 2 processor still trails by about 25% the Intel P4 processor in integer performance, albeit P4 runs at more than three times Itanium's clock frequency. However, its floating-point performance clearly leads in the SPEC CFP2000 charts, and its rating is about 25% higher than that of the P4 processor. While the general features of the Itanium architecture such as large register sets, predication, speculation, and support for explicit parallelism [1] have been presented in several papers, books, and mainstream college textbooks [2], its floating-point architecture has been less publicized. Two books, [3] and [4], cover well this area. The present paper focuses on the floating-point architecture of the Itanium processor family, and points out a few remarkable features suitable to be the focus of a lecture, lab session, or project in a computer architecture class.
[1991] Proceedings 10th IEEE Symposium on Computer Arithmetic
Table-lookup algorithms for calculating elementary functions offer superior speed and accuracy wh... more Table-lookup algorithms for calculating elementary functions offer superior speed and accuracy when compared with more traditional algorithms. \Vith careful design, we show that it is feasible to implement table-lookup algorithms in hardware. Furthermore, we present a uniform approach to carry out tight error analysis for such implementations. I _ I or reproduce the publir, hedothers (ormto do Ofso, thiS[for[ /contribution, or allow | ' _'" / II c; Gnvernrnenl purposes. ..__]
This paper presents a generalization of incremental condition estimation, a technique for trackin... more This paper presents a generalization of incremental condition estimation, a technique for tracking the extremal singular values of a triangular matrix. While the original approach allowed for the estimation of the largest or smallest singular value, the generalized scheme allows for the estimation of any number of extremal singular values. For example, we can derive estimates for the three smallest singular values and the corresponding singular vectors at the same time. When estimating k singular values at the same time, the cost of one step of our generalized scheme on an n n matrix is Onk 2. Experimental results show that the resulting estimator does a good job of estimating the extremal singular values of triangular matrices and that, in particular, it leads to an inexpensive, y et very accurate and robust condition estimator.
We consider the solutions of Hermitian Toeplitz systems where the Toeplitz matrices are generated... more We consider the solutions of Hermitian Toeplitz systems where the Toeplitz matrices are generated by nonnegative functions f. The preconditioned conjugate gradient method with wellknown circulant preconditioners fails in the case when f has zeros. In this paper, we employ Toeplitz matrices of xed band-width as preconditioners. Their generating functions g are trigonometric polynomials of xed degree and are determined by minimizing the maximum relative error jj(f ?g)=fjj1. We show that the condition number of systems preconditioned by the band-Toeplitz matrices are O(1) for f with or without zeros. When f is positive, our preconditioned systems converge at the same rate as other well-known circulant preconditioned systems. We also give an a priori bound of the number of iterations required for convergence.
The 64-bit Intel® Itanium® architecture is designed for high-performance scientific and enterpris... more The 64-bit Intel® Itanium® architecture is designed for high-performance scientific and enterprise computing, and the Itanium processor is its first silicon implementation. Features such as extensive arithmetic support, predication, speculation, and explicit parallelism can be used to provide a sound infrastructure for supercomputing. A large number of high-performance computer companies are offering Itanium® -based systems, some capable of peak performance exceeding 50 GFLOPS. In this paper we give an overview of the most relevant architectural features and provide illustrations of how these features are used in both low-level and high-level support for scientific and engineering computing, including transcendental functions and linear algebra kernels.
18th IEEE Symposium on Computer Arithmetic (ARITH '07), 2007
Most implementations of the modular exponentiation, M E mod N , computation in cryptographic algo... more Most implementations of the modular exponentiation, M E mod N , computation in cryptographic algorithms employ Montgomery multiplication, ABR −1 mod N , instead of modular multiplication, AB mod N , even the former requires some transformational overheads. This is so because a state-of-the-art Montgomery multiplication implementation has a performance advantage over direct modular multiplication based on the Barrett algorithm that more than compensates for the overhead. In this paper, we present a direct modular multiplication method that is comparable in speed to Montgomery multiplication. One consequence is that when the exponent in small, direct computation (which does not incur the transformational overhead) using the modular multiplication algorithm presented here results in practical performance gain. For the exponent 17, for instance, which requires five modular multiplication, a saving of up to 40% can be achieved.
The Fast Fourier Transform (FFT) algorithm that calculates the Oiscrete Fourier Transform (OFT) i... more The Fast Fourier Transform (FFT) algorithm that calculates the Oiscrete Fourier Transform (OFT) is one of the major breakthrough in scientific computing and is now an indispensable tool in a vast number of fields. Unfortunately, software that provide fast computation of OFT via FFT differ vastly in functionality as well as uniformity. A widely accepted Applications Programmer Interface (API) for OFT would advance the field of scientific computing significantly. In this paper, we formulate an API for OFT computation that encompasses all the functionality that are offered by a number of popular packages combined, allows easy porting from existing codes, and exhibits a systematic naming convention with relatively short calling sequences.
SIAM Journal on Matrix Analysis and Applications, 1990
This paper presents an improved version of incremental condition estimation, a technique for trac... more This paper presents an improved version of incremental condition estimation, a technique for tracking the extremal singular values of a triangular matrix as it is being constructed one column at a time. We present a new motivation for this estimation technique using orthogonal projections. The paper focuses on an implementation of this estimation scheme in an accurate and consistent fashion. In particular, we address the subtle numerical issues arising in the computation of the eigensystem of a symmetric rank-one perturbed diagonal 22 matrix. Experimental results show that the resulting scheme does a good job in estimating the extremal singular values of triangular matrices, independent of matrix size and matrix condition number, and that it performs qualitatively in the same fashion as some of the commonly used nonincremental condition estimation schemes.
We propose a new algorithm for finding best minimax polynomial approximations in the complex plan... more We propose a new algorithm for finding best minimax polynomial approximations in the complex plane. The algorithm is the first satisfactory generalization of the well-known Remez algorithm for real approximations. Among all available algorithms, ours is the only quadratically convergent one. Numerical examples are presented to illustrate rapid convergence.
The IEEE Standard 754-1985 for Binary Floating-Point Arithmetic [1] was revised [2], and an impor... more The IEEE Standard 754-1985 for Binary Floating-Point Arithmetic [1] was revised [2], and an important addition is the definition of decimal floating-point arithmetic. This is intended mainly to provide a robust, reliable framework for financial applications that are often subject to legal requirements concerning rounding and precision of the results, because the binary floating-point arithmetic may introduce small but unacceptable errors. Using binary floating-point calculations to emulate decimal calculations in order to correct this issue has led to the existence of numerous proprietary software packages, each with its own characteristics and capabilities. IEEE 754R decimal arithmetic should unify the ways decimal floating-point calculations are carried out on various platforms. New algorithms and properties are presented in this paper which are used in a software implementation of the IEEE 754R decimal floatingpoint arithmetic, with emphasis on using binary operations efficiently. The focus is on rounding techniques for decimal values stored in binary format, but algorithms for the more important or interesting operations of addition, multiplication, division, and conversions between binary and decimal floating-point formats are also outlined. Performance results are included for a wider range of operations, showing promise that our approach is viable for applications that require decimal floating-point calculations.
high performance computing symposium, Apr 13, 2014
Calculating portions of eigenvalues and eigenvectors of matrices or matrix pencils has many appli... more Calculating portions of eigenvalues and eigenvectors of matrices or matrix pencils has many applications. An approach to this calculation for Hermitian problems based on a density matrix has been proposed in 2009 and a software package called FEAST has been developed. The density-matrix approach allows FEAST's implementation to exploit a key strength of modern computer architectures, namely, multiple levels of parallelism. Consequently, the software package has been well received and subsequently commercialized. A detailed theoretical analysis of Hermitian FEAST has also been established very recently. This paper generalizes the FEAST algorithm and theory, for the first time, to tackle non-Hermitian problems. Fundamentally, the new algorithm is basic subspace iteration or Bauer bi-iteration, except applied with a novel accelerator based on Cauchy integrals. The resulting algorithm retains the multi-level parallelism of Hermitian FEAST, making it a valuable new tool for large-scale computational science and engineering problems on leading-edge computing platforms.
A detailed new upgrade of the FEAST eigensolver targeting non-Hermitian eigenvalue problems is pr... more A detailed new upgrade of the FEAST eigensolver targeting non-Hermitian eigenvalue problems is presented and thoroughly discussed. It aims at broadening the class of eigenproblems that can be addressed within the framework of the FEAST algorithm. The algorithm is ideally suited for computing selected interior eigenvalues and their associated right/left bi-orthogonal eigenvectors, located within a subset of the complex plane. It combines subspace iteration with efficient contour integration techniques that approximate the left and right spectral projectors. We discuss the various algorithmic choices that have been made to improve the stability and usability of the new non-Hermitian eigensolver. The latter retains the convergence property and multi-level parallelism of Hermitian FEAST, making it a valuable new software tool for the scientific community.
New microprocessor architectures often require software support for basic arithmetic operations s... more New microprocessor architectures often require software support for basic arithmetic operations such as divide, or square root. The Intel R XScale processor, designed for low power mobile devices, provides no hardware support for floating-point. We show that an efficient software imple- mentation of the basic operations and math library routines can achieve competitive performance, and effectively hide the lack of
Table-driven techniques can be used to test highly accurate implementation of EXP LOG. The larges... more Table-driven techniques can be used to test highly accurate implementation of EXP LOG. The largest error observeed in EXP and LOG accurately to within 1/500 unit in the last place are reported in our tests. Methods to verify the tests' reliability are discussed. Results of applying the tests to our own as well as to a number of other implementations of EXP and LOG are presented.
The IA-64 architecture provides new opportunities and challenges for implementing an improved set... more The IA-64 architecture provides new opportunities and challenges for implementing an improved set of transcendental functions. Using several novel polynomial-based table-driven techniques, we are able to provide new algorithms for the transcendental functions. Major improvements include an accuracy level of about 0.6 ulps (units in the last place) and forward trigonometric functions that have a period of 2 π. The accuracy enhancements are achieved at improved speed, yet without an increase in the table size. In this paper, we highlight the key IA-64 architectural features that influenced our designs and explain the main ideas used in our new algorithms.
Uploads
Papers by Ping Tang