LAPACK Routines Research Articles

A Set of Batched Basic Linear Algebra Subprograms and LAPACK Routines

This article describes a standard API for a set of Batched Basic Linear Algebra Subprograms (Batched BLAS or BBLAS). The focus is on many independent BLAS operations on small matrices that are grouped together and processed by a single routine, called a Batched BLAS routine. The matrices are grouped together in uniformly sized groups, with just one group if all the matrices are of equal size. The aim is to provide more efficient, but portable, implementations of algorithms on high-performance many-core platforms. These include multicore and many-core CPU processors, GPUs and coprocessors, and other hardware accelerators with floating-point compute facility. As well as the standard types of single and double precision, we also include half and quadruple precision in the standard. In particular, half precision is used in many very large scale applications, such as those associated with machine learning.

SLASs: A fully automatic auto-tuned linear algebra library based on OpenMP extensions implemented in OmpSs (LASs Library)

In this work we have implemented a novel Linear Algebra Library on top of the task-based runtime OmpSs-2. We have used some of the most advanced OmpSs-2 features; weak dependencies and regions, together with the final clause for the implementation of auto-tunable code for the BLAS-3 trsm routine and the LAPACK routines npgetrf and npgesv. All these implementations are part of the first prototype of sLASs library, a novel library for auto-tunable codes for linear algebra operations based on LASs library. In all these cases, the use of the OmpSs-2 features presents an improvement in terms of execution time against other reference libraries such as, the original LASs library, PLASMA, ATLAS and Intel MKL. These codes are able to reduce the execution time in about 18% on big matrices, by increasing the IPC on gemm and reducing the time of task instantiation. For a few medium matrices, benefits are also seen. For small matrices and a subset of medium matrices, specific optimizations that allow to increase the degree of parallelism in both, gemm and trsm tasks, are applied. This strategy achieves an increment in performance of up to 40%.

BLASFEO

Basic Linear Algebra Subroutines for Embedded Optimization (BLASFEO) is a dense linear algebra library providing high-performance implementations of BLAS- and LAPACK-like routines for use in embedded optimization and small-scale high-performance computing, in general. A key difference with respect to existing high-performance implementations of BLAS is that the computational performance is optimized for small- to medium-scale matrices, i.e., for sizes up to a few hundred. BLASFEO comes with three different implementations: a high-performance implementation aimed at providing the highest performance for matrices fitting in cache, a reference implementation providing portability and embeddability and optimized for very small matrices, and a wrapper to standard BLAS and LAPACK providing high performance on large matrices. The three implementations of BLASFEO together provide high-performance dense linear algebra routines for matrices ranging from very small to large. Compared to both open-source and proprietary highly tuned BLAS libraries, for matrices of size up to about 100, the high-performance implementation of BLASFEO is about 20--30% faster than the corresponding level 3 BLAS routines and two to three times faster than the corresponding LAPACK routines.

Multiple Response Variables Regression Models in R: The mcglm Package

This article describes the R package mcglm implemented for fitting multivariate covariance generalized linear models (McGLMs). McGLMs provide a general statistical modeling framework for normal and non-normal multivariate data analysis, designed to handle multivariate response variables, along with a wide range of temporal and spatial correlation structures defined in terms of a covariance link function and a matrix linear predictor involving known symmetric matrices. The models take non-normality into account in the conventional way by means of a variance function, and the mean structure is modeled by means of a link function and a linear predictor. The models are fitted using an estimating function approach based on second-moment assumptions. This provides a unified approach to a wide variety of different types of response variables and covariance structures, including multivariate extensions of repeated measures, time series, longitudinal, genetic, spatial and spatio-temporal structures. The mcglm package allows a flexible specification of the mean and covariance structures, and explicitly deals with multivariate response variables, through a user friendly formula interface similar to the ordinary glm function. Illustrations in this article cover a wide range of applications from the traditional one response variable Gaussian mixed models to multivariate spatial models for areal data using the multivariate Tweedie distribution. Additional features, such as robust and bias-corrected standard errors for regression parameters, residual analysis, measures of goodness-of-fit and model selection using the score information criterion are discussed through six worked examples. The mcglm package is a full R implementation based on the Matrix package which provides efficient access to BLAS (basic linear algebra subroutines), Lapack (dense matrix), TAUCS (sparse matrix) and UMFPACK (sparse matrix) routines for efficient linear algebra in R.

OptQC v1.3: An (updated) optimized parallel quantum compiler

We present a revised version of the OptQC program of Loke et al. (2014) [1]. We have removed the simulated annealing process in favour of a descending random walk. We have also introduced a new method for iteratively generating permutation matrices during the random walk process, providing a reduced total cost for implementing the quantum circuit. Lastly, we have also added a synchronization mechanism between threads, giving quicker convergence to more optimal solutions. New version program summaryProgram title: OptQC v1.3Catalogue identifier: AEUA_v1_3Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEUA_v1_3.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 240903No. of bytes in distributed program, including test data, etc.: 632395Distribution format: tar.gzProgramming language: Fortran, MPI.Computer: Any computer with Fortran compiler (not gfortran4.9 or earlier) and MPI library.Operating system: Linux.Classification: 4.15.Catalogue identifier of previous version: AEUA_v1_3Journal reference of previous version: Comput. Phys. Comm. 185(2014)3307External routines: Intel MKL LAPACK routines and MPI routines.Does the new version supersede the previous version?: YesNature of problem:It aims to minimize the number of quantum gates required to implement a given unitary operation.Solution method:It utilizes a descending random walk to select permutation matrices P and Q for a given unitary matrix U such that the number of gates in the quantum circuit of U=QTPTU′PQ is minimized, where U′ is equivalent to U up to a permutation. The decomposition of a unitary operator is performed by recursively applying the cosine–sine decomposition.Reasons for new version:Simulated annealing process was found to give suboptimal results compared to a normal descending random walk. Computation time was also bloated by the necessity of running the CS decomposition thrice (for U′,P and PT) for each iteration of the optimization process.Summary of revisions:•Simulated annealing process was replaced by a descending random walk (equivalent to setting the threshold value β to 0), due to poor convergence of the simulated annealing process, and due to the lack of pathological instances of local minima in this particular search space.•Introduced an iterative method for generating the general permutation matrix P, in which we start with P=I and build up the quantum circuit (and modify P correspondingly) by adding gates onto the existing circuit. This removes the requirement of running the CS decomposition on P and PT to find the quantum circuit implementation, since it is already known by construction.•Added a synchronization mechanism between threads (after some prescribed number of iterations) in which the current state of the top 10% processes with the fittest solutions is copied over to the remaining 90%. This works so as to discard the less fit solutions and focuses the searching algorithm in the state space with the fittest solutions.Additional comments:The program contains some Fortran2003 features and will not compile with gcc4.9 or earlier.Running time:As before, running time increases with the size of the unitary matrix, as well as the prescribed maximum number of iterations for qubit permutation selection and the descending random walk. All simulation results presented in this paper are obtained from running the program on the Fornax supercomputer managed by iVEC@UWA with Intel Xeon X5650 CPUs. A comparison of running times is also given in Table 1.

MolSOC: A spin–orbit coupling code

A description of the theoretical background, the numerical algorithm and the input/output files structures of the computer program MolSOC is presented. This code provides an efficient tool for the spin–orbit (SO) coupling treatment in many-electron systems by a non-iterative method implementation. Preliminary versions of the MolSOC program have been already tested and successfully applied in the framework of density functional theory (DFT). SO calculations can be performed either employing a full SO operator or a reduced screened nuclear-charge operator. With the second option, the required computational time can be significantly reduced. However, it has been previously demonstrated that the full SO operator can be successfully applied to extended molecular systems. Finally, practical instructions for the program use are also included. Program summaryProgram title: MolSOC v0.1Catalogue identifier: AERK_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AERK_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 13333No. of bytes in distributed program, including test data, etc.: 2301930Distribution format: tar.gzProgramming language: Fortran 90.Computer: All Linux based workstations.Operating system: Linux.Classification: 16.1.External routines: Subset of BLAS and LAPACK routines (included in the distribution file).Nature of problem:An accurate description of the electronic structure of heavy many electron systems requires, very often, the inclusion of the spin–orbit coupling contributions.Solution method:Spin–Orbit matrix elements can be computed between optimized states at the density functional theory level using the full Breit–Pauli operator or a reduced screened nuclear-charge operator.Running time:It depends on the size of the chemical system, on the basis set used and on the computational method used. The tests included only take seconds or a few minutes to run.

CMIstark: Python package for the Stark-effect calculation and symmetry classification of linear, symmetric and asymmetric top wavefunctions in dc electric fields

The Controlled Molecule Imaging group (CMI) at the Center for Free Electron Laser Science (CFEL) has developed the CMIstark software to calculate, view, and analyze the energy levels of adiabatic Stark energy curves of linear, symmetric top and asymmetric top molecules. The program exploits the symmetry of the Hamiltonian to generate fully labeled adiabatic Stark energy curves.CMIstark is written in Python and easily extendable, while the core numerical calculations make use of machine optimized BLAS and LAPACK routines. Calculated energies are stored in HDF5 files for convenient access and programs to extract ASCII data or to generate graphical plots are provided. Program summaryProgram title: CMIstarkCatalogue identifier: AEQS_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEQS_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: GNU General Public License version 3 or later with amendments. See code for details.No. of lines in distributed program, including test data, etc.: 4091No. of bytes in distributed program, including test data, etc.: 62805Distribution format: tar.gzProgramming language: Python (version 2.6.x, 2.7.x).Computer: Any Macintosh, PC, or Linux/UNIX workstations with a modern Python distribution.Operating system: Tested on Mac OS X and a variety of Linux distributions.RAM: 2 GB for typical calculationsClassification: 16.1.External routines:Python packages numpy and scipy; utilizes (optimized) LAPACK and BLAS through scipy. All packages available under open-source licenses.Nature of problem:Calculation of the Stark effect of asymmetric top molecules in arbitrarily strong dc electric fields in a correct symmetry classification and using correct labeling of the adiabatic Stark curves.Solution method:We set up the full M matrices of the quantum-mechanical Hamiltonian in the basis set of symmetric top wavefunctions and, subsequently, Wang transform the Hamiltonian matrix. We separate, as far as possible, the sub-matrices according to the remaining symmetry, and then diagonalize the individual blocks. This application of the symmetry consideration to the Hamiltonian allows an adiabatic correlation of the asymmetric top eigenstates in the dc electric field to the field-free eigenstates. This directly yields correct adiabatic state labels and, correspondingly, adiabatic Stark energy curves.Restrictions:The maximum value of J is limited by the available main memory. A modern desktop computer with 16 GB of main memory allows for calculations including all Js up to a values larger than 100 even for the most complex cases of asymmetric tops.Running time:Typically 1 s–1 week on a single CPU or equivalent on multi-CPU systems (depending greatly on system size and RAM); parallelization through BLAS/LAPACK. For instance, calculating all energies up to J=25 of indole (vide infra) for one field strength takes 1 CPU-s on a current iMac.

A Cellwise Block-Gauss-Seidel Iterative Method for Multigroup SN Transport on a Hybrid Parallel Computer Architecture

A Fourier analysis is conducted in two-dimensional (2-D) geometry for the discrete ordinates (SN) approximation of the neutron transport problem solved with Richardson iteration (source iteration) using the cellwise block-Jacobi (bJ) and block-Gauss-Seidel (bGS) algorithms. The results of the Fourier analysis show that convergence of bJ can degrade, leading to a spectral radius ρ equal to 1, in problems containing optically thin cells. For problems containing cells that are optically thick, instead, ρ tends to 0. Hence, in the optically-thick-cell regime, bJ is rapidly convergent even for scattering-dominated problems. Similar conclusions hold for bGS, except bGS approaches the asymptotic, thick-cell regime at convergence rates higher than bJ. Hence, we have implemented the bGS algorithm on the Roadrunner hybrid, parallel computer architecture. A compute node of this massively parallel machine comprises AMD Opteron cores that are linked to a Cell Broadband Engine (Cell/B.E.). LAPACK routines have been ported to the Cell/B.E. in order to make use of its parallel synergistic processing elements (SPEs). The bGS algorithm is based on the LU factorization and solution of a linear system that couples the fluxes for all SN angles and energy groups on a mesh cell. For every cell of a mesh that has been parallel decomposed on the higher-level Opteron processors, a linear system is transferred to the Cell/B.E. and the parallel LAPACK routines are used to compute a solution, which is then transferred back to the Opteron, where the rest of the SN transport computations take place. Compared to standard parallel machines, a one-hundred-fold speedup of the bGS was observed on Roadrunner. Numerical experiments with strong and weak parallel scaling demonstrate that the bGS method is viable and compares favorably to full parallel transport sweeps (FPS) on 2-D unstructured meshes when it is applied to optically thick, multimaterial problems. Specifically, the strong parallel efficiency of accelerated bGS on Roadrunner can achieve 73% at 512 processors, compared with 32 processors, while efficiency is 34% for the (Opteron-only) implementation of FPS. The weak parallel efficiency of bGS is 58% while it reaches 10% for FPS. As expected, however, bGS is not as efficient as FPS in optically thin problems.

Scaling LAPACK panel operations using parallel cache assignment

In LAPACK many matrix operations are cast as block algorithms which iteratively process a panel using an unblocked algorithm and then update a remainder matrix using the high performance Level 3 BLAS. The Level 3 BLAS have excellent scaling, but panel processing tends to be bus bound, and thus scales with bus speed rather than the number of processors ( p ). Amdahl's law therefore ensures that as p grows, the panel computation will become the dominant cost of these LAPACK routines. Our contribution is a novel parallel cache assignment approach to the panel factorization which we show scales well with p . We apply this general approach to the QR, QL, RQ, LQ and LU panel factorizations. We show results for two commodity platforms: an 8-core Intel platform and a 32-core AMD platform. For both platforms and all twenty implementations (five factorizations each of which is available in 4 types), we present results that demonstrate that our approach yields significant speedup over the existing state of the art.

Level-3 Cholesky Factorization Routines Improve Performance of Many Cholesky Algorithms

Four routines called DPOTF3i, i = a,b,c,d, are presented. DPOTF3i are a novel type of level-3 BLAS for use by BPF ( B locked P acked F ormat) Cholesky factorization and LAPACK routine DPOTRF. Performance of routines DPOTF3i are still increasing when the performance of Level-2 routine DPOTF2 of LAPACK starts decreasing. This is our main result and it implies, due to the use of larger block size nb , that DGEMM, DSYRK, and DTRSM performance also increases! The four DPOTF3i routines use simple register blocking. Different platforms have different numbers of registers. Thus, our four routines have different register blocking sizes. BPF is introduced. LAPACK routines for POTRF and PPTRF using BPF instead of full and packed format are shown to be trivial modifications of LAPACK POTRF source codes. We call these codes BPTRF. There are two variants of BPF: lower and upper. Upper BPF is “identical” to Square Block Packed Format (SBPF). “LAPACK” implementations on multicore processors use SBPF. Lower BPF is less efficient than upper BPF. Vector inplace transposition converts lower BPF to upper BPF very efficiently. Corroborating performance results for DPOTF3i versus DPOTF2 on a variety of common platforms are given for n ≈ nb as well as results for large n comparing DBPTRF versus DPOTRF.

Empirical Modelling of Linear Algebra Shared-Memory Routines

In this work the behavior of the multithreaded implementation of some LAPACK routines on PLASMA and Intel MKL is analyzed. The main goal is to develop a methodology for the installation and modelling of shared-memory linear algebra routines so that some decisions to reduce the execution time can be taken at running time. Typical decisions are: the number of threads to use, the block or tile size in algorithms by blocks or tiles, and the routine to use when there are several algorithms or implementations to solve the problem available. Experiments carried out with PLASMA and Intel MKL show that decisions can be taken automatically and satisfactory execution times are obtained.

The Software Architecture for Performing Scientific Computation with the JLAPACK Libraries in ScalaLab

Although LAPACK is a powerful library its utilization is difficult. JLAPACK, a Java translation obtained automatically from the Fortran LAPACK sources, retains exactly the same difficult to use interface of LAPACK routines. The MTJ library implements an object oriented Java interface to JLAPACK that hides many complicated details. ScalaLab exploits the flexibility of the Scala language to present an even more friendly and convenient interface to the powerful but complicated JLAPACK library. The article describes the interfacing of the low-level JLAPACK routines within the ScalaLab environment. This is performed rather easily by exploiting well suited features of the Scala language. Also, the paper demonstrates the convenience of using JLAPACK routines for linear algebra operations from within ScalaLab.

Solvers for the verified solution of parametric linear systems

We present a newly developed version of our solvers for the verified solution of dense parametric linear systems, i.e. linear systems whose system matrix and right-hand side depend affine-linearly on parameters that vary inside prescribed intervals. The solvers use our C++ class library for reliable computing, C-XSC. The C-XSC library provides many features, especially easy to handle data types for dense and sparse matrices and vectors and the ability to compute dot products and dot product expressions in arbitrary precision. The new solvers can use either sparse or dense matrices as the coefficient matrices for the parameters. The use of sparse coefficient matrices can result in huge improvements in both performance and memory consumption. BLAS and LAPACK routines are used where applicable, and OpenMP is used for the parallelization on multi-core and multi-processor systems. The solvers also provide the ability to compute not only an outer but also a componentwise inner enclosure of the solution set of the system and to choose between two versions of the algorithm, one being very fast and one giving sharp results and extending the range of solvable systems. We give some examples for parametric linear systems (also from real world examples such as worst-case tolerance analysis of linear electric circuits), give performance measurements of our solvers and also demonstrate that they scale very well when using multiple cores or processors.

EPW: A program for calculating the electron–phonon coupling using maximally localized Wannier functions

EPW ( Electron– Phonon coupling using Wannier functions) is a program written in Fortran90 for calculating the electron–phonon coupling in periodic systems using density-functional perturbation theory and maximally localized Wannier functions. EPW can calculate electron–phonon interaction self-energies, electron–phonon spectral functions, and total as well as mode-resolved electron–phonon coupling strengths. The calculation of the electron–phonon coupling requires a very accurate sampling of electron–phonon scattering processes throughout the Brillouin zone, hence reliable calculations can be prohibitively time-consuming. EPW combines the Kohn–Sham electronic eigenstates and the vibrational eigenmodes provided by the Quantum ESPRESSO package (see Giannozzi et al., 2009 [1]) with the maximally localized Wannier functions provided by the wannier90 package (see Mostofi et al., 2008 [2]) in order to generate electron–phonon matrix elements on arbitrarily dense Brillouin zone grids using a generalized Fourier interpolation. This feature of EPW leads to fast and accurate calculations of the electron–phonon coupling, and enables the study of the electron–phonon coupling in large and complex systems. Program summary Program title: EPW Catalogue identifier: AEHA_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEHA_v1_0.html Program obtainable from: CPC Program Library, Queen's University, Belfast, N. Ireland Licensing provisions: GNU Public License No. of lines in distributed program, including test data, etc.: 304 443 No. of bytes in distributed program, including test data, etc.: 1 487 466 Distribution format: tar.gz Programming language: Fortran 90 Computer: Any architecture with a Fortran 90 compiler Operating system: Any environment with a Fortran 90 compiler Has the code been vectorized or parallelized?: Yes, optimized for 1 to 64 processors RAM: Heavily system dependent, as small as a few MB Supplementary material: A copy of the “EPW/examples” directory containing the phonon binary files can be downloaded Classification: 7 External routines: MPI, Quantum-ESPRESSO package [1], BLAS, LAPACK, FFTW. (The necessary Blas, Lapack and FFTW routines are included in the Quantum-ESPRESSO package [1].) Nature of problem: The calculation of the electron–phonon coupling from first-principles requires a very accurate sampling of electron–phonon scattering processes throughout the Brillouin zone; hence reliable calculations can be prohibitively timeconsuming. Solution method: EPW makes use of a real-space formulation and combines the Kohn–Sham electronic eigenstates and the vibrational eigenmodes provided by the Quantum-ESPRESSO package with the maximally localized Wannier functions provided by the wannier90 package in order to generate electron–phonon matrix elements on arbitrarily dense Brillouin zone grids using a generalized Fourier interpolation. Running time: Single processor examples typically take 5–10 minutes.

Rectangular full packed format for cholesky's algorithm

We describe a new data format for storing triangular, symmetric, and Hermitian matrices called Rectangular Full Packed Format (RFPF). The standard two-dimensional arrays of Fortran and C (also known as full format ) that are used to represent triangular and symmetric matrices waste nearly half of the storage space but provide high performance via the use of Level 3 BLAS. Standard packed format arrays fully utilize storage (array space) but provide low performance as there is no Level 3 packed BLAS. We combine the good features of packed and full storage using RFPF to obtain high performance via using Level 3 BLAS as RFPF is a standard full-format representation. Also, RFPF requires exactly the same minimal storage as packed the format. Each LAPACK full and/or packed triangular, symmetric, and Hermitian routine becomes a single new RFPF routine based on eight possible data layouts of RFPF. This new RFPF routine usually consists of two calls to the corresponding LAPACK full-format routine and two calls to Level 3 BLAS routines. This means no new software is required. As examples, we present LAPACK routines for Cholesky factorization, Cholesky solution, and Cholesky inverse computation in RFPF to illustrate this new work and to describe its performance on several commonly used computer platforms. Performance of LAPACK full routines using RFPF versus LAPACK full routines using the standard format for both serial and SMP parallel processing is about the same while using half the storage. Performance gains are roughly one to a factor of 43 for serial and one to a factor of 97 for SMP parallel times faster using vendor LAPACK full routines with RFPF than with using vendor and/or reference packed routines.

Scaling LAPACK panel operations using parallel cache assignment

In LAPACK many matrix operations are cast as block algorithms which iteratively process a panel using an unblocked algorithm and then update a remainder matrix using the high performance Level 3 BLAS. The Level~3 BLAS have excellent weak scaling, but panel processing tends to be bus bound, and thus scales with bus speed rather than the number of processors ( p ). Amdahl's law therefore ensures that as p grows, the panel computation will become the dominant cost of these LAPACK routines. Our contribution is a novel parallel cache assignment approach which we show scales well with p . We apply this general approach to the QR and LU panel factorizations on two commodity 8-core platforms with very different cache structures, and demonstrate superlinear panel factorization speedups on both machines. Other approaches to this problem demand complicated reformulations of the computational approach, new kernels to be tuned, new mathematics, an inflation of the high-order flop count, and do not perform as well. By demonstrating a straight-forward alternative that avoids all of these contortions and scales with p , we address a critical stumbling block for dense linear algebra in the age of massive parallelism.

Linear‐scaling Cholesky decomposition

We present linear-scaling routines for the calculation of the Cholesky decomposition of a symmetric positive-definite matrix and its inverse. As an example, we consider the inversion of the overlap matrix of DNA and amylose fragments as well as of linear alkanes, where the largest system corresponds to a 21,442 x 21,442 matrix. The efficiency and the scaling behavior are discussed and compared to standard LAPACK routines. Our Cholesky routines are publicly available on the web.

POTHMF: A program for computing potential curves and matrix elements of the coupled adiabatic radial equations for a hydrogen-like atom in a homogeneous magnetic field

A FORTRAN 77 program is presented which calculates with the relative machine precision potential curves and matrix elements of the coupled adiabatic radial equations for a hydrogen-like atom in a homogeneous magnetic field. The potential curves are eigenvalues corresponding to the angular oblate spheroidal functions that compose adiabatic basis which depends on the radial variable as a parameter. The matrix elements of radial coupling are integrals in angular variables of the following two types: product of angular functions and the first derivative of angular functions in parameter, and product of the first derivatives of angular functions in parameter, respectively. The program calculates also the angular part of the dipole transition matrix elements (in the length form) expressed as integrals in angular variables involving product of a dipole operator and angular functions. Moreover, the program calculates asymptotic regular and irregular matrix solutions of the coupled adiabatic radial equations at the end of interval in radial variable needed for solving a multi-channel scattering problem by the generalized R-matrix method. Potential curves and radial matrix elements computed by the POTHMF program can be used for solving the bound state and multi-channel scattering problems. As a test desk, the program is applied to the calculation of the energy values, a short-range reaction matrix and corresponding wave functions with the help of the KANTBP program. Benchmark calculations for the known photoionization cross-sections are presented.

Large-scale ab initio calculations based on three levels of parallelization

We suggest and implement a parallelization scheme based on an efficient multiband eigenvalue solver, called the locally optimal block preconditioned conjugate gradient ( lobpcg) method, and using an optimized three-dimensional (3D) fast Fourier transform (FFT) in the ab initio plane-wave code abinit. In addition to the standard data partitioning over processors corresponding to different k-points, we introduce data partitioning with respect to blocks of bands as well as spatial partitioning in the Fourier space of coefficients over the plane waves basis set used in abinit. This k-points-multiband-FFT parallelization avoids any collective communications on the whole set of processors relying instead on one-dimensional communications only. For a single k-point, super-linear scaling is achieved for up to 100 processors due to an extensive use of hardware-optimized blas, lapack and scalapack routines, mainly in the lobpcg routine. We observe good performance up to 200 processors. With 10 k-points our three-way data partitioning results in linear scaling up to 1000 processors for a practical system used for testing.

A precise boundary element method for macromolecular transport properties.

A very precise boundary element numerical solution of the exact formulation of the hydrodynamic resistance problem with stick boundary conditions is presented. BEST, the Fortran 77 program developed for this purpose, computes the full transport tensors in the center of resistance or the center of diffusion for an arbitrarily shaped rigid body, including rotation-translation coupling. The input for this program is a triangulation of the solvent-defined surface of the molecule of interest, given by Connolly's MSROLL or other suitable triangulator. The triangulation is prepared for BEST by COALESCE, a program that allows user control over the quality and number of triangles to describe the surface. High numerical precision is assured by effectively exact integration of the Oseen tensor over triangular surface elements, and by scaling the hydrodynamic computation to the precise surface area of the molecule. Efficiency of computation is achieved by the use of public domain LAPACK routines that call BLAS Level 3 hardware-optimized subroutines available for most processors. A protein computation can be done in less than 10 min of CPU time in a modern Pentium IV processor. The present work includes a complete analysis of the sources of error in the numerical work and techniques to eliminate these errors. The operation of BEST is illustrated with applications to ellipsoids of revolution, and Lysozyme, a small protein. The typical numerical accuracy achieved is 0.05% compared to analytical theory. The numerical precision for a protein is better than 1%, much better than experimental errors in these quantities, and more than 10 times better than traditional bead-based methods.