Search

This paper analyzes the process of merging groups of many Gaussian basis functions into a single basis function in vortex simulations. Analysis of the equations governing this process yields fundamental parameters and uniform estimates of the merger-induced error. In a merging event, the uniform error bound depends only on relationships between nearby computational elements permitting fast and effective merging schemes. In this paper, one such algorithm is proposed and demonstrated. [ABSTRACT FROM AUTHOR]

This paper presents some guidance for choosing multigrid preconditioners for systems of partial differential equations (PDEs). Rather than developing a black-box algorithm for arbitrary systems of PDEs, the strategy in this paper is to classify systems with multigrid preconditioners. From this classification, a new preconditioner is developed, and it is shown to be suitable for a particular class of strongly cross-coupled systems, solving the whole cross-coupled system almost as effectively as solving just the diagonal system. Also, this paper gives an explanation of why block diagonal preconditioners are somewhat effective for certain strongly cross-coupled systems that at first glance appear to be intractable for multigrid preconditioning. [ABSTRACT FROM AUTHOR]

In this paper, we propose iterative algorithms for solving image restoration problems. The iterative algorithms are based on decoupling of deblurring and denoising steps in the restoration process. In the deblurring step, an efficient deblurring method using fast transforms can be employed. In the denoising step, effective methods such as the wavelet shrinkage denoising method or the total variation denoising method can be used. The main advantage of this proposal is that the resulting algorithms can be very efficient and can produce better restored images in visual quality and signalto-noise ratio than those by the restoration methods using the combination of a data-fitting term and a regularization term. The convergence of the proposed algorithms is shown in the paper. Numerical examples are also given to demonstrate the effectiveness of these algorithms. [ABSTRACT FROM AUTHOR]

The Schwarz alternating method can be used to solve linear elliptic boundary value problems on domains which consist of two or more overlapping subdomains. The solution is approx- imated by an infinite sequence of functions which result from solving a sequence of elliptic boundary value problems in each of the subdomains. This paper considers four Schwarz alternating methods for the N-dimensional, steady, viscous, incompressible Navier­Stokes equations, N ≤ 4. It is shown that the Schwarz sequences converge to the true solution provided that the Reynolds number is sufficiently small. [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]

In this paper, we consider the unilateral frictional contact problem of a hyperelastic body in the case of large displacements and small strains. In order to retain the linear elasticity framework, we decompose the deformation into a large global rotation and a small elastic displacement. This corotational approach is combined with a primal-dual active set strategy to tackle the contact problem. The resulting algorithm preserves both energy and angular momentum. [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]

We describe an asynchronous parallel derivative-free algorithm for linearly constrained optimization. Generating set search (GSS) is the basis of our method. At each iteration, a GSS algorithm computes a set of search directions and corresponding trial points and then evaluates the objective function value at each trial point. Asynchronous versions of the algorithm have been developed in the unconstrained and bound-constrained cases which allow the iterations to continue (and new trial points to be generated and evaluated) as soon as any other trial point completes. This enables better utilization of parallel resources and a reduction in overall run time, especially for problems where the objective function takes minutes or hours to compute. For linearly constrained GSS, the convergence theory requires that the set of search directions conforms to the nearby boundary. This creates an immediate obstacle for asynchronous methods where the definition of nearby is not well defined. In this paper, we develop an asynchronous linearly constrained GSS method that overcomes this difficulty and maintains the original convergence theory. We describe our implementation in detail, including how to avoid function evaluations by caching function values and using approximate lookups. We test our implementation on every CUTEr test problem with general linear constraints and up to 1000 variables. Without tuning to individual problems, our implementation was able to solve 95% of the test problems with 10 or fewer variables, 73% of the problems with 11-100 variables, and nearly half of the problems with 100-1000 variables. To the best of our knowledge, these are the best results that have ever been achieved with a derivative-free method for linearly constrained optimization. Our asynchronous parallel implementation is freely available as part of the APPSPACK software. [ABSTRACT FROM AUTHOR]

We discuss in this paper the numerical approximation of fluid-structure interaction (FSI) problems dealing with strong added-mass effect. We propose new semi-implicit algorithms based on inexact block-LU factorization of the linear system obtained after the space-time discretization and linearization of the FSI problem. As a result, the fluid velocity is computed separately from the coupled pressure-structure velocity system at each iteration, reducing the computational cost. We investigate explicit-implicit decomposition through algebraic splitting techniques originally designed for the FSI problem. This approach leads to two different families of methods which extend to FSI the algebraic pressure correction method and the Yosida method, two schemes that were previously adopted for pure fluid problems. Furthermore, we have considered the inexact factorization of the fluid-structure system as a preconditioner. The numerical properties of these methods have been tested on a model problem representing a blood-vessel system. [ABSTRACT FROM AUTHOR]

In this paper we describe a simple algorithm for the approximation of an invariant curve T of a planar map f. We are particularly interested in the case where T is attracting but not smooth. Results for the delayed logistic map illustrate the performance of the algorithm. [ABSTRACT FROM AUTHOR]

In this paper we present a new ILU decomposition which is based on an existing sparse direct solver. In contrast to many incomplete LU decompositions this ILU incorporates information about the inverse factors L-1 and U-1 which have direct influence on the dropping strategy. We demonstrate in several large scale examples that this implementation constructs a robust preconditioner. [ABSTRACT FROM AUTHOR]

In this paper, a discrete singular convolution (DSC) algorithm is introduced to solve Fisher’s equation, to which obtaining an accurate and reliable traveling wave solution is a challenging numerical problem. Two novel numerical treatments, a moving frame scheme and a new asymptotic scheme, are designed to overcome the subtle difficulties involved in the numerical solution of Fisher’s equation. The moving frame scheme is proposed to provide accurate and efficient spatial resolution of the problem. The new asymptotic scheme is introduced to facilitate the correct prediction of the wave speed after long-time integrations. The resulting DSC algorithm is able to correctly predict long-time traveling wave behavior. The performance of the present algorithm is demonstrated through extensive numerical experiments as well as through a comparison with three other standard approaches, i.e., methods of the Fourier pseudospectral, the accurate spatial derivatives, and the Crank­Nicolson schemes. The spatial and temporal accuracies of all four methods are analyzed and numerically verified. [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]

The use of approximate projection methods for modeling low Mach number flows avoids many of the numerical complications associated with exact projection methods, but introduces additional design choices in developing a robust algorithm. In this paper we first explore these design choices in the setting of inviscid incompressible flow using several computational examples. We then develop a framework for analyzing the behavior of the different design variations and use that analysis to explain the features observed in the computations. As part of this work we introduce a new variation of the approximate projection algorithm that combines the advantages of several of the previous versions. [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]

In this paper we consider a matrix-free path following algorithm for nonlinear parameter-dependent compact fixed point problems. We show that if these problems are discretized so that certain collective compactness and strong convergence properties hold, then this algorithm can follow smooth folds and capture simple bifurcations in a mesh-independent way. [ABSTRACT FROM AUTHOR]

Multilevel algorithms are a successful class of optimization techniques which addresses the mesh partitioning problem. They usually combine a graph contraction algorithm together with a local optimization method which refines the partition at each graph level. In this paper we present an enhancement of the technique which uses imbalance to achieve higher quality partitions. We also present a formulation of the Kernighan­Lin partition optimization algorithm which incorporates load-balancing. The resulting algorithm is tested against a different but related state-of-the-art partitioner and shown to provide improved results. [ABSTRACT FROM AUTHOR]

A successive continuation method for locating connecting orbits in parametrized systems of autonomous ODEs was considered in [Numer. Algorithms, 14 (1997), pp. 103­124]. In this paper we present an improved algorithm for locating and continuing connecting orbits, which includes a new algorithm for the continuation of invariant subspaces. The latter algorithm is of independent interest and can be used in contexts different than the present one. [ABSTRACT FROM AUTHOR]

Low-rank approximation of large and/or sparse matrices is important in many applications, and the singular value decomposition (SVD) gives the best low-rank approximations with respect to unitarily-invariant norms. In this paper we show that good low-rank approximations can be directly obtained from the Lanczos bidiagonalization process applied to the given matrix without computing any SVD. We also demonstrate that a so-called one-sided reorthogonalization process can be used to maintain an adequate level of orthogonality among the Lanczos vectors and produce accurate low-rank approximations. This technique reduces the computational cost of the Lanczos bidiagonalization process. We illustrate the efficiency and applicability of our algorithm using numerical examples from several applications areas. [ABSTRACT FROM AUTHOR]

Multidimensional positive definite advection transport algorithm (MPDATA) is an iterative process for approximating the advection equation, which uses a donor cell approximation to compensate for the truncation error of the originally specified donor cell scheme. This step may be repeated an arbitrary number of times, leading to successively more accurate solutions to the advection equation. In this paper, we show how to sum the successive approximations analytically to find a single antidiffusive velocity that represents the effects of an arbitrary number of passes. The analysis is first done in one dimension to illustrate the method and then is repeated in two dimensions. The existence of cross terms in the truncation analysis of the two-dimensional equations introduces an extra complication into the calculation. We discuss the implementation of our new antidiffusive velocities and provide some examples of applications, including a third-order accurate scheme. [ABSTRACT FROM AUTHOR]

In this paper we present several new fast and reliable algorithms for solving rank-deficient linear least squares problems by means of the complete orthogonal decomposition. Experimental comparison of our algorithms with the existing implementations in LAPACK on a wide range of platforms shows considerable performance improvements. Moreover, for full-rank matrices the performance of the new algorithm is very close to that of the fast method based on QR factorization, thus providing a valuable general tool for full-rank matrices, rank-deficient matrices, and those matrices with unknown rank. [ABSTRACT FROM AUTHOR]

In this paper we describe and analyze a class of nonlinear multiresolution schemes for the multiresolution setting which corresponds to discretization by local averages with respect to the hat function. These schemes are based on the essentially-non-oscillatory (ENO) interpolatory procedure described in [A. Harten, B. Engquist, S. Osher, and S. Chakravarthy, J. Comput. Phys., 71 (1987), pp. 231–302]. We show that by allowing the approximation to fit the local nature of the data, one can improve the compression capabilities of the multiresolution algorithms. The question of stability for nonlinear (data-dependent) reconstruction techniques is also addressed. [ABSTRACT FROM AUTHOR]

When performing sparse matrix factorization, the ordering of matrix rows and columns has a dramatic impact on the factorization time. This paper describes an approach to the reordering problem that produces significantly better orderings than prior methods. The algorithm is a hybrid of nested dissection and minimum degree ordering, and combines an assortment of different algorithmic advances. New or improved algorithms are described for graph compression, multilevel partitioning, and separator improvement. When these techniques are combined, the resulting orderings average 39% better than minimum degree over a suite of test matrices, while requiring roughly 2.7 times the run time of Liu's multiple minimum degree. [ABSTRACT FROM AUTHOR]

An adaptive remeshing procedure based on a cell volume deformation method is presented. Starting with an initial grid, this method offers direct cell volume control through the specification of the transformation Jacobian. Grid points are moved with appropriate grid velocities so that the specified cell volume distribution can be achieved at the end of the grid movement without adding or removing grid points. The grid velocities are determined by solving a scalar Poisson equation. This method is applied to solving the compressible Euler equations. Computational test cases of transonic flow over an airfoil are presented and demonstrate the desired control of grid sizes across shock waves. [ABSTRACT FROM AUTHOR]

A simple parallel algorithm for the evaluation of polynomials written in the Chebyshev form is introduced. By this method only 2 [log2(p - 2)] + 4 [log2p] + 4 [N/p] - 7 steps on p processors are needed to evaluate a Chebyshev series of degree N. Theoretical analysis of the efficiency is performed and some numerical examples on a CRAY T3D are shown. [ABSTRACT FROM AUTHOR]

Many problems in fluid dynamics require the representation of complicated internal or external boundaries of the flow. Here we present a method for calculating time-dependent incompressible inviscid flow which combines a projection method with a \Cartesian grid" approach for representing geometry. In this approach, the body is represented as an interface embedded in a regular Cartesian mesh. The advection step is based on a Cartesian grid algorithm for compressible flow, in which the discretization of the body near the flow uses a volume-of-fluid representation. A redistribution procedure is used to eliminate time-step restrictions due to small cells where the boundary intersects the mesh. The projection step uses an approximate projection based on a Cartesian grid method for potential flow. The method incorporates knowledge of the body through volume and area fractions along with certain other integrals over the mixed cells. Convergence results are given for the projection itself and for the time-dependent algorithm in two dimensions. The method is also demonstrated on flow past a half-cylinder with vortex shedding. [ABSTRACT FROM AUTHOR]

We develop a very simple algorithm that permits to construct compact, high order schemes for steady first order Hamilton-Jacobi equations. The algorithm relies on the blending of a first order scheme and a compact high order scheme. The blending is conducted in such a way that the scheme is formally high order accurate. A convergence proof without error estimate is given. We provide several numerical illustrations that demonstrate the effective accuracy of the scheme. The numerical examples use triangular unstructured meshes, but our method may be applied to other kind of meshes. Several implementation remarks are also given. [ABSTRACT FROM AUTHOR]

In this article, we propose new parallel algorithms for the construction and 2:1 balance refinement of large linear octrees on distributed memory machines. Such octrees are used in many problems in computational science and engineering, e.g., object representation, image analysis, unstructured meshing, finite elements, adaptive mesh refinement, and N-body simulations. Fixed-size scalability and isogranular analysis of the algorithms using an MPI-based parallel implementation was performed on a variety of input data and demonstrated good scalability for different processor counts (1 to 1024 processors) on the Pittsburgh Supercomputing Center's TCS-1 AlphaServer. The results are consistent for different data distributions. Octrees with over a billion octants were constructed and balanced in less than a minute on 1024 processors. Like other existing algorithms for constructing and balancing octrees, our algorithms have O(N logN) work and O(N) storage complexity. Under reasonable assumptions on the distribution of octants and the work per octant, the parallel time complexity is O( N/np log( N/np ) + np log np), where N is the size of the final linear octree and np is the number of processors. [ABSTRACT FROM AUTHOR]

The VFRoe scheme has been recently introduced by Buffard, Gallouët, and Hérard [Comput. Fluids, 29 (2000), pp. 813-847] to approximate the solutions of the shallow water equations. One of the main interests of this method is to be easily implemented. As a consequence, such a scheme appears as an interesting alternative to other more sophisticated schemes. The VFRoe methods perform approximate solutions in good agreement with the expected ones. However, the robustness of this numerical procedure has not been proposed. Following the ideas introduced by Jin and Xin [Comm. Pure Appl. Math., 45 (1995), pp. 235-276], a relevant relaxation method is derived. The interest of this relaxation scheme is twofold. In the first hand, the relaxation scheme is shown to coincide with the considered VFRoe scheme. In the second hand, the robustness of the relaxation scheme is established, and thus the nonnegativity of the water height obtained involving the VFRoe approach is ensured. Following the same idea, a family of relaxation schemes is exhibited. Next, robust high order slope limiter methods, known as MUSCL reconstructions, are proposed. The final scheme is obtained when considering the hydrostatic reconstruction to approximate the topography source terms. Numerical experiments are performed to attest the interest of the procedure. [ABSTRACT FROM AUTHOR]

The difficulties encountered in sampling of systems with rough energy landscapes using present methodology significantly limit the impact of simulation on molecular biology, in particular protein folding and design. Here, we present a major methodological development based on a promising new technique, the reference potential spatial warping algorithm (REPSWA) [Z. Zhu et al., Phys. Rev. Lett., 88 (2002), pp. 100201-100204], and present applications to several realistic systems. REPSWA works by introducing a variable transformation in the classical partition function that reduces the volume of phase space associated with a priori known barrier regions while increasing that associated with attractive basins. In this way, the partition function is preserved so that enhanced sampling is achieved without the need for reweighting phase-space averages. Here, a new class of transformations, designed to overcome the barriers induced by intermolecular/nonbonded interactions, whose locations are not known a priori, is introduced. The new transformations are designed to work in synergy with transformations originally introduced for overcoming intramolecular barriers. The new transformation adapts to the fluctuating local environment and is able to handle barriers that arise "on the fly." Thus, the new method is referred to as dynamic contact REPSWA (DC-REPSWA). In addition, combining hybrid Monte Carlo (HMC) with DC-REPSWA allows more aggressive sampling to take place. The combined DC-REPSWA-HMC method and its variants are shown to substantially enhance conformational sampling in long molecular chains composed of interacting single beads and beads with branches. The latter topologies characterize the united residue and united side chain representation of protein structures. [ABSTRACT FROM AUTHOR]

A new algorithm for the computation of dominant poles of transfer functions of largescale second-order dynamical systems is presented: the quadratic dominant pole algorithm (QDPA). The algorithm works directly with the system matrices of the original system, so no linearization is needed. To improve global convergence, the QDPA uses subspace acceleration, and deflation of found dominant poles is implemented in a very efficient way. The dominant poles and corresponding eigenvectors can be used to construct structure-preserving modal approximations and also to improve reduced-order models computed by Krylov subspace methods, as is illustrated by numerical results. [ABSTRACT FROM AUTHOR]

We present a method for creating large two-dimensional graded Delaunay meshes on parallel machines. The method decouples the mesh generation procedure for each subdomain and thus eliminates the communication between processors, while the final mesh is of the same guaranteed quality as in the sequential case. This work extends our previous result on the uniform Delaunay decoupling method [L. Linardakis and N. Chrisochoides, SIAM J. Sci. Comput., 27 (2006), pp. 1394- 1423], by allowing the element size to be governed by a sizing function, or a background mesh. The graded Delaunay decoupling method can produce in parallel very large graded meshes efficiently, demonstrating high speedup and scalability. Ten billion elements can be created in five minutes with small overrefinement (about 2%). The decoupling parallelization approach is effective, using off-the-shelf, well-tested, and fine-tuned sequential mesh generation libraries without modification. [ABSTRACT FROM AUTHOR]

The problem of variational data assimilation for a nonlinear evolution model is formulated as an optimal control problem to find the initial condition function (analysis). The equation for the analysis error is derived through the errors of the input data (background and observation errors). This equation is used to show that in a nonlinear case the analysis error covariance operator can be approximated by the inverse Hessian of an auxiliary data assimilation problem which involves the tangent linear model constraints. The inverse Hessian is constructed by the quasi-Newton BFGS algorithm when solving the auxiliary data assimilation problem. A fully nonlinear ensemble procedure is developed to verify the accuracy of the proposed algorithm. Numerical examples are presented. [ABSTRACT FROM AUTHOR]

We present a two-phase algorithm for solving large-scale quadratic programs (QPs). In the first phase, gradient-projection iterations approximately minimize a bound-constrained augmented Lagrangian function and provide an estimate of the optimal active set. In the second phase, an equality-constrained QP defined by the current active set is approximately minimized in order to generate a second-order search direction. A filter determines the required accuracy of the subproblem solutions and provides an acceptance criterion for the search directions. The resulting algorithm is globally and finitely convergent. The algorithm is suitable for large-scale problems with many degrees of freedom, and provides an alternative to interior-point methods when iterative methods must be used to solve the underlying linear systems. Numerical experiments on a subset of the CUTEr QP test problems demonstrate the effectiveness of the approach. [ABSTRACT FROM AUTHOR]

We describe the implementation of a parallel algorithm that solves a computational geometry problem arising in meshfree methods. We solve the following problem: Given a collection of N d-rectangles Siand Q points xkin ℝd, compute Ik = {i | xk ε Si}. This task is necessary for computing values of the meshfree basis functions Ψi(xk) and can be used to efficiently compute the mass and stiffness matrices for the meshfree Galerkin method. [ABSTRACT FROM AUTHOR]

We describe a reordering particularly suited to the PageRank problem, which reduces the computation of the PageRank vector to that of solving a much smaller system and then using forward substitution to get the full solution vector. We compare the theoretical rates of convergence of the original PageRank algorithm to that of the new reordered PageRank algorithm, showing that the new algorithm can do no worse than the original algorithm. We present results of an experimental comparison on five datasets, which demonstrate that the reordered PageRank algorithm can provide a speedup of as much as a factor of 6. We also note potential additional benefits that result from the proposed reordering. [ABSTRACT FROM AUTHOR]

We propose a new strategy for choosing the optimal parameters for the Gegenbauer reconstruction method based on Chebyshev spectral coefficients under different assumptions on the smoothness of the function f. These parameters are optimal in the sense that the bounds on the truncation and regularization errors were forced to be equal. This strategy is independent on the number of terms N in the Chebyshev expansion of the function f and guarantees exponential convergence as N → ∞ of the Gegenbauer series to f on the intervals of smoothness. The effectiveness of this strategy and exponential convergence are confirmed by numerical examples for functions with varying degrees of smoothness. [ABSTRACT FROM AUTHOR]

The multilevel fast multipole algorithm is based on the multipole expansion, which has numerical error sources such as truncation of the addition theorem, numerical integration, and interpolation/anterpolation. Of these, we focus on the truncation error and discuss its control precisely. The conventional selection rule fails when the buffer size is small compared to the desired numerical accuracy. We propose a new approach and show that the truncation error can be controlled and predicted regardless of the number of buffer sizes. [ABSTRACT FROM AUTHOR]

We consider a large-scale convex minimization problem with nonnegativity constraints that arises in astronomical imaging. We develop a cost functional which incorporates the statistics of the noise in the image data and Tikhonov regularization to induce stability. We introduce an efficient hybrid gradient projection-reduced Newton (active set) method. By “reduced Newton” we mean taking Newton steps only in the inactive variables. Due to the large size of our problem, we compute approximate reduced Newton steps using conjugate gradient (CG) iteration. We also introduce a highly effective sparse preconditioner that dramatically speeds up CG convergence. A numerical comparison between our method and other standard large-scale constrained minimization algorithms is presented. [ABSTRACT FROM AUTHOR]

In a recent work Ng, Chan, and Tang introduced reflecting (Neumann) boundary conditions (BCs) for blurring models and proved that the resulting choice leads to fast algorithms for both deblurring and detecting the regularization parameters in the presence of noise. The key point is that Neumann BC matrices can be simultaneously diagonalized by the (fast) cosine transform DCT III. Here we propose antireflective BCs that can be related to τ structures, i.e., to the algebra of the matrices that can be simultaneously diagonalized by the (fast) sine transform DST I. We show that, in the generic case, this is a more natural modeling whose features are (a) a reduced analytical error since the zero (Dirichlet) BCs lead to discontinuity at the boundaries, the reflecting (Neumann) BCs lead to C0 continuity at the boundaries, while our proposal leads to C¹ continuity at the boundaries; (b) fast numerical algorithms in real arithmetic for both deblurring and estimating regularization parameters. Finally, simple yet significant 1D and 2D numerical evidence is presented and discussed. [ABSTRACT FROM AUTHOR]

We present and analyze several simple algorithms for accurately computing the sum of n floating point numbers using a wider accumulator. Let f and F be the number of significant bits in the summands and the accumulator, respectively. Then assuming gradual underflow, no overflow, and round-to-nearest arithmetic, up to approximately 2F-f numbers can be added accurately by simply summing the terms in decreasing order of exponents, yielding a sum correct to within about 1.5 units in the last place (ulps). We apply this result to the floating point formats in the IEEE floating point standard. For example, a dot product of single precision vectors of length at most 33 computed using double precision and sorting is guaranteed correct to nearly 1.5 ulps. If double-extended precision is used, the vector length can be as large as 65,537. We also investigate how the cost of sorting can be reduced or eliminated while retaining accuracy. [ABSTRACT FROM AUTHOR]

A propagating interface can develop corners and discontinuities as it advances. Level set algorithms have been extensively applied for the problems in which the solution has advancing fronts. One of the most popular level set algorithms is the so-called fast marching method (FMM), which requires total O(N log2 N) operations, where N is the number of grid points. The article is concerned with the development of an O(N) level set algorithm called the group marching method (GMM). The new method is based on the narrow band approach as in the FMM. However, it is incorporating a correction-by-iteration strategy to advance a group of grid points at a time, rather than sorting the solution in the narrow band to march forward a single grid point. After selecting a group of grid points appropriately, the GMM advances the group in two iterations for the cost of slightly larger than one iteration. Numerical results are presented to show the efficiency of the method, applied to the eikonal equation in two and three dimensions. [ABSTRACT FROM AUTHOR]

We present one-sweep parallel algorithms for the inversion of general and symmetric positive definite matrices. The algorithms feature simple programming and performance optimization while maintaining the same arithmetic cost and numerical properties of conventional inversion algorithms. Our experiments on a Cray T3E-600 and a Beowulf cluster demonstrate high performance of implementations for distributed memory parallel computers. [ABSTRACT FROM AUTHOR]

Trigonometric polynomials are widely used for the approximation of a smooth function from a set of nonuniformly spaced samples. If the samples are perturbed by noise, a good choice for the polynomial degree of the trigonometric approximation becomes an essential issue to avoid overfitting and underfitting of the data. Standard methods for trigonometric least squares approximation assume that the degree for the approximating polynomial is known a priori, which is usually not the case in practice. We derive a multilevel algorithm that recursively adapts to the least squares solution of suitable degree. We analyze under which conditions this multilevel approach yields the optimal solution. The proposed algorithm computes the solution in at most O(rM + M2) operations (M being the polynomial degree of the approximation and r being the number of samples) by solving a family of nested Toeplitz systems. It is shown how the presented method can be extended to multivariate trigonometric approximation. We demonstrate the performance of the algorithm by applying it in echocardiography to the recovery of the boundary of the left ventricle of the heart. Key words: trigonometric approximation; Toeplitz matrix; Levinson algorithm; multilevel method [ABSTRACT FROM AUTHOR]

We study numerical approximations of bifurcating solution curves of the von Kŕrmαn equations with simply supported and clamped boundary conditions, respectively. Of special interest here is the splitting of a double bifurcation point into two simple bifurcation points, and the tracing of the associated secondary solution branches, which corresponds to the phenomenon of mode jumping in the buckling of a rectangular plate. A continuation-unsymmetric Lanczos algorithm is utilized for curve-tracking. Our numerical results verify the theoretical prediction of Schaeffer and Golubitsky, namely, mode jumping occurs when the boundary is partially clamped, but not if it is merely simply supported. [ABSTRACT FROM AUTHOR]

A new class of efficient stochastic algorithms for the numerical treatment of co-agulation processes is proposed. The algorithms are based on the introduction of fictitious jumps combined with an acceptance-rejection technique for distributions depending on particle size. The increased efficiency is demonstrated by numerical experiments. In particular, gelation phenomena are studied. [ABSTRACT FROM AUTHOR]

We give an overview of fast algorithms for solving least squares problems with Toeplitz structure, based on generalization of the classical Schur algorithm, and discuss their stability properties. In order to obtain more accurate triangular factors of a Toeplitz matrix as well as accurate solutions for the least squares problems, methods based on corrected seminormal equations (CSNE) can be used. We show that the applicability of the generalized Schur algorithm is considerably enhanced when the algorithm is used in conjunction with CSNE. Several numerical tests are reported, where different variants of the generalized Schur algorithm and CSNE are compared for their accuracy and speed. [ABSTRACT FROM AUTHOR]

The application of nonparametric probability density function estimation for the purpose of data analysis is well established. More recently, such methods have been applied to fluid flow calculations since the density of the fluid plays a crucial role in determining the flow. Furthermore, when the calculations involve directional or axial data, the domain of interest falls on the surface of the sphere. Accurate and fast estimation of probability density functions is crucial for these calculations since the density estimation is performed at each iteration during the computation. In particular the values fn(X1),fn(X2),…,fn(Xn) of the density estimate at the sampled points Xi are needed to evolve the system. Usual nonparametric estimators make use of kernel functions to construct fn. We propose a special sequence of weight functions for nonparametric density estimation that is especially suitable for such applications. The resulting method has a computational advantage over kernel methods in certain situations and also parallelizes easily. Conditions for convergence turn out to be similar to those required for kernel-based methods. We also discuss experiments on different distributions and compare the computational efficiency of our method with kernel based estimators. [ABSTRACT FROM AUTHOR]

The problem of minimizing a sum of Euclidean norms dates from the 17th century and may be the earliest example of duality in the mathematical programming literature. This nonsmooth optimization problem arises in many different kinds of modern scientific applications. We derive a primal-dual interior-point algorithm for the problem, by applying Newton's method directly to a system of nonlinear equations characterizing primal and dual feasibility and a perturbed complementarity condition. The main work at each step consists of solving a system of linear equations (the Schur complement equations). This Schur complement matrix is not symmetric, unlike in linear programming. We incorporate a Mehrotra-type predictor-corrector scheme and present some experimental results comparing several variations of the algorithm, including, as one option, explicit symmetrization of the Schur complement with a skew corrector term. We also present results obtained from a code implemented to solve large sparse problems, using a symmetrized Schur complement. This has been applied to problems arising in plastic collapse analysis, with hundreds of thousands of variables and millions of nonzeros in the constraint matrix. The algorithm typically finds accurate solutions in less than 50 iterations and determines physically meaningful solutions previously unobtainable. [ABSTRACT FROM AUTHOR]

Variable time-stepping algorithms for initial value ordinary differential equations are traditionally designed to solve a problem for a fixed initial condition and over a finite time. It can be shown that these algorithms may perform poorly for long time computations with initial conditions that lie in a small neighborhood of a fixed point. In this regime there are orbits that are bounded in space but unbounded in time, and the classical error-per-step or error-per-unit-step philosophy may be improved upon. A new error criterion is introduced that essentially bounds the truncation error at each step by a fraction of the solution arc length over the corresponding time interval. This new control can be incorporated within a standard algorithm as an additional constraint at negligible additional computational cost. It is shown that this new criterion has a positive effect on the linear stability properties and hence improves behavior in the neighborhood of stable fixed points. Furthermore, spurious fixed points and period two solutions are prevented. The new criterion is shown to be admissible in the sense that it can always be satisfied with nonzero stepsizes. Implementation details and numerical results are given. [ABSTRACT FROM AUTHOR]