Search

A framework, PhaseLift, was recently proposed to solve the phase retrieval problem. In this framework, the problem is solved by optimizing a cost function over the set of complex Hermitian positive semidefinite matrices. This approach to phase retrieval motivates a more general consideration of optimizing cost functions on semidefinite Hermitian matrices where the desired minimizers are known to have low rank. This paper considers an approach based on an alternative cost function defined on a union of appropriate manifolds. It is related to the original cost function in a manner that preserves the ability to find a global minimizer and is significantly more efficient computationally. A rank-based optimality condition for stationary points is given and optimization algorithms based on state-of-the-art Riemannian optimization and dynamically reducing rank are proposed. Empirical evaluations are performed using the PhaseLift problem. The new approach is shown to be an effective method of phase retrieval with computational efficiency increased substantially compared to a state-of-the-art algorithm, the Wirtinger flow algorithm. A preliminary version of this paper can be found in [W. Huang, K. A. Gallivan, and X. Zhang, Procedia Comput. Sci., 80 (2016), pp. 1125-1134]. [ABSTRACT FROM AUTHOR]

Let A ∈ ℂN x N be a nonsingular complex matrix and b and c be complex vectors of length N. The goal of this paper is to investigate approaches for efficient approximations of the bilinear form c*A-1b. Equivalently, we wish to approximate the scalar value c*x, where x solves the linear system Ax = b. Here the matrix A can be very large or its elements can be too costly to compute so that A is not explicitly available and it is used only in the form of the matrix-vector product. Therefore a direct method is not an option. For A Hermitian positive definite, b*A-1b can be efficiently approximated as a by-product of the conjugate-gradient iterations, which is mathematically equivalent to the matching moment approximations computed via the Gauss-Christoffel quadrature. In this paper we propose a new method using the biconjugate gradient iterations which is applicable to the general complex case. The proposed approach will be compared with existing ones using analytic arguments and numerical experiments. [ABSTRACT FROM AUTHOR]

This paper concerns numerical methods for computing complex geometrical optics (CGO) solutions to the conductivity equation ∇ ⋅ σ∇u (⋅, k) = 0 in ℝ² for piecewise smooth conductivities σ, where k is a complex parameter. The key is to solve an R-linear singular integral equation defined in the unit disk. Recently, Astala et al. [Appl. Comput. Harmon. Anat., 29 (2010), pp. 2-17] proposed a complicated method for numerical computation of CGO solutions by solving a periodic version of the ℝ-linear integral equation in a rectangle containing the unit disk. In this paper, based on the fast algorithms in [P. Daripa and D. Mashat, Numer. Algorithms, 18 (1998), pp. 133-157] for singular integral transforms, we propose a simpler numerical method which solves the ℝ-linear integral equation in the unit disk directly. For the resulting ℝ-linear operator equation, a minimal residual iterative method is proposed. Numerical examples illustrate the accuracy and efficiency of the new method. [ABSTRACT FROM AUTHOR]

As a solution to the frequently reported problem of the high frequencies reconstructed sometimes many iterations after the low frequencies in tomography, we proposed in a recent paper [V. Israel-Jost, Ph. Choquet, and A. Constantinesco, Int. J. Biomed. Imaging, 2006, paper 34043.] a general adaptation that can apply to many algorithms: the incomplete backprojection. We stressed the fact that this type of frequency adapted (FA) algorithm works only when the linear problem of tomography to be solved involves deblurring, i.e., when the point spread functions (PSFs) associated with each voxel are large. Thus, depending on the desired level of detail, a parameter is set to achieve the reconstruction within a small number of iterations. The aim of this paper is to give a mathematical analysis of both the acceleration obtained and the convergence of an FA algorithm derived from the simultaneous algebraic reconstruction technique (SART) algorithm of Andersen and Kak, the so-called FA-SART algorithm. [ABSTRACT FROM AUTHOR]

This paper presents numerical results for the asynchronous version of the fast adaptive composite-grid algorithm (AFACx). These results confirm the level-independent convergence bounds established theoretically in a companion paper. These numerical results include the case of AFACx applied to first-order system least-squares finite element discretizations of the stationary Stokes equations on curvilinear adaptive mesh refinement grids. [ABSTRACT FROM AUTHOR]

In an earlier paper, we developed a new formulation of the shallow-water equations on the sphere, based on the semi-Lagrangian advection of potential vorticity. An important feature of our two-time-level model is that forward extrapolation of the nonlinear terms is avoided. In this paper, we develop a multigrid algorithm for solving the resulting system of three coupled nonlinear equations. We include a discussion of basic multigrid methods, along with some more advanced techniques required for the development of the multigrid scheme for this system. Results are presented showing that usual optimal multigrid efficiency is obtained. [ABSTRACT FROM AUTHOR]

Dynamical systems can often be decomposed into loosely coupled subsystems. The system of ordinary differential equations (ODEs) modelling such a problem can then be partitioned corresponding to the subsystems, and the loose couplings can be exploited by special integration methods to solve the problem using a parallel computer or just solve the problem more efficiently than by standard methods. This paper presents accuracy analysis of methods for the numerical integration of stiff partitioned systems of ODEs. The discretization formulas are based on the implicit Euler formula and the second order implicit backward differentiation formula (BDF2). Each subsystem of the partitioned problem is discretized independently, and the couplings to the other subsystems are based on solution values from previous time steps. Applied this way, the discretization formulas are called decoupled. The stability properties of the decoupled implicit Euler formula are well understood. This paper presents error bounds and asymptotic error expansions to be used in controlling step size, relaxation between subsystems and the validity of the partitioning. The decoupled BDF2 formula is analyzed within the same framework. Finally, the analysis is used in the design of a decoupled numerical integration algorithm with variable step size to control the local error and adaptive selection of partitionings. Two versions of the algorithm with decoupled implicit Euler and BDF2, respectively, are used in examples where a realistic problem is solved. The examples compare the results from the decoupled implicit Euler and BDF2 formulas, and compare them with results from the corresponding classical formulas. [ABSTRACT FROM AUTHOR]

Stochastic models of chemical systems are often analyzed by solving the corresponding Fokker--Planck equation, which is a drift-diffusion partial differential equation for the probability distribution function. Efficient numerical solution of the Fokker--Planck equation requires adaptive mesh refinements. In this paper, we present a mesh refinement approach which makes use of a stochastic simulation of the underlying chemical system. By observing the stochastic trajectory for a relatively short amount of time, the areas of the state space with nonnegligible probability density are identified. By refining the finite element mesh in these areas, and coarsening elsewhere, a suitable mesh is constructed and used for the computation of the stationary probability density. Numerical examples demonstrate that the presented method is competitive with existing a posteriori methods. [ABSTRACT FROM AUTHOR]

Mesh generation by Delaunay refinement is a widely used technique for constructing guaranteed quality triangular and tetrahedral meshes. The quality guarantees are usually provided in terms of the bounds on circumradius-to-shortest-edge ratio and on the grading of the resulting mesh. Traditionally circumcenters of skinny elements and middle points of boundary faces and edges are used for the positions of inserted points. However, recently variations of the traditional algorithms are being proposed that are designed to achieve certain optimization objectives by inserting new points in neighborhoods of the center points. In this paper we propose a general approach to the selection of point positions by defining one-, two-, and three-dimensional selection regions such that any point insertion strategy based on these regions is automatically endowed with the theoretical guarantees proven here. In particular, for the input models defined by planar linear complexes under the assumption that no input angle is less than 90°, we prove the termination of the proposed generalized algorithm, as well as the fidelity and good grading of the resulting meshes. [ABSTRACT FROM AUTHOR]

Modeling atomic and molecular systems requires computation-intensive quantum mechanical methods such as, but not limited to, density functional theory [R. A. Friesner, Proc. Natl. Acad. Sci. USA, 102 (2005), pp. 6648-6653]. These methods have been successful in predicting various properties of chemical systems at atomistic scales. Due to the inherent nonlocality of quantum mechanics, the scalability of these methods ranges from O(N³)to O (N7) depending on the method used and approximations involved. This significantly limits the size of simulated systems to a few thousand atoms, even on large scale parallel platforms. On the other hand, classical approximations of quantum systems, although computationally (relatively) easy to implement, yield simpler models that lack essential chemical properties such as reactivity and charge transfer. The recent work of van Duin et al. [J. Phys. Chem. A, 105 (2001), pp. 9396-9409] overcomes the limitations of nonreactive classical molecular dynamics (MD) approximations by carefully incorporating limited nonlocality (to mimic quantum behavior) through an empirical bond order potential. This reactive classical MD method, called ReaxFF, achieves essential quantum properties, while retaining the computational simplicity of classical MD, to a large extent. Implementation of reactive force fields presents significant algorithmic challenges. Since these methods model bond breaking and formation, efficient implementations must rely on complex dynamic data structures. Charge transfer in these methods is accomplished by minimizing electrostatic energy through charge equilibration. This requires the solution of large linear systems (108 degrees of freedom and beyond) with shielded electrostatic kernels at each time-step. Individual time-steps are themselves typically in the range of tenths of femtoseconds, requiring optimizations within and across time-steps to scale simulations to nanoseconds and beyond, where interesting phenomena may be observed. In this paper, we present implementation details of sPuReMD (serial Purdue reactive molecular dynamics program), a unique reactive classical MD code. We describe various data structures, and the charge equilibration solver at the core of the simulation engine. This Krylov subspace solver relies on a preconditioner based on incomplete LU factorization with thresholds (ILUT), specially targeted to our application. We comprehensively validate the performance and accuracy of sPuReMD on a variety of hydrocarbon systems. In particular, we show excellent per-time-step time, linear time scaling in system size, and a low memory footprint. sPuReMD is a freely distributed software with GPL and is currently being used to model diverse systems ranging from oxidative stress in biomembranes to strain relaxation in Si-Ge nanorods. [ABSTRACT FROM AUTHOR]

Anisotropic diffusion methods are well known for giving good qualitative results for image denoising. This paper gives a review of the anisotropic diffusion methodology and its application to image denoising. We propose a fixed-point iteration using a multigrid solver to solve a regularized anisotropic diffusion equation, which is not only well-posed, but also has a nontrivial steady-state solution. A new regularization parameter-choice method (Brent-NCP), combining Brent's method and the normalized cumulative periodogram information of the misfit, is also introduced. We test our algorithm on several common test images with different noise levels. The experimental results demonstrate the effectiveness of the anisotropic diffusion with a multigrid approach and the broad applicability of the Brent-NCP parameter-choice algorithm. [ABSTRACT FROM AUTHOR]

We present two iterative algorithms for solving real nonsymmetric and complex nonHermitian linear systems of equations and that were developed from variants of the nonsymmetric Lanczos method. In this paper, we give the theoretical background of the two iterative methods and discuss their main computational aspects. Using a large number of numerical experiments, we analyze their convergence properties, and we also compare them with other popular nonsymmetric iterative solvers in use today. [ABSTRACT FROM AUTHOR]

Linear advection of a scalar quantity by a specified velocity field arises in a number of different applications. Of particular interest here is the transport of species and energy in low Mach number models for combustion, atmospheric flows, and astrophysics, as well as contaminant transport in Darcy models of saturated subsurface flow. An important characteristic of these problems is that the velocity field is not known analytically. Instead, an auxiliary equation is solved to compute averages of the velocities over faces in a finite volume discretization. In this paper, we present a customized three-dimensional finite volume advection scheme for this class of problems that provides accurate resolution for smooth problems while avoiding undershoot and overshoot for nonsmooth profiles. The method is an extension of an algorithm by Bell, Dawson, and Shubin (BDS), which was developed for a class of scalar conservation laws arising in porous media flows in two dimensions. The original BDS algorithm is a variant of unsplit. higher-order Godunov methods based on construction of a limited bilinear profile within each computational cell. Here we present a three-dimensional extension of the original BDS algorithm that is based on a limited trilinear profile within each cell. We compare this new method to several other unsplit approaches, including piecewise linear methods, piecewise parabolic methods, and wave propagation schemes. [ABSTRACT FROM AUTHOR]

This paper presents an efficient multiscale stochastic framework for uncertainty quantification in modeling of flow through porous media with multiple rock types. The governing equations are based on Darcy's law with nonstationary stochastic permeability represented as a sum of local Karhunen-Loève expansions. The approximation uses stochastic collocation on either a tensor product or a sparse grid, coupled with a domain decomposition algorithm known as the multiscale mortar mixed finite element method. The latter method requires solving a coarse scale mortar interface problem via an iterative procedure. The traditional implementation requires the solution of local fine scale linear systems on each iteration. We employ a recently developed modification of this method that precomputes a multiscale flux basis to avoid the need for subdomain solves on each iteration. In the stochastic setting, the basis is further reused over multiple realizations, leading to collocation algorithms that are more efficient than the traditional implementation by orders of magnitude. Error analysis and numerical experiments are presented. [ABSTRACT FROM AUTHOR]

In this paper we show that recently introduced quantics tensor train (QTT) decomposition call be considered as an algebraic wavelet transform with adaptively determined filters. The main algorithm for obtaining QTT decomposition can be reformulated as a method seeking "good subspaces" or "good bases" and considered as a parameterized transformation of an initial tensor into a sparse tensor. This interpretation allows us to introduce a modification of the tensor train-SVD (TT-SVD) algorithm to make it work in cases where the original algorithm does not work; it results in the new wavelet-like transforms called wavelet tensor train (WTT) transform. Properties of WTT transforms are studied numerically, and a theoretical conjecture on the number of vanishing moments is proposed. It is shown that WTT transforms are orthogonal by construction, and the efficiency of WTT is compared with and often outperforms Daubechies wavelet transforms on certain classes of function-related vectors and matrices. [ABSTRACT FROM AUTHOR]

This paper is concerned with the numerical simulations for the dynamics of the molecular beam epitaxy (MBE) model. The numerical simulations of the MBE models require long time computations, and therefore large time-stepping methods become necessary. In this work, we consider some unconditionally energy stable finite difference schemes, which will be used in the time adaptivity strategies. It is found that the use of the time adaptivity cannot only resolve the steady-state solutions but also the dynamical changes of the solution accurately and efficiently. The adaptive time step is selected based on the energy variation or the change of the roughness of the solution. The numerical experiments demonstrated that the CPU time is significantly saved for long time simulations. [ABSTRACT FROM AUTHOR]

High order accurate weighted essentially nonoscillatory (WENO) schemes are usually designed to solve hyperbolic conservation laws or to discretize the first derivative convection terms in convection dominated partial differential equations. In this paper we discuss a high order WENO finite difference discretization for nonlinear degenerate parabolic equations which may contain discontinuous solutions. A porous medium equation (PME) is used as all example to demonstrate the algorithm structure and performance. By directly approximating the second derivative term using a conservative flux difference, the sixth order and eighth order finite difference WENO schemes are constructed. Numerical examples are provided to demonstrate the accuracy and nonoscillatory performance of these schemes. [ABSTRACT FROM AUTHOR]

A new algorithm is developed for computing etAB, where A is an n x n matrix and B is n x n0 with n0 « n. The algorithm works for any A, its computational cost is dominated by the formation of products of A with n x n0 matrices, and the only input parameter is a backward error tolerance. The algorithm can return a single matrix etA B or a sequence Due to image rights restrictions, multiple line equation(s) cannot be graphically displayed. on an equally spaced grid of points tk. It uses the scaling part of the scaling and squaring method together with a truncated Taylor series approximation to the exponential. It determines the amount of scaling and the Taylor degree, using the recent analysis of Al-Mohy and Higham [SIAM J. Matrix Anal. Appl., 31 (2009), pp. 970 989], which provides sharp truncation error bounds expressed in terms of the quantities ||Ak

This paper presents a numerical method to simulate an incompressible fluid through an immersed porous interface. The interface is modeled by a surface measure term in the Navier-Stokes equations and it is characterized by a resistance parameter. This approach can be used, for example, to model valves or to simulate blood flood through an immersed stent. Starting from a monolithic formulation proposed recently, a fractional step algorithm is derived. The difficult point is that this formulation is singular when the resistance vanishes, which can be a serious issue in some applications. We show that an appropriate Nitsche treatment of the interface condition fixes this problem and ensures uniform energy stability in time for any nonnegative value of the resistance. The theoretical stability and convergence results are illustrated with numerical experiments. [ABSTRACT FROM AUTHOR]

In this paper, we propose and study the use of alternating direction algorithms for several ℓ1-norm minimization problems arising from sparse solution recovery in compressive sensing, including the basis pursuit problem, the basis pursuit denoising problems of both unconstrained and constrained forms, and others. We present and investigate two classes of algorithms derived from either the primal or the dual form of ℓ1-problems. The construction of the algorithms consists of two main steps: (1) to reformulate an ℓ1-problem into one having blockwise separable objective functions by adding new variables and constraints; and (2) to apply an exact or inexact alternating direction method to the augmented Lagrangian function of the resulting problem. The derived alternating direction algorithms can be regarded as first-order primal-dual algorithms because both primal and dual variables are updated at every iteration. Convergence properties of these algorithms are established or restated when they already exist. Extensive numerical experiments are performed, using randomized partial Walsh Hadamard sensing matrices, to demonstrate the versatility and effectiveness of the proposed approach. Moreover, we present numerical results to emphasize two practically important but perhaps overlooked points: (i) that algorithm speed should be evaluated relative to appropriate solution accuracy; and (ii) that when erroneous measurements possibly exist, the ℓ1-fidelity should generally be preferable to the ℓ2-fidelity. [ABSTRACT FROM AUTHOR]

Numerical algorithms have two kinds of costs: arithmetic and communication, by which we mean either moving data between levels of a memory hierarchy (in the sequential case) or over a network connecting processors (in the parallel case). Communication costs often dominate arithmetic costs, so it is of interest to design algorithms minimizing communication. In this paper we first extend known lower bounds on the communication cost (both for bandwidth and for latency) of conventional (O(n3)) matrix multiplication to Cholesky factorization, which is used for solving dense symmetric positive definite linear systems. Second, we compare the costs of various Cholesky decomposition implementations to these lower bounds and identify the algorithms and data structures that attain them. In the sequential case, we consider both the two-level and hierarchical memory models. Combined with prior results in [J. Demmel et al., Communication-optimal Parallel and Sequential QR and LU Factorizations, Technical report EECS-2008-89, University of California, Berkeley, CA, 2008], [J. Demmel et al., Implementing Communication-optimal Parallel and Sequential QR and LU Factorizations, SIAM. J. Sci. Comp., submitted], and [J. Demmel, L. Grigori, and H. Xiang, Communication-avoiding Gaussian Elimination, Proceedings of the 2008 ACM/IEEE Conference on Supercomputing, 2008] this gives a set of communication-optimal algorithms for O(n³) implementations of the three basic factorizations of dense linear algebra: LU with pivoting, QR, and Cholesky. But it goes beyond this prior work on sequential LU by optimizing communication for any number of levels of memory hierarchy. [ABSTRACT FROM AUTHOR]

In this paper we proposed quasi-Newton and limited memory quasi-Newton methods for objective functions defined ell Grassmannians or a product of Grassmannians. Specifically we defined BFGS and limited memory BFGS updates in local and global coordinates oil Grassmannians or a product, of these. We proved that, when local coordinates are used, our BFGS updates on Grassmaunians share the same optimality property as the usual BFGS updates on Euclidean spaces. When applied to the best multilinear rank approximation problem for general and symmetric tensors, our approach yields fast, robust, and accurate algorithms that exploit the special Grassmannian structure of the respective problems and which work on tensors of large dimensions and arbitrarily high order. Extensive numerical experiments are included to substantiate our claims. [ABSTRACT FROM AUTHOR]

We develop in this paper some efficient algorithms which are essential to implementations of spectral methods on the sparse grid by Smolyak's construction based on a nested quadrature. More precisely, we develop a fast algorithm for the discrete transform between the values at the sparse grid and the coefficients of expansion in a hierarchical basis; and by using the aforementioned fast transform, we construct two very efficient sparse spectral-Galerkin methods for a model elliptic equation. In particular, the Chebyshev Legendre-Galerkin method leads to a sparse matrix with a much lower number of nonzero elements than that of low-order sparse grid methods based on finite elements or wavelets, and can be efficiently solved by a suitable sparse solver. Ample numerical results are presented to demonstrate the efficiency and accuracy of our algorithms. [ABSTRACT FROM AUTHOR]

Parallel simulations at extreme scale require that the mesh is distributed across a large number of processors with equal work load and minimum interpart communications. A number of algorithms have been developed to meet these goals, e.g., graph/hypergraph and coordinate- based methods. However, the global implementation of current approaches Call fail oil very large core counts, which is resolved by combining global and local partitioning using multiple parts per processor. The other limitation of graph/hypergraph-based partitioning is that it uses one type of mesh entity as graph nodes; thus, the balance of other mesh entities may not be optimal. In the case of three-dimensional (3-D) linear finite element analysis, it is common to select mesh regions (elements) as partition objects. In current examples, the regions axe well balanced up to 163,840 parts for a 1.07 billion element mesh, while the vertices have an imbalance which is as high as 19.52%. Two methods are developed that work in conjunction with graph/hypergraph-based procedures to provide improved partitions. Example computations executed on an IBM Blue Gene/P system using up to 163,840 cores demonstrate the usefulness of the procedures, particularly for time-critical calculations where individual cores may be lightly loaded in terms of the number of mesh entities per core. The algorithms presented in this paper reduced the vertex imbalance from 17.8% to 4.97% for a partition with 131,072 parts and accelerated the equation solution phase of the finite element analysis by 10.4%. [ABSTRACT FROM AUTHOR]

Sparse direct solvers, and in particular multifrontal methods, are methods of choice to solve the large sparse systems of linear equations arising in certain simulation problems. However, they require a large amount of memory (e.g., in comparison to iterative methods). In this context, out-of-core solvers may be employed: disks are used when the required storage exceeds the available physical memory. In this paper, we show how to process the task dependency graph of multifrontal methods in a way that minimizes the input/output (I/O) requirements. From a theoretical point of view, we show that minimizing the storage requirement can lead to a huge volume of I/O compared to directly minimizing the I/O volume. Then experiments on large real-world problems also show that applying standard algorithms to minimize the storage is not always efficient at reducing the volume of I/O and that significant gains can be obtained with the use of our algorithms to minimize I/O. We finally show that efficient memory management algorithms can be applied to all the variants proposed. [ABSTRACT FROM AUTHOR]

One approach to solve optimal control problems by direct methods is the so-called sequential approach or single shooting. Only the control variables are discretized resulting in a nonlinear program (NLP) which can be solved with SQP or interior point methods. This paper presents a new methodology to efficiently provide the Hessian of the Lagrangian of that resulting NLP. The algorithm is based on the second-order adjoint method and introduces the novel concept of composite adjoints to reduce the computational effort of a Hessian evaluation. Though this contribution is for the sake of simplicity restricted to single shooting, the same methodology can also be easily applied to multiple shooting. [ABSTRACT FROM AUTHOR]

In this paper we derive two a posteriori upper bounds for the heat equation. A continuous, piecewise linear finite element discretization in space and the Crank-Nicolson method for the time discretization are used. The error due to the space discretization is derived using anisotropic interpolation estimates and a postprocessing procedure. The error due to the time discretization is obtained using two different continuous, piecewise quadratic time reconstructions. The first reconstruction is developed following G. Akrivis, C. Makridakis, and R. H. Nochetto [Math. Comp., 75 (2006), pp. 511-531], while the second one is new. Moreover, in the case of isotropic meshes only, upper and lower bounds are provided as in [R. Verfürth, Calcolo, 40 (2003), pp. 195- 212]. An adaptive algorithm is developed. Numerical studies are reported for several test cases and show that the second error estimator is more efficient than the first one. In particular, the second error indicator is of optimal order with respect to both the mesh size and the time step when using our adaptive algorithm. [ABSTRACT FROM AUTHOR]

This paper presents an algorithm for finding a piecewise linear band which bounds given two-dimensional data. The band is of constant height, and the distance between adjacent slope discontinuities in opposite directions is restricted by two user-specified parameters. [ABSTRACT FROM AUTHOR]

Image restoration problems are often solved by finding the minimizer of a suitable objective function consisting of a data-fitting term and a regularization term. In this paper, we consider the data-fitting term measured in the ℓ1 norm to handle non-Gaussian additive noise and the regularization term given by the total variation (TV) to restore image edges. We propose a new algorithm for this image restoration problem by making use of new variables to modify the data-fitting term and the TV regularization term. An alternating minimization method based on the new formulation is employed to restore blurred and noisy images. Our experimental results show that the quality of restored images by the proposed method is competitive with those restored by the other tested methods. We also show the convergence of the alternating minimization algorithm and demonstrate that the proposed algorithm is very efficient. [ABSTRACT FROM AUTHOR]

In this paper, we propose a new phase-flow method for Hamiltonian systems with discontinuous Hamiltonians. In the original phase-flow method introduced by Ying and Candès [J. Comput. Phys., 220 (2006), pp. 184-215], the phase map should be smooth to ensure the accuracy of the interpolation. Such an interpolation is inaccurate if the phase map is nonsmooth, for example, when the Hamiltonian is discontinuous. We modify the phase-flow method using a discontinuous Hamiltonian solver and establish the stability (for piecewise constant potentials) of such a solver. This extends the applicability of the highly efficient phase-flow method to singular Hamiltonian systems, with a mild increase of algorithm complexity. Such a particle method can be useful for the computation of high frequency waves through interfaces. [ABSTRACT FROM AUTHOR]

In Part II of this paper we first refine the analysis of error-free vector transformations presented in Part I. Based on that we present an algorithm for calculating the rounded-to-nearest result of s := Σ pi for a given vector of floating-point numbers pi, as well as algorithms for directed rounding. A special algorithm for computing the sign of s is given, also working for huge dimensions. Assume a floating-point working precision with relative rounding error unit eps. We define and investigate a K-fold faithful rounding of a real number r. Basically the result is stored in a vector Resν of K nonoverlapping floating-point numbers such that Σ Resν approximates r with relative accuracy epsK, and replacing ResK by its floating-point neighbors in Σ Resν forms a lower and upper bound for r. For a given vector of floating-point numbers with exact sum s, we present an algorithm for calculating a K-fold faithful rounding of s using solely the working precision. Furthermore, an algorithm for calculating a faithfully rounded result of the sum of a vector of huge dimension is presented. Our algorithms are fast in terms of measured computing time because they allow good instruction-level parallelism, they neither require special operations such as access to mantissa or exponent, they contain no branch in the inner loop, nor do they require some extra precision. The only operations used are standard floating-point addition, subtraction, and multiplication in one working precision, for example, double precision. Certain constants used in the algorithms are proved to be optimal. [ABSTRACT FROM AUTHOR]

In this paper, we propose a novel multiresolution adaptive mesh refinement algorithm for solving initial-boundary value problems (IBVP) for evolution PDEs. The proposed algorithm dynamically adapts the grid to any existing or emerging irregularities in the solution, by refining the grid only at those places where the solution exhibits sharp features. The main advantage of the proposed grid adaptation method is that it results in a grid with a fewer number of nodes when compared to adaptive grids generated by existing multiresolution-based mesh refinement techniques. Several examples show the robustness and stability of the proposed algorithm. [ABSTRACT FROM AUTHOR]

In this paper, we present a multigrid method and its implementation on parallel computers to solve the bidomain equation that appears in excitation propagation analysis of the human heart with the torso. The bidomain equation is discretized with the finite element method on a composite mesh composed of a fine voxel mesh around the heart and a coarse voxel mesh covering the torso. The extracellular potential problem on the torso is formulated as a variational problem with a constraint at the interface of the fine and coarse meshes. We show that this formulation naturally satisfies the conservation property of the electric currents and fits into the multilevel adaptive solution technique framework. We also present our special treatment of the Purkinje fiber network in the multigrid algorithm where it is modeled as multiway branching lines connected to the nodes in the voxel mesh of the heart. A parallel implementation of the proposed multigrid algorithm on distributed memory computers is presented and the performance is evaluated using real-life applications. [ABSTRACT FROM AUTHOR]

In this paper we present a new incomplete factorization of a square matrix into triangular factors in which we get standard LU or LDLT factors (direct factors) and their inverses (inverse factors) at the same time. Algorithmically, we derive this method from the approach based on the Sherman-Morrison formula [R. Bru, J. Cerdán, J. Marín, and J. Mas, SIAM J. Sci. Comput., 25 (2003), pp. 701-715]. In contrast to the robust incomplete decomposition (RIF) algorithm [M. Benzi and M. Tůma, Numer. Linear Algebra Appl., 10 (2003), pp. 385-400] the direct and inverse factors here directly influence each other throughout the computation. Consequently, the algorithm to compute the approximate factors may mutually balance dropping in the factors and control their conditioning in this way. For the symmetric positive definite case, we derive the theory and present an algorithm for computing the incomplete LDLT factorization, and we discuss experimental results. We call this new approximate LDLT factorization the balanced incomplete factorization (BIF). Our experimental results confirm that this factorization is very robust and may be useful in solving difficult ill conditioned problems by preconditioned iterative methods. Moreover, the internal coupling of the computation of direct and inverse factors results in much shorter setup times (times to compute approximate decomposition) than RIF, a method of a similar and very high level of robustness. We also derive and present the theory for the general nonsymmetric case, but do not discuss its implementation. [ABSTRACT FROM AUTHOR]

An effective 2-phase method for removing impulse noise was recently proposed. Its phase 1 identifies noisy pixel candidates by using median-type filters. Then in phase 2, it restores only the noisy pixel candidates by some variational methods. The resulting method can handle saltand- pepper noise and random-valued impulse noise at a level as high as 90% and 60%, respectively. The aim of this paper is to generalize a fast multilevel method for Gaussian denoising to solving the minimization problem arising in phase 2 of the 2-phase method. The multilevel algorithm gives better images than standard optimization methods such as the Newton method or the conjugate gradient method. Also it can handle more general regularization functionals than the smooth ones previously considered. Supporting numerical experiments on two-dimensional gray-scale images are presented. [ABSTRACT FROM AUTHOR]

Standing (solitary)-wave/steady-state solutions in many nonlinear wave motions and Schrödinger flows lead to nonlinear eigenproblems. In [X. Yao and J. Zhou, SIAM J. Sci. Comput., 29 (2007), pp. 1355-1374], a Rayleigh-local minimax method is developed to solve iso-homogeneous eigenproblems. In this subsequent paper, a unified method in Banach spaces is developed for solving non-iso-homogeneous and even nonhomogeneous eigenproblems and applied to solve two models: the Gross-Pitaevskii problem in the Bose-Einstein condensate and the p-Laplacian problem in non-Newtonian flows/materials. First a new active Lagrange functional is formulated to establish a local minimax characterization. A local minimax method is then devised and implemented to solve the model problems. Numerical results are presented. Convergence results of the algorithm and an order of eigensolutions computed by the algorithm are also established. [ABSTRACT FROM AUTHOR]

In this paper, we investigate the effects of pollution error on the performance of the parallel adaptive finite element technique proposed by Bank and Holst in 2000 [R. E. Bank and M. Holst, SIAM J. Sci. Comput., 22 (2000), pp. 1411-1443]. In particular, we consider whether the performance of the algorithm as it was originally proposed can be improved through the use of certain dual functions which give an indication of the global influences on subdomain errors. [ABSTRACT FROM AUTHOR]

In this paper a new stochastic algorithm for the solution of population balance equations is presented. The population balance equations have the form of extended Smoluchowski equations which include linear and source terms. The new algorithm, called the linear process deferment algorithm (LPDA), is used for solving a detailed model describing the formation of soot in premixed laminar flames. A measure theoretic formulation of a stochastic jump process is developed and the corresponding generator presented. The numerical properties of the algorithm are studied in detail and compared to the direct simulation algorithm and various splitting techniques. LPDA is designed for all kinds of population balance problems where nonlinear processes cannot be neglected but are dominated in rate by linear ones. In those cases the LPDA is seen to reduce run times for a population balance simulation by a factor of up to 1000 with a negligible loss of accuracy. [ABSTRACT FROM AUTHOR]

We propose a new algebraic multigrid (AMG) method for solving the eddy current approximations to Maxwell's equations. This AMG method has its roots in an algorithm proposed by Reitzinger and Schöberl. The main focus in the Reitzinger and Schöberl method is to maintain null-space properties of the weak ∇ × ∇ × operator on coarse grids. While these null-space properties are critical, they are not enough to guarantee h-independent convergence rates of the overall multigrid scheme. We present a new strategy for choosing intergrid transfers that not only maintains the important null-space properties on coarse grids but also yields significantly improved multigrid convergence rates. This improvement is related to those we explored in a previous paper, but is fundamentally simpler, easier to compute, and performs better with respect to both multigrid operator complexity and convergence rates. The new strategy builds on ideas in smoothed aggregation to improve the approximation property of an existing interpolation operator. By carefully choosing the smoothing operators, we show how it is sometimes possible to achieve h-independent convergence rates with a modest increase in multigrid operator complexity. Though this ideal case is not always possible, the overall algorithm performs significantly better than the original scheme in both iterations and run time. Finally, the Reitzinger and Schöberl method, as well as our previous smoothed method, are shown to be special cases of this new algorithm. [ABSTRACT FROM AUTHOR]

This paper presents a multigrid method for solving variable coefficient Maxwell's equations. The novelty in this method is the use of interpolation operators that do not produce multilevel commutativity complexes that lead to multilevel exactness. Rather, the effects of multilevel exactness are built into the level equations themselves—on the finest level using a discrete T – V formulation and on the coarser grids through the Galerkin coarsening procedure of a T – V formulation. These built-in structures permit the levelwise use of an effective hybrid smoother on the curl-free near-nullspace components and permit the development of interpolation operators for handling the curl-free and divergence-free error components separately. The resulting block-diagonal interpolation operator does not satisfy multilevel commutativity but has good approximation properties for both of these error components. Applying operator-dependent interpolation for each of these error components leads to an effective multigrid scheme for variable coefficient Maxwell's equations, where multilevel commutativity-based methods can degrade. Numerical results are presented to verify the effectiveness of this new scheme. [ABSTRACT FROM AUTHOR]

In this paper we present optimization algorithms for image restoration based on the total variation (TV) minimization framework of Rudin, Osher, and Fatemi (ROF). Our approach formulates TV minimization as a second-order cone program which is then solved by interior-point algorithms that are efficient both in practice (using nested dissection and domain decomposition) and in theory (i.e., they obtain solutions in polynomial time). In addition to the original ROF minimization model, we show how to apply our approach to other TV models, including ones that are not solvable by PDE-based methods. Numerical results on a varied set of images are presented to illustrate the effectiveness of our approach. [ABSTRACT FROM AUTHOR]

In this paper two noniterative algorithms for locating small electromagnetic inclusions from the scattering amplitude at a fixed frequency are developed. In particular, a variety of numerical results is presented to highlight their potential and their limitations. [ABSTRACT FROM AUTHOR]

This paper is concerned with the solution of a linear system with multiple right-hand sides. Such problems arise from nonlinear, time-dependent, inverse, or optimization problems. In order to solve these problems efficiently we use block variants of the conjugate gradient method and combine them with cache aware algorithms. Further, we describe the acceleration of algebraic multigrid (AMG) methods by using such cache aware algorithms. Finally, numerical studies are presented that show the high efficiency of the proposed methods. [ABSTRACT FROM AUTHOR]

Higher order panel methods are used to solve the Laplace equation in the presence of complex geometries. These methods are useful when globally accurate velocity or potential fields are desired as in the case of vortex based fluid flow solvers. This paper develops a fast multipole algorithm to compute velocity fields due to higher order, two-dimensional vortex panels. The technique is applied to panels having a cubic geometry and a linear distribution of vorticity. The results of the present method are compared with other available techniques. [ABSTRACT FROM AUTHOR]

In this work we study a dynamical system with a complicated nonlinearity, which describes oscillation of a turbine rotor, and give an algorithm to compute Melnikov functions for analysis of its chaotic behavior. We first derive the rotor model whose nonlinear term brings difficulties to investigating the distribution and qualitative properties of its equilibria. This nonlinear model provides a typical example of a system for which the homoclinic and heteroclinic orbits cannot be analytically determined. In order to apply Melnikov's method to make clear the underlying conditions for chaotic motion, we present a generic algorithm that provides a systematic procedure to compute Melnikov functions numerically. Substantial analysis is done so that the numerical approximation precision at each phase of the computation can be guaranteed. Using the algorithm developed in this paper, it is straightforward to obtain a sufficient condition for chaotic motion under damping and periodic external excitation, whenever the rotor parameters are given. [ABSTRACT FROM AUTHOR]

Parameterization of unstructured surface meshes is of fundamental importance in many applications of digital geometry processing. Such parameterization approaches give rise to large and exceedingly ill-conditioned systems which are difficult or impossible to solve without the use of sophisticated multilevel preconditioning strategies. Since the underlying meshes are very fine to begin with, such multilevel preconditioners require mesh coarsening to build an appropriate hierarchy. In this paper we consider several strategies for the construction of hierarchies using ideas from mesh simplification algorithms used in the computer graphics literature. We introduce two novel hierarchy construction schemes and demonstrate their superior performance when used in conjunction with a multigrid preconditioner. [ABSTRACT FROM AUTHOR]

In this paper we investigate the accuracy of two numerical procedures commonly used to solve 2D advection problems: spectral and finite difference (FD) schemes. These schemes are widely used, simulating, e.g., neutral and plasma flows. FD schemes have long been considered fast, relatively easy to implement, and applicable to complex geometries, but are somewhat inferior in accuracy compared to spectral schemes. Using two study cases at high Reynolds number, the merging of two equally signed Gaussian vortices in a periodic box and dipole interaction with a no-slip wall, we will demonstrate that the accuracy of FD schemes can be significantly improved if one is careful in choosing an appropriate FD scheme that reflects conservation properties of the nonlinear terms and in setting up the grid in accordance with the problem. [ABSTRACT FROM AUTHOR]

In this paper, we study some efficient numerical methods for parameter identifications in elliptic systems. The proposed numerical methods are conducted iteratively and each iteration involves only solving positive definite linear algebraic systems, although the original inverse problems are ill-posed and highly nonlinear. The positive definite systems can be naturally preconditioned with their corresponding block diagonal matrices. Numerical experiments are presented to illustrate the efficiency of the proposed algorithms. [ABSTRACT FROM AUTHOR]

This paper presents the results of a scalability study for a three-dimensional semi- coarsening multigrid solver on a distributed memory computer. In particular, we are interested in the scalability of the solver-how the solution time varies as both problem size and number of processors are increased. For an iterative linear solver, scalability involves both algorithmic issues and implementation issues. We examine the scalability of the solver theoretically by constructing a simple parallel model and experimentally by results obtained on an IBM SP. The results are compared with those obtained for other solvers on the same computer. [ABSTRACT FROM AUTHOR]

This paper formulates several algorithms for the direct computation of periodic orbits as solutions of boundary value problems. The algorithms emphasize the use of coarse meshes and high orders of accuracy. Convergence theorems are given in the limit of increasing order with a fixed mesh. The algorithms are implemented with the use of MATLAB and ADOL-C, a software package for automatic differentiation. Automatic differentiation enables accurate computation of high-order derivatives of functions without the truncation errors inherent in finite difference calculations. We embed the algorithms in a continuation framework and extend them to compute saddle-node bifur- cations of periodic orbits directly. We present data from numerical studies of four test problems, making some comparisons with other methods for computing periodic orbits. These results demonstrate that high-order methods based upon automatic differentiation are capable of high precision with small meshes. [ABSTRACT FROM AUTHOR]