SCIENTIFIC PROGRAMS AND ACTIVITIES

March 28, 2024

Seventh IMACS International Symposium on Iterative Methods in Scientific Computing

May 5-8, 2005
The Fields Institute for Research in the Mathematical Sciences
and the University of Toronto, Ontario, Canada

Abstracts of talks

[Invited talks] [Contributed talks] [Minisymposia talks] [Student papers' talks]


Invited talks



Duality-Based Iterative Methods for Total Variation Minimization


Tony Chan
Mathematics Department
UCLA
Los Angeles, CA 90095-1438, U.S.A.

Regularization based on total variation, and related higher order variants based on curvature of level sets, have been used very successfully in image processing. However, their non-smooth nature poses major challenge for numerical methods designed for their iterative numerical minimization. Primal based methods have to be regularized themselves, and even then are often only linearly convergent. Consequently, there has been interest in developing duality-based iterative methods for their solution. The advantage of these dual methods is that the objective function is smooth (often quadratic), but accompanied by many constraints, which may or may not be active at the solution. In this talk, I'll review some of these developments, starting with the PhD thesis of J. Carter, then reviewing the primal-dual method of Chan-Golub-Mulet (related to work by Conn and Overton), to a recent major advance by A. Chambolle for a all-dual formulation. I'll finish by presenting some recent work of mine (in collaboration with J.F. Aujol, S. Esedoglu and F. Park) on extending Chambolle's approach to higher order variants of TV and their applications to computing decomposition of images into geometric, texture and noise components.



Minimizing VaR, CVaR and Hedging Issues for a Portfolio of Derivatives.


Thomas F Coleman
Computer Science and Applied Mathematics
Director, Cornell Theory Center
New York, NY 10004-2501, U.S.A.

Increasingly, financial institutions are concerned with portfolios of derivative instruments. Basic operations include designing effective hedging strategies and constructing investment rebalancing schemes to minimize value-at-risk (VaR) or conditional VaR. Unfortunately, for a portfolio of financial derivatives these problems are typically either ill-posed or ill-conditioned. We illustrate why this is the case, the negative computing implications of this ill-conditioning, and propose problem modifications to yield better-conditioned problems. In addition, we illustrate how some very expensive CVaR computations can, with a problem reformulation, be solved by efficient numerical procedures.

This work is joint work with colleagues Yuying Li, Katharyn Boyle, and Siddharth Alexander.



Derivative Free Optimization -- Some New Results


Andrew R. Conn
IBM

Derivative free optimization methods have been extensively developed in the past decade. To be able to employ the well developed convergence theory of derivative based methods, the models used (whether in a trust region or a line search approach) have to satisfy Taylor-like error bounds. We will present a unified framework to study these error bounds, initially for the case of polynomial interpolation. These bounds depend on the geometry of the interpolation set. We will analyze this geometry and present generalisations which include a viable way of measuring the quality of the geometry. We will also present extensions to least squares and minimum norm methods. Practical techniques of ensuring that the geometric requirement is satisfied and of improving the geometry of a sample set, as well as connections with other results, will also be included.

Joint work with Katya Scheinberg and Luis Vicente.



Spectral Element Multigrid for the Incompressible Navier Stokes Equations


Paul F. Fischer
Mathematics and Computer Science Division
Argonne National Laboratory
Argonne, IL 60439, U.S.A.

We consider application of recently developed spectral element multigrid (SEMG) methods to the study of unsteady incompressible flows. Our simulation approach is based on a semi-implicit formulation, in which the space-discretized Navier-Stokes equations are split into independent convective, viscous, and pressure substeps. Our multigrid solver is developed for the consistent Poisson problem governing the pressure arising from the lP_N - lP_{N-2} spectral element discretization of the Stokes problem. The pressure is associated with the fastest timescale and is therefore the leading contributor to stiffness and computational complexity in the Navier-Stokes time advancement.

The spectral element method (SEM) is a high-order weighted residual discretization similar to p-type finite element methods (FEMs). The SEM, however, is generally restricted to hexahedral elements (quadrilateral in lR^2) with nodal bases located on tensor products of Gauss-Legendre (GL) or Gauss-Lobatto-Legendre (GLL) quadrature points. The GL and GLL bases provide stability and allow efficient evaluation of quadratures in the weighted residual formulation, while the tensor product forms allow for efficient local evaluation of the operators. For a discretization involving E elements of order N (i.e., ~ E N^3 gridpoints in lR^3), the SEM operator evaluation and storage complexities are respectively only O(E N^4) and O(E N^3), which is in marked contrast to the O(E N^6) costs incurred by the standard p-type FEM. Typical values of N for the SEM are in the range N = 4 ­to 15.

Following the early work of Rønquist and Patera, and of Maday and Munoz, our SEMG method is based on intra-element restriction and prolongation, with coarse-grid spaces derived from on a reduction of N within a fixed mesh (i.e., fixed E). The earlier SEMG approaches were based on pointwise-Jacobi smoothing. This was shown numerically and theoretically to yield order-independent convergence rates in lR^1 with a multigrid convergence factor of rho = .75 for the Poisson equation. In lR^2, however, the performance deteriorates to rho = O(N^{1/2}) because the tensor-product GL and GLL bases yield high-aspect ratio cells within each element that are not amenable to pointwise smoothing.

Here, we develop an approach that extends to two and three dimensions. We avoid the cell aspect-ratio difficulty by using Schwarz-based smoothing in which one solves directly, using fast tensor-product solvers, local Poisson problems over subdomains that are taken to be extensions of each element. The coarsest level problem is derived from a linear finite element discretization based on a triangulation of the quadrilateral (or hexahedral, in 3D) spectral element mesh. We demonstrate that effective smoothing requires the local Schwarz solutions to be weighted by the inverse of the diagonal counting matrix that represents the number of subdomains sharing a given vertex. We compare this hybrid Schwarz-multigrid approach with other multi-level solvers for the Poisson equation, the Helmholtz equation, and for the Schur complement (consistent Poisson) matrix that governs the pressure in our unsteady Navier-Stokes computations. We demonstrate for several incompressible flow simulations that SEMG provides a two-fold reduction over our earlier two-level additive-Schwarz based approach.

Joint work with James W. Lottes.



Recent Advances in Algebraic Multigrid


Van Emden Henson
Center for Applied Scientific Computing
Lawrence Livermore National Laboratory

Over the past decade, Algebraic Multigrid (AMG) has become a prominent tool for the rapid solution of unstructured-grid problems. For nearly 20 years after its inception, AMG was largely ignored, a victim of the overhead and setup costs necessary to create the components of the algorithm. So long as geometrically structured grids were the norm, this condition persisted. In recent years, however, as unstructured grids became more common and as the problem sizes grew to demand large-scale massively parallel computers, AMG came into its own. In this talk we describe AMG at a broad-brush level, indicating its relation to geometric multigrid. We discuss some pertinent facts of the history of AMG, mention the underlying philosophy of the method, and describe briefly the major families of AMG algorithms and their underpinnings. We then concentrate on the more recent advances in AMG, and in particular how an intense period of research primarily into perfecting the prolongation operators has led to a contemporary focus on coarse-grid selection, compatible relaxation, and adaptive AMG algorithms. We describe some of the recent advances in AMG theory. We conclude the presentation with some general observations and remarks designed to set the stage for the AMG presentations in minisymposia (M2 and M3) that follow.



Analysis and Computation of Google's PageRank


Ilse Ipsen
Department of Mathematics
North Carolina State University
Raleigh, NC, U.S.A.

The PageRank of a web page is a major ingredient in how the search engine Google determines the ranking of web pages. Computing the PageRank amounts to computing the stationary distribution of a stochastic matrix whose size is now in the billions.

We discuss the effect on the PageRank when links are added and removed, and when the personalization vector or damping factor are changed. We also discuss different approaches for computing PageRank.



Continuation Algorithms for Parameter Dependent
Compact Fixed Point Problems


Tim Kelley
Department of Mathematics, Box 8205
North Carolina State University
Raleigh, NC 27695-8205 USA

In this talk we show how compactness can be used to design and analyze Newton-Krylov and multilevel algorithms for parameter-dependent compact fixed point problems. Such problems arise not only as integral equations, but also from indirect integrations such as source iteration in radiative transfer calculations and time-stepping approaches to dynamic analysis of parabolic partial differential equations.

The methods we propose are easy to implement and require no detailed knowledge of the internal structure of the fixed point map. We will illustrate the ideas with an example from computational chemistry.



A Novel Multigrid Based Preconditioner For
Heterogeneous Helmholtz Problems


Cornelis W. Oosterlee
Delft University of Technology,
Delft Institute of Applied Mathematics (DIAM),
the Netherlands.

An iterative solution method, in the form of a preconditioner for a Krylov subspace method, is presented for the Helmholtz equation. The preconditioner is based on a Helmholtz type differential operator with a complex term. A multigrid iteration is used for approximately inverting the preconditioner. The choice of multigrid components for the corresponding preconditioning matrix with a complex diagonal is made with the help of Fourier analysis. Multigrid analysis results are verified by numerical experiments. High wavenumber Helmholtz problems in heterogeneous media are solved indicating the performance of the preconditioner.



Fast Solvers for Incompressible Flow


Andy Wathen
Oxford Universit
United Kingdom

The Navier-Stokes equations describing the flow of an incompressible viscous fluid lead via discretization and lineaization to non-symmetric matrices with eigenvalues in both halves of the complex plane. Related subproblems are the Stokes equations which are symmetric but indefinite and the convection-diffusion equations which are non-symmetric but have eigenvalues with only positive real parts. The discrete Poisson operator represented by a symmetric and positive definite matrix is a component of all of these problems.

We will use this sequence of problems to describe how fast preconditioned iterative solvers can be built for the more and more challenging problems described by these various models of incompressible fluid flow. In particular we will show how good preconditioners respect the natural norms associated with the originating PDE problems and thus lead to iterative convergence in the norms relevant for the underlying problem.

This work is joint with Howard Elman and David Silvester.


Contributed talks



High order evaluation of European and American options
in Levy markets of finite activity


Ariel Almendral (email: a.almendral@ewi.tudelft.nl)
Delft Technical University, Delft, The Netherlands
Co-authors: Cornelis W. Oosterlee

Vanilla option prices satisfy an Integro-Differential equation, when the underlying asset follows a Levy process. We propose two different approaches to compute the numerical solution in the finite variation case.In the first approach, a graded mesh is used in combination with Simpson's rule. In the second approach, a connection with Volterra integral equations with weakly singular kernels is established, opening a possibilty for an application of known solution methods from the theory of Volterra Integro-Differential equations. Both approaches produce a second order convergence rate.



Energy-Minimizing Pulse Design for Magnetic Resonance Imaging
Using Interior Point Methods, Remez Exchange, Surrogate Models,
and Symbolic Computation in Lie Groups


Christopher Kumar Anand (email: anandc@mcmaster.ca)
McMaster University
Co-authors: Andrew Thomas Curtis

Unlike x-ray and CT imaging, the radio-frequency radiation used in Magnetic Resonance Imaging is low-energy and not thought to pose any long-term health risk. However, it does carry enough energy to heat the patient. To prevent patient discomfort the energy is limited by government regulation. With advances in other aspects of MR Imagers, these limits constrain the types and rapidity of clinical imaging, including the majority of body imaging in high-end 3 Tesla imagers. A new class of energy-minimizing (rf,gradient) pulse pairs is needed to remove this limitation. This has previously been formulated as an optimal control problem, but the simultaneous optimization of both controls has not been practical with current tools. So, we are developing a dedicated iterative optimizer.

In this talk, we will explain the structure of the control problem as a continuously-parametrized family of ODEs with a common set of piecewise-constant controls. We have developed an exact method of integration for such systems using symbolic calculation in Lie Groups, and a computational technique for obtaining gradients in O(N) time, with measured computation time much less than double the time for simple function evaluation. The overall optimization framework (1) approximates the continuous family of terminal constraints using a Remez-like algorithm; (2) approximates the penalized constraints with a low-dimensional quadratic surrogate; (3) defines an unconstrained problem incorporating both penalties and barriers for the different types of path constraints; and (4) solves the resulting problem using an interleaved Interior-Point/Remez algorithm.



A domain decomposition method to parallelize a constrained 3D Maxwell solver


Franck Assous (email: franckassous@netscape.net)
Research Institute, College of Judea and Samaria, Ariel & Dpt of Math.
and Stat., Bar-Ilan Univer. Ramat-Gan - Israel
Co-authors: Jacques Segre, Eric Sonnendrucker

In order to be able to simulate complex devices, we need in some cases a full three-dimensional code for the resolution of the Vlasov-Maxwell equations. A three-dimensional code has been written for this purpose and has already been used for many applications. However this code is sequential and on several occasions, it appears worthwhile to use multiple processors. This makes it indispensable to develop a version of the code that will run fast and reliably on parallel machines. In order to develop a parallel version that scales well when using many processors we choose to use a domain decomposition approach moreover this option allows us to reuse a large part of the sequential code for the resolution on each subdomain.

Let us now recall the different kind of Maxwell solvers and their features with respect to parallelization and point out the specificities of ours. The most used and fastest Maxwell solvers are based on Yee's scheme on uniform grids. These are completely explicit and straightforward to parallelize, at least when the continuity equation ¶t r+ div J=0 is numerically verified. However in complicated geometries or for some other reason it might be preferable to use unstructured grids. Several kind of methods have been developed in this case : Delaunay-Voronoi Finite Volume methods have the same nice properties concerning parallelization as those on structured grids. Their drawback is that they need Delaunay-Voronoi meshes which are difficult to generate in 3 dimensions. Other types of methods include traditional vertex centered or cell centered finite volume methods which are explicit and thus can be easily parallelized in the same conditions as finite difference methods. And finite element methods based either on edge elements or as in our case on Taylor-Hood elements. Concerning edge element solvers, the mass matrix needs to be inverted, so the parallelization problems are alike those arising for elliptic problems. In our case, i.e. Taylor-Hood elements, the mass matrix is diagonal. But the formulation is based on a constrained formulation of the Maxwell equations, which makes it impossible to have a completely explicit formula for the time advance and thus makes the parallelization challenging.

In our talk, we first shall recall the constrained wave equations formulation of Maxwell's equations which is used in the code. Then we shall introduce a variational formulation adapted to the given domain decomposition, the continuity at the interfaces being imposed by duality using Lagrange multipliers. After that we shall describe the discretization of this variational formulation and derive a linear system suitable for multiprocessor resolution. The preconditioned Uzawa algorithm used for the resolution of this system shall be described. And finally we shall present numerical results demonstrating the efficiency of the method.



Peano Curves and Cache Oblivious Multiplication of Full and Sparse Matrices


Michael Bader (email: bader@in.tum.de)
Institut für Informatik, TU München, Boltzmannstr. 3, 85748 Garching, Germany
Co-authors: Christoph Zenger

Space filling curves have become a valuable tool in many applications, where data locality is an issue. A well known example in scientific computing is the parallelization of data during the discretization of partial differential equations. In this context, space filling curves have quasi-optimal locality properties. We would like to present several approaches that try to benefit from these locality properties in the context of cache efficient algorithms.

For the multiplication of two matrices we present an approach, where the matrix elements are mapped to the linear memory space via an index function that resembles iterations of the Peano space filling curve. The respective multiplication algorithm can then be designed such that temporal and spatial locality in the data access is optimized. With respect to temporal locality, the algorithm guarantees that during any k3 operations only O(k2) elements are accessed. Every element will be reused approximately k times. For spatial locality, the algorithm totally avoids jumps in memory space: matrix elements will either be directly reused, or the next accessed element will be one of its direct neighbours in memory. This is not only advantageous for cache reuse, but may also be an important property for hardware implementations.

The presented approach can also be generalized to store sparse matrices. Then, an adaptively refined iteration of the Peano curve is used for the index mapping. Locality properties are retained to a large extent.

We will also discuss a related approach to deal with sparse matrices, and especially with the corresponding systems of equations that result from discretizations of partial differential equations. In that case, optimal locality can be retained, if we allow a collection of stacks that hold intermediate results. This can even be extended to certain kinds of adaptive or multilevel discretizations.



Numerical solution for a type of interface problems


Pallav Kumar Baruah (email: baruahpk@yahoo.com)
DMACS, Sri Sathya Sai Institute of HIgher Learning, INDIA

We come across the type of interface problems where two different ordinary differential equations are defined on two different intervals with a common interface point and the solution satisfies certain matching condition at the interface. The analytical solution for this kind of problem have been studied. But there is no systematic study to find numerical solution of this type of problems. Here we like to present analysis and results to show that exising methods for IVPs and BVPs can be extended to find solution to these problems.



Computational Method for Solving Two Point Boundary Value Problems
Using Parametric Cubic Spline


Rajesh K. Bawa (email: rajesh_k_bawa@yahoo.com)
Department of computer science, Punjabi University, Patiala, INDIA

In this paper, a computational method using parametric cubic spline for general linear second order two point boundary value problems is analysed. The method is shown to have second and fourth order convergence for specific choice of parameters present. Some test examples are given to support the theoretical estimates.



Generalized Thermo-elastic Problem of a Plate


Sukhwinder Kaur Bhullar (email: sbhullar_2000@hotmail.com)
Centre for Advanced Studies in Mathematics, Panjab University, Chandigarh, India

When a thermo-elastic body is suddenly heated or cooled from out side, heat flow arises resulting in change in the temperature distribution in the body and development of stress field. In many engineering problems, thermoelastic effects due to a temperature distribution in the structure have to be considered. This paper aims to study the developments of these fields in a homogeneous, isotropic elastic plate of thickness ?H?, initially has temperature T0 and is in stress free state .The variation in temperature field and stress field occurs owing to the action of external loadings. Temperature and stress distribution have been determined using appropriate initial and boundary conditions. The numerical computations are carried out for a plate of stainless steel to display the variation of temperature and stress fields. It is concluded that the temperature distribution is symmetrically distributed and variation in stress component (tau)xy is more in left side and less in right side whereas stress component


Peano Curves and Cache Oblivious Multiplication of Full and Sparse Matrices

Michael Bader (email: bader@in.tum.de)
Institut für Informatik, TU München, Boltzmannstr. 3, 85748 Garching, Germany
Coauthors: Christoph Zenger

Space filling curves have become a valuable tool in many applications, where data locality is an issue. A well known example in scientific computing is the parallelization of data during the discretization of partial differential equations. In this context, space filling curves have quasi-optimal locality properties. We would like to present several approaches that try to benefit from these locality properties in the context of cache efficient algorithms.

For the multiplication of two matrices we present an approach, where the matrix elements are mapped to the linear memory space via an index function that resembles iterations of the Peano space filling curve. The respective multiplication algorithm can then be designed such that temporal and spatial locality in the data access is optimized. With respect to temporal locality, the algorithm guarantees that during any k3 operations only O(k2) elements are accessed. Every element will be reused approximately k times. For spatial locality, the algorithm totally avoids jumps in memory space: matrix elements will either be directly reused, or the next accessed element will be one of its direct neighbours in memory. This is not only advantageous for cache reuse, but may also be an important property for hardware implementations.

The presented approach can also be generalized to store sparse matrices. Then, an adaptively refined iteration of the Peano curve is used for the index mapping. Locality properties are retained to a large extent.

We will also discuss a related approach to deal with sparse matrices, and especially with the corresponding systems of equations that result from discretizations of partial differential equations. In that case, optimal locality can be retained, if we allow a collection of stacks that hold intermediate results. This can even be extended to certain kinds of adaptive or multilevel discretizations.



A-Posteriori Finite Element Bound Method devised by
an Adaptive Refinement, the Direct Equilibration and
a Parallel Computing Strategies for the
Multi-physical, Multi-scale and Multi-Dimensional
Partial Differential Equations


Hae-Won Choi (email: haewon@mie.utoronto.ca)
Department of Mechanical Engineering, University of Toronto
Co-authors: Marius Paraschivoiu (Concordia University)

As numerical techniques widely are applied to multi-physical, multi-scale and multi-dimensional engineering applications, the ability to compute solutions and furthermore the ability to quantify the accuracy of underlying mathematical models seem to be significant importance for robust and reliable numerical simulations. Over last two decades, a number of endeavors have focused on the challenge for development and improvement of notable error estimation and adaptive refinement methods. Significant progress has been dedicated in recent years by the use of the advanced a-posteriori finite element error estimation technique termed the bound method. This novel technique provides fast yet robust, reliable and efficient bounds to the `output' of underlying partial differential equations (PDEs) as well as naturally yields a local error indicator, which estimates numerical error for given mathematical models, leading to construct economical and optimized meshes for computations. The bound method in this work is applied to the multi-physical, micro-scale and multi-dimensional PDEs, in particular to facilitate engineering applications governed by the incompressible Navier-Stokes and Energy equations. Furthermore, the bound method devised by advanced numerical strategies, such as an adaptive refinement, the direct equilibration and a parallel computing techniques, will address inter-disciplinary and emerging applications: heat and mass transport phenomena in an array of electronic chip devices and heat exchangers; as well as fluidic motion flow and species transport phenomena of the electro-osmotic flow in microchannels for integrated microfluidic systems such as biochip and Lap-on-chip devices.



Advanced Numerical-Analytical Methods for Path Integral Calculation
and Its Application to Some Famous Problems of 3-D Turbulence Theory


Jaykov Foukzon (email: advancedguidance@list.ru)
Israel, Tel-Aviv, st.Rambam 7a2

Numerical-analytical study of the three-dimensional nonlinear stochastic partial differential equations, analogous to that proposed by V. N. Nikolaevskii [Recent Advances in Engineering Science (Springer - Verlag, Berlin. 1989)] to describe longitudinal seismic waves, is presented. The equation has a threshold of short-wave instability and symmetry, providing long-wave dynamics. The results of computation are in a sharp contradiction with Tribelsky M.I and K. Tsuboi's work (Phys. Rev. Lett. 76 1631 (1996)), in which the influence of the thermal fluctuations was not taken into account and a principally erroneous scheme of numerical modeling of chaos on the lattice was used. The proposed new mechanism for quantum chaos generates nonlinear dynamical systems. The hypothesis is that physical turbulence could be identified with quantum chaos of considered type.



A preconditioned fast method for higher-order accurate quantification
of perturbation effects in nuclear systems


R. van Geemert (email: rene.vangeemert@framatome-anp.com)
Framatome ANP GmbH, Freyeslebenstrasse 1, 91058 Erlangen (Germany)

Adequate quantification of the effects of nuclear system parameter variations, as well as of material composition perturbations, on the lowest mode solution (the three-dimensional flux shape) of the neutronics eigenvalue equation is of importance in various areas in nuclear reactor technology. Examples are nuclear uncertainty analysis, sensitivity analysis for reactor transients and design optimization (i.e. fuel lattice design optimization, control rod movement planning, nuclear shielding design or in-core fuel management optimization).

A numerical development, utilizing physically meaningful, iteratively obtained preconditioners, is presented which is aimed essentially at enabling a fast higher order accurate treatment of the flux shape change caused by localised material composition changes (such as material choice variations, burnable poison density variations at lattice level, or complete fuel assembly permutations at core level). This fast method is particularly useful when the effects of large numbers of variations are to be assessed.

According to previously pursued studies in the numerical reactor physics community, acquisition of higher-order accurate results for the perturbed flux shape either entailed the use of preconditioned iterative methods, or required the availability of an extensive set of higher reference eigensolutions of the system equations. Both of these methodologies excluded the option to obtain a higher-order accurate estimate of a perturbed flux defined in only one or a limited number of localised spatial subregions in the system, without the need to compute the perturbed fluxes in all other positions. For example, the exact change in power density in one specific fuel element in a reactor, due to a material perturbation imposed elsewhere in the system (for example, due to control rod movements or burnable poison concentration variations) could follow only from an iterative procedure defined for the entire flux distribution.

The methodology presented here enables the setup of what could be interpreted as a polynomial form. This form features the property that the effect on a specific, only locally defined flux, caused by system variations elsewhere that are localised in space as well, can be expanded polynomially up to higher order accuracy, with the specific imposed variations themselves as functional arguments. For specific cases of interest, the application of this method leads to substantial savings in computational effort, while obtaining similarly accurate results. Numerical investigations, showing the accuracy and computational efficiency of the method, are reported, and possible application areas are discussed.



A general framework for recursions for Krylov space solvers


Martin H. Gutknecht (email: mhg@math.ethz.ch)
ETH Zurich

Krylov space methods for solving linear systems of equations come in many flavors and use various types of recursions to generate iterates xn (approximate solutions of Ax = b), corresponding residuals rn : = b - A xn, and direction vectors (or search directions) vn : = (xn+1 - xn) / wn. Starting from a general definition for a Krylov space solver we give necessary and sufficient conditions for the existence of the the various types of recursions, and we recall the relations that exist between the matrix representations of these recursions. Much of this is more or less well known, but there are also some new, perhaps even surprising aspects.



Interpolation of numerical solutions of PDEs at off mesh points
using iterative collocation


Samir Hamdi (email: samir.hamdi@utoronto.ca)
1145 Hunt Club Road, Suite 500, Ottawa, Canada, K1V OY3
Co-authors: W. H. Enright J. J. Gottlieb, and W. E. Schiesser

In this study, we first present a review of standard interpolating algorithms (based on splines and Hermite approximations) and then we introduce a new interpolation procedure proposed recently by Enright (2000) as an alternative to traditional interpolation methods. This new interpolation, called Differential Equation Interpolation (DEI), is based on associating a set of collocation points with each mesh element and requiring that the local approximation interpolates the meshpoint data and almost satisfies the PDE at the collocation points. We show, both theoretically and through numerical experiments, that the accuracy of this interpolation is higher than the accuracy of traditional interpolation and generally consistent with the meshpoint accuracy of the PDE solution. For most linear PDEs, the DEI is obtained by exact collocation and the associated linear system is solved analytically using Maple symbolic computation. For nonlinear PDEs, the DEI is based on almost collocation and the resulting nonlinear system is solved using an iterative Newton method. Examples are given for off mesh interpolation of numerical solutions of evolution PDEs such as the Korteweg-de Vries equation and regularized long wave equation and also for the integration of several invariants of motion related to the conservation of mass, momentum and energy.



Dispatching Buses in a Depot Minimizing Mismatches


Mohamed Hamdouni (email: mohamed.hamdouni@polymtl.ca)
Département de Mathématiques et génie Industriel, École polytechnique & GERAD
Co-authors: François Soumis, Guy Desaulniers

We consider the problem of assigning parking slots to buses of different types as soon as their arrivals such that each bus can leave the depot next morning without maneuvers (i.e, rearrangement of buses within lanes), this problem is recognized in literature by "Dispatching Problem". We introduce a class of robust solutions for this problem where each depot lane is partitioned at most two blocks, lanes with one-block (containing one bus type) and lanes with two-block (containing two bus types gathered together, one block for each type). Each block have to contains only the buses of the right type, and the problem may not have a solution. So, the assignment an arrival bus to block of wrong type yield incompatible arrival, likewise, the assignment a bus that parked in block of a given type to departure that required a different bus type(each departure require a bus type) yield incompatible departure. Defining a mismatches by the assignment of a bus to a departure of wrong type. First, we formulate a model in witch the depot lanes arefilled according to specific patterns one-block patterns and two-block patterns, that minimize the number of incompatible arrivals and incompatibledepartures. A computational results are presented. Second, we extend the previous model to a model that minimize the number of mismatches. We propose Bender's decomposition as resolution approche, and a computational study is performed.



Iteration Method for Integro-Differential Equations


K. Ivaz (email: Ivaz@tabrizu.ac.ir)
Department of Mathematics,
Shabestar Islamic Azad University, Shabestar, Tabriz, Iran

In this paper we give an iteration method for solving non-linear integro-differential equations.



Overlapping level for the restricted version of the
overlapping Schur complement preconditioner


Zhongze Li (email: lzz@lsec.cc.ac.cn)
Institute of Computational Mathematics, Chinese Academy of Sciences, P.R. China

[2] presents a preconditioner based on solving approximate Schur complement systems with overlapping restricted additive Schwarz methods (RAS) [1]. The construction of SchurRAS is as simple as in the standard RAS. The communication requirements for each application of SchurRAS are similar with those of the standard RAS approach. For a tested MHD problem, the number of iterations required by SchurRAS exhibits moderate growth with the increase in the number of processors. In the tests, the residual norm reduction by 10-6 was achieved by flexible FGMRES(60). In preconditioning, ILUT with lfil = 70 and the dropping tolerance 10-3 was used as a choice of an incomplete LU factorization. The tests were performed on 1.7Ghz IBM p690 at MSI. The SchurRAS takes about 234 steps to converge on 16 processors. However, iteration counts required by the RAS is 3431. Because the computational costs of SchurRAS are similar to those of RAS, the reduction in the iteration number also leads to the reduction in the execution time: RAS takes 654.39 seconds to converge on 16 processors, SchurRAS takes only 47.88 seconds.

In [2], the overlapping level is 1. This talk will investigate how the overlapping level affects the performance of RAS and SchurRAS respectively.

References
[1] X.-C. Cai and M. Sarkis. A restricted additive Schwarz preconditioner for general sparse linear systems. SIAM Journal on Scientific Computing, 21(2):792-797, 1999.
[2] Z. Li and Y. Saad. SchurRAS: A restricted version of the overlapping Schur complement preconditioner. Technical Report 2004/77, Minnesota Supercomputing Institute, 2004.



Spectral-Viscosity Method for Dynamic Traffic Flow Simulation


Shih-Ching Lo (email: sclo@nchc.org.tw)
National Center for High-Performance Computing

The behavior of road using is one of the complexities of living in today's world. As the trend of increasing travel demand, planning, design, prediction, control and management of the transportation system become more and more important. Traffic flow theory provides analytical techniques and the description of the fundamental traffic flow characteristics, such as flow, density and speed. This new science has addressed questions related to understand traffic processes and to optimize these processes through proper design and control. To make traffic management strategies (such as ramp metering, congestion tolls, entrance/exit control and so on) become more efficient, dynamic traffic flow model is necessary. There are four main methodologies of modeling dynamic traffic flow: car-following model, kinetic model, Boltzmann-like model and cellular automation. The kinetic model is the simplest of the methodologies and has advantages of having analyzable equations and rapidly computational properties. In 1955, Lighthill, Whitham and Richards firstly proposed a macroscopic kinetic traffic flow model (LWR model), which is a continuity equation. The model is derived from the conservation of vehicle numbers. With different boundary conditions, the LWR model can be considered to describe traffic flow under different control strategies. Periodical boundary conditions are suggested for describing signalized intersections or ramp metering interchanges. Furthermore, the LWR model is employed to predict the formation or dispersion of platoon with different initial conditions. Although the LWR model is simple, it is still a nonlinear model. Also, the formation of shock wave is a problem. In the viewpoint of numerical approximation, to handle the nonlinearity of the model and the discontinuity of the solutions are the difficult parts in simulation. Spectral methods are considered good approximations of nonlinear conservation laws. However, using spectral methods for nonlinear conservation laws will result in the formation of the Gibbs phenomenon once spontaneous shock discontinuities appear in solution. In this work, an enhanced spectral viscosity method is applied to stabilize the spectral method by adding a spectrally small amount of high-frequencies diffusion carried out in the dual space and a post-processing method, which removes the spurious oscillations at the discontinuities. An enhanced shock location detector is developed to ensure the post-processing method successful. In addition, the shock location detector can be used to detect the location of traffic incident (such as location of traffic jam or accidents), where discontinuity of density occurs.



Decoupled and Iterative Method for Numerical solution of
Three-Dimensional Density-Gradient Model in Semiconductor Devices Simulation


Shih-Ching Lo (email: sclo@nchc.org.tw)
National Center for High-Performance Computing
Co-authors: Yiming Li

Numerical solution of a set of semiconductor partial differential equations (PDEs) has been of great interest in these years. Nowadays, the dimension of semiconductor devices, such as MOSFETs is in the regime of nanometer. Quantum mechanical confinement effect becomes an important factor in numerical simulation of semiconductor devices. Among approaches, density-gradient (DG) approach has successfully been proposed to account for the effect of quantum mechanical confinement. A set of DG model involves the Poisson equation, the electron current continuity equation, and the electron quantum corrected density equation to be solved for a specified device domain. With the continuous decrease of device dimensions, one- and two-dimensional simulations are not accurate enough to describe and explore the physical transportation phenomena for electrons and holes in a semiconductor device.

In this work, we iteratively solve the 3D DG model for a given nanoscale semiconductor device, where several physical quantities and convergence properties are discussed. The solution procedure is described briefly as follows. (1) the stop criteria, mesh, output variables and simulated devices are chosen, (2) solving the Poisson equation with density-gradient correction modified potential iteratively, (3) after the Poisson equation is convergent, continuity equations are solved, (4) then, we'll check the whole system converges or not, (5) and if the whole system converges, then stop computing, otherwise, the Poisson equation and continuity equations should be solved again until the whole system equations converge. All decoupled PDEs to be solved are approximated with the finite element method over 3D domain. This scheme makes sure the solution will be self-consistent. By calculating several essential characteristics of devices, such as the threshold voltage (VTH) and the drain current (ID), the device physical quantities are compared with the results of 1D and 2D simulations by varying channel lengths. Also, the convergence behavior in the inner and outer iteration loops is discussed among 1D, 2D, and 3D simulations.



A Lanczos method for solving symmetric indefinite systems


Roummel Marcia (email: marcia@math.wisc.edu)
Department of Biochemistry, University of Wisconsin-Madison

The Lanczos method is an iterative method for solving the linear system $Ax = b$, where $A$ is large, sparse, symmetric, and positive definite. This method breaks down if $A$ is indefinite and the $LDL^T$ factorization of the tridiagonal matrix $T$ arising from the underlying Lanczos process does not exist. In this talk, we present a pivoting strategy for a \textsl{block} $LDL^T$ factorization of $T$. We demonstrate normwise backwards stability and how this factorization can be applied to the Lanczos method.



Simulation of Phase Combinations in SMA Patches with
Hybrid Optimization Methods


Roderick Melnik (email: rmelnik@wlu.ca)
Wilfrid Laurier University
Co-authors: Linxiang Wang,
Mads Clausen Institute, University of Southern Denmark

Microstructures in Shape Memory Alloys (SMA) hold the key to the unique properties of these materials. However, the numerical simulation of the SMA microstructures is quite involved due to non-convexity of the associated free energy functional, coupling between mechanical and thermal fields, and strong nonlinearities.

In this presentation, we analyse numerically phase combinations in a SMA patch by reformulating the original problem as an energy minimization problem. Our main emphasis is given to square-to-rectangular phase transformations, induced in a low-temperature SMA patch by applying either thermal or mechanical loadings.

The free energy density of the SMA patch is constructed on the basis of the modified Ginzburg-Landau theory for a specific case of square-to-rectangular transformations. Due to non-convexity of the Ginzburg-Landau free energy functional, we have to deal with many local minima.

First, we apply the multiple population genetic algorithm to search for an estimation of the global minimizer. Then, by using this estimation as an initial guess, we apply the sequential quadratic programming based on the quasi-Newton iteration method to refine the global minimum with a higher accuracy. By this hybrid optimization method, we are able to model laminar microstructures in a SMA patch with and without external loadings, and several specific examples will be discussed in this talk.



Analysis of Krylov Subspace Recycling for Sequences of Linear Systems


Michael L. Parks (email: mlparks@sandia.gov)
Sandia National Laboratories
Co-authors: Eric de Sturler (University of Illinois at Urbana Champaign)

Many problems from engineering and the sciences require the solution of sequences of linear systems where the matrix and right-hand side change from one system to the next, and the linear systems are not available simultaneously. We review a class of Krylov subspace methods for sequences of linear systems, which can significantly reduce the cost of solving the next system in the sequence by "recycling" subspace information from previous systems. These methods have been successfully applied to sequences of linear systems arising from several different application areas. We analyze a particular method, GCRO-DR, that recycles approximate invariant subspaces, and establish residual bounds that suggest a convergence rate similar to one obtained by removing select eigenvector components from the initial residual. We review implications of this analysis, which suggests problem classes where we expect this technique to be particularly effective. From this analysis and related numerical experiments we also demonstrate that recycling the invariant subspace corresponding to the eigenvalues of smallest absolute magnitude is often not the best choice, especially for nonsymmetric problems, and that GCRO-DR will, in practice, select better subspaces. These results suggest possibilities for improvement in the subspace selection process.



Mathematical Modeling of Large Forest Fire Initiation


Valeriy Perminov (email: pva@belovo.kemsu.ru)
Belovo Branch of Kemerovo State University

A large technogeneous or space catastrophe, as a rule, is known to accompanied by the initiation of large forest fires. In this paper attention is given to questions of description of the initial stage in the development of a mass forest fires initiated by high altitude radiant energy source (for example Tunguska celestial body fall and etc.). Mathematical model of forest fire was based on an analysis of known experimental data and using concept and methods from reactive media mechanics. Within the framework of this model, the forest and combustion products are considered as a homogeneous two temperatures reacting medium. The temperature of condensed (solid) phase includes a dry organic substance, moisture (water in the liquid-drop state), condensed pyrolysis and combustion products and mineral part of forest fuels. In the gaseous phase we separate out only the components necessary to describe reactions of combustion (oxygen, combustible products of pyrolysis of forest fuels and the rest inert components). It is considered that 1) the flow has a developed turbulent nature, molecular transfer being neglected, 2) gaseous phase density doesn't depend on the pressure because of the low velocities of the flow in comparison with the velocity of the sound, 3) forest canopy is supposed to be non-deformed porous medium. The research is done by means of mathematical modeling of physical processes. It is based on numerical solution of Reynolds equations for chemical components and equations of energy conservation for gaseous and condensed (for canopy) phases. To describe the transfer of energy by radiation the diffusion approximation is used. The boundary-value problem is solved numerically using the method of splitting according to physical processes. In the first stage, the hydrodynamic pattern of flow and distribution of scalar functions are calculated. The system of ordinary differential equations of chemical kinetics obtained as a result of splitting is then integrated. The chemical time step is selected automatically. A discrete analog for equations is obtained by means of the control volume method using the SIMPLE algorithm. Calculation method and program have been check. As a result of mathematical modeling the distributions of temperatures, mass concentrations of components of gaseous phase, volume fractions of components of solid phase, as well as fields of velocity at different instants of time are obtained. It allows to investigate dynamics of forest fire initiation under influence of various external conditions: a) meteorology conditions (air temperature, wind velocity etc.), b) type (various kinds of forest combustible materials) and their state (load, moisture etc.). A great deal of final and intermediate gaseous and dispersed combustion products of forest fuels is known to be exhausted into the atmosphere during forest fires. The knowledge of these kinds of ejections enables a full estimate of the damage from forest fires to be made. The numerical results obtained by the proposed method show satisfactory agreement with experimental measured data.



Optimized Schwarz methods with an overset grid system for the Shallow-Water Equations


Abdessamad Qaddouri (email: Abdessamad.qaddouri@ec.gc.ca)
Recherche en prévision numérique, Meteorological Service of Canada
Co-authors: Jean Côté, Martin Gander and Lahcen Laayouni

The overset grid system nicknamed "Ying-Yang" grid is singularity free and has quasi-uniform grid spacing. It is composed of two identical latitude/longitude orthogonal grid panels that are combined to cover the sphere with partial overlap on their boundaries. The system of Shallow-Water equation (SWEs) is a hyperbolic system at the core of many models of the atmosphere. In this paper SWEs are solved on the Ying-Yang grid by using a semi-implicit and semi-Lagrangian discretization on a staggered mesh. The scalar elliptic equation is solved trough Schwarz-type domain decomposition methods, known as optimized Schwarz methods, which gives better performance than the classical Schwarz methods by using specific Dirichlet- or Robin- type interface conditions.



An Efficient Wavelet-Based Solution of Electromagnetic Field Problems


Yotka Rickard (email: yotka@ece.mcmaster.ca)
McMaster University

This work is an application of numerical analysis to the solution of open problems in computational electromagnetics. More specifically, modeling of integrated antennas by the Method of Moment (MoM) is considered. We propose using an efficient wavelet approximate method to solve the electrical field integral equation (EFIE). The method is applied on irregular grids (triangular sub-domains) and is applicable to arbitrary antenna shapes, which cannot be accomplished by using FFT, for example.

Instead of searching new discretizations, we propose the discrete wavelet transform (DWT) to be applied to the discretized EFIE via the MoM, and then to solve the resulting matrix equation in the wavelet space, followed by the inverse DWT of the solution vector. This approach allows the available software for discretizing EFIE to be used with greater success, especially when solving matrix equations of very large size. The proposed procedure involves O(N^2) operations and leads to a matrix problem, which can be sparsified. The latter opportunity is based on the fact that the MoM matrix is fully populated but has many small off-main-diagonal elements, which can be neglected without losing valuable information. To that goal, suitable threshold techniques are used to create sparse matrix while not compromising the accuracy. This is possible because the wavelet expansions create matrices with condition numbers orders of magnitude better that the conventional MoM. We solve the matrix problem using iterative methods (e.g., by the BiCG method) with suitable preconditioners in O(N^2) operations.

Several numerical examples are considered and it is shown that the CPU time is reduced by a factor of 10 for problems of N > 1000.



Use of near breakdowns in block Arnoldi method to solve Sylvester equations


Miloud Sadkane (email: sadkane@univ-brest.fr)
University of Brest, France
Co-authors: Mickael Robbé

The Sylvester equation

A X + X B = CDT, (1)

where A and B are large square matrices and C and D are rectangular matrices plays a significant role in linear control theory. The most popular methods construct low rank approximate solutions to (1) of the form

Xj = Vj Yj WjT, (2)

where Yj solves a small order Sylvester equation and where Vj and Wj are rectangular matrices whose columns form orthonormal bases of the block Krylov spaces

Kj (A, C ) = range[ C, A C, ..., Aj-1 C ],

Kj (BT, D ) = range[ D, BTD, ..., (BT)j-1 D ].

We propose an algorithm to compute low rank approximations of the type (2), where the matrices Vj and Wj are constructed with a recent version of the block Arnoldi algorithm. A version which takes into account the near breakdowns, i.e., the near singularities in the Krylov matrices

[ C, A C, ..., Aj-1 C ] and [ D, BTD, ..., (BT)j-1 D ].

This leads to matrices Vj and Wj whose columns are selected from the spaces Kj (A, C) and Kj (BT, D) in a way appropriate for convergence.

We show how to detect and use the near breakdowns and discuss the theory and numerical aspects of the proposed algorithm.



Symplectic Householder transformations, an algebraic and geometric approaches


Ahmed Salam (email: Ahmed.Salam@lmpa.univ-littoral.fr)
Laboratoire de Mathématiques Pures et Appliquées,
Université du Littoral-Côte d'Opale, France
Co-authors: Anas ELFAROUK

The aim of this paper is to introduce and to study symplectic Householder transformations, based on an algebraic and geometric approaches. The construction is so that the parallel between orthogonal Householder and symplectic Householder becomes easy and natural. A block form is given. A structure-preserving algorithm for computing SR factorization via symplectic Householder transformations is obtained. Numerical tests are given. Keywords: Skew-symmetric inner product, symplectic geometry, symplectic transvection, symplectic Householder transformation, symplectic SR factorization, structure-preserving methods.

AMS subject classifications. 65F15, 65F50



Dynamically Adapted Inexact Additive Schwarz Preconditioner


Daniel B Szyld (email: szyld@math.temple.edu)
Temple University
Co-authors: Marcus Sarkis, IMPA, Rio de Janeiro, and Worcester Polytechnic Inst.

Consider the solution of A u = f using additive Schwarz (right) preconditioning, i.e., A M-1 w = f, with Mu = w. The preconditioner is M-1 = åi=1p RiT Ai-1 Ri , where Ai = Ri A RiT. We also consider a restricted additive Schwarz preconditioner [1]. When the local problems Ai zi = vi are hard to solve, one considers approximations to their solution, and uses some inexact local solver [A\tilde]i, for example, one may use a (local or inner) Krylov method with some (fixed) inner stopping criterion. We apply the recently developed theory of inexact Krylov subspace methods [2], or in this case of inexact preconditioning, by which we allow the inexactness of the local solver to change from step to step. We provide a computable criterion to allow the local (or inner) stopping criterion to be relaxed as the (outer) iteration progresses while maintaining overall convergence to the desired (outer) tolerance. Numerical examples illustrate the potential of this dynamically chosen local inexact solvers.

References:
[1] X.-C. Cai and M. Sarkis. A restricted additive Schwarz preconditioner for general sparse linear systems. SIAM J. Sci. Comput., 21:792-797, 1999.
[2] V. Simoncini and D. B. Szyld. Theory of inexact Krylov subspace methods and applications to scientific computing. SIAM J. Sci. Comput., 25:454-477, 2003.



Convergence of Krylov subspace methods when using non-orthogonal bases


Daniel B Szyld (email: szyld@math.temple.edu)
Department of Mathematics, Temple University, Philadelphia
Co-authors: Valeria Simocini, Dept. Matematica, Univerista di Bologna, and
IMATI-CNR, Pavia, Italy

There are many examples where non-orthogonality of a basis for Krylov subspace methods arises naturally. These methods usually require less storage or computational effort per iteration than methods using an orthonormal basis (öptimal" methods), but the convergence may be "delayed". Truncated Krylov subspace methods and other examples of "non-optimal" methods have been shown to converge in many situations, often with small delay, but not in others. We explore the question of what is the effect of having a non-optimal basis. We prove certain identities for the relative residual gap, i.e., the relative difference between the residuals of the optimal and non-optimal methods. These identities and related bounds provide insight into when the delay is small and convergence is achieved. Further understanding is gained by using a general theory of superlinear convergence recently developed. Our analysis confirms the observed fact that the orthogonality of the basis is not important, only the need to maintain linear independence is. Numerical examples illustrate our theoretical results.



A modified Newton interior point method for nonlinear programming


G. Tanoh (email: gtanoh@dim.uchile.cl)
Centro de Modelamiento Matematico, Universidad de Chile

We analyze an iterative method for minimizing general nonlinear function with equality and inequality constraints by interior point method. Applying the barrier penalty technique, we solve a sequence of unconstrained minimization problems with primal-dual merit function. We develop an iterative method to compute the Newton steps. Descent direction and direction of negative curvature are found by using alternatively Conjugate Gradient and Lanczos methods. Various computational enhancements to improve these directions are provided. Our strategy is well adapted to deal with indefiniteness and can be applied to large scale nonlinear optimization problems. Promising preliminary numerical results are presented.



Screening of the patients with arrhythmia based on fractal dimension
of heart rate variability by an artificial neural network


Iman Tavassoly (email: iman_tavassoly@hotmail.com)
Bioinformatics and Biomathematics Unit, SRC,
Mazandaran University of Medical Sciences, Sary, IRAN
Co-authors: Omid Tavassoly, Mohammad soltany

Objectives: The heart rate variability has been shown to have fractal behaviour during the time. Fractals are self similar shapes with fractional dimension and we can use them to find order in the biologic systems which seem to have non predictive and irregular behaviour. The exact relationship between fractal behaviour of heart rate variability and the condition of the cardiovascular system and its changes with age and diseases have not well studied.

We used an artificial neural network for screening the patient with arrhythmia as a pathologic state based on fractal dimension of heart rate variability.

Methods: We used 200 digital rhythmograms which included 100 healthy subjects and 100 patients with arrhythmia. We calculated the fractal dimension of heart rate variability using Nlyser software (version 3.5). We designed an artificial neural network by Neurosollutions software (version 4) and we used the age, gender and fractal dimension of heart rate variability of 150 rhythmograms (75 healthy rhythmograms and 75 ones with arrhythmia) as input data and the condition of being healthy or having arrhythmia as target for training the artificial neural network. We used the data of 20 rhythmograms (10 healthy rhythmograms and 10 ones with arrhythmia) as validation set. Finally after training the network it was tested using the data of 30 rhythmograms.

Results: The designed artificial neural network could screen the healthy persons from arrhythmia patients based on the fractal dimension of heart rate variability with 0.2 mean squared errors. The network was 100% successful in finding healthy persons and it was 90.14% successful in finding arrhythmia patients. It has fault for discriminating the arrhythmia patients in 98.6% of cases.

Conclusions: Fractal dimension which is used as a factor showing fractal characteristics of heart rate variability can be used as a diagnostic tool to screen the patients with arrhythmia by artificial neural networks. The results will depend on the size of our training set for artificial neural network. If we develop the training set the network will be more successful in screening the patients.



Two Uses for Updating the Partial Singular Value Decomposition
in Latent Semantic Indexing


Jane E. Tougas (email: tougas@cs.dal.ca)
Dalhousie University
Co-authors: Henry Stern (Dalhousie University), Raymond J.Spiteri (University of Saskatchewan)

Latent Semantic Indexing (LSI) is an information retrieval (IR) method that unites IR with numerical linear algebra by representing a dataset as a term-document matrix. Because of the tremendous size of modern databases, such a term-by-document matrix can potentially be very large, with millions of entries. The partial singular value decomposition (PSVD) is a matrix factorization that captures the salient features of a matrix, while using much less storage. We look at two challenges posed by this PSVD data compression process in LSI. Traditional methods of computing the PSVD are very expensive; most of the processing time in LSI is spent in calculating the PSVD of the term-by-document matrix [1]. Given the potentially huge size of this matrix, the first challenge is calculating the PSVD efficiently, in terms of computational and memory requirements. The second challenge is efficiently updating the PSVD when the matrix is altered slightly. In a rapidly expanding environment, such as the Internet, the term-by-document matrix is altered often as new documents and terms are added. Updating the PSVD of this matrix is much more efficient than recalculating it after each change. We investigate the use of SVD updating methods proposed by Zha and Simon [2] to meet both of these challenges. These updating methods for adding documents or terms to the term-by-document matrix use the previously calculated PSVD of the term-by-document matrix to compute the PSVD of the updated matrix. Results will be presented illustrating that updating in this manner provides a tremendous saving in computation time, with little or no significant reduction in accuracy. Having shown that the updating methods are both fast and accurate, an algorithm for iteratively computing the PSVD of an extremely large matrix using the document updating method will then be presented. This iterative method, suggested by Zha and Zhang [3], provides a means of calculating the PSVD for matrices so large that the computation would be infeasible using traditional methods. Again, results will be given showing that this method provides savings in computational time and memory resources without compromising the accuracy of the results.

[1] M.W. Berry, S.T. Dumais, and G.W. O'Brien. Using Linear Algebra for Intelligent Information Retrieval. SIAM Rev., 37(4):573-595, 1995.
[2] H. Zha and H.D. Simon. Timely Communication on Updating Problems in Latent Semantic Indexing. SIAM J. Sci. Comput., 21(2):782-791, 1999.
[3] H. Zha and Z. Zhang. Matrices with Low-rank-plus-shift Structure: Partial SVD and Latent Semantic Indexing. SIAM J. Matrix Anal. Appl., 21(2):522-536, 1999.



Trust region GMRES methods for systems of nonlinear equations


Ming-yan Wang (email: mywangworld@yahoo.com)
Institute of Computational Mathematics and ScientificEngineering Computing

In this paper we propose a new trust region method for systems of nonlinear equarions. We consider this problem only in a subspace. Global convergece is got using GMRES strategy. Numerical experiences show our algorithm is efficient.



A high-performance method for the biharmonic Dirichlet problem


Jingrui Zhang (email: jingrui@cs.toronto.edu)
Department of Computer Science, University of Toronto
Co-authors: Christina Christara

The biharmonic Dirichlet problem appears in many applications. Therefore, there is a lot of interest in developing high-performance methods for its solution. Two important components of performance are accuracy and efficiency. In this research, we combine a sixth order discretization method and a Fast Fourier Transform (FFT)-based solver for the solution of the biharmonic Dirichlet problem.

The discretization of the problem uses the collocation methodology and quartic splines. Quartic Spline Collocation (Quartic SC) methods that produce sixth order approximations have been recently developed for the solution of fourth order Boundary Value Problems (BVPs). In this paper, we apply an adjustment to the quartic spline basis functions so that the Quartic SC matrix has more favorable properties. We study the properties of two Bi-Quartic SC matrices, more specifically, one that corresponds to Dirichlet and Neumann conditions (the standard harmonic problem) and another that corresponds to Dirichlet and second derivative conditions. The latter is solvable by FFTs and used as preconditioner for the first which is solved by the preconditioned GMRES method. Numerical experiments verify the effectiveness of the solver and preconditioner.


Minisymposia talks


Minisymposium: Progress in Eigenvalue Methods (M1)

Jacobi-Davidson techniques for the Hamiltonian eigenvalue problem


Michiel Hochstenbach (email: michiel.hochstenbach@case.edu)
Department of Mathematics
Case Western Reserve University

Structure preserving methods for eigenvalue problems have received quite some attention recently. We present several alternative numerical methods for the solution of an eigenvalue problem with Hamiltonian structure. These approaches are based on Jacobi-Davidson type techniques.



New ideas for computing complex eigenvalues of an
asymmetric matrix applied to metastable states


Tucker Carrington (email: Tucker.Carrington@umontreal.ca)
Department of Chemistry
Université de Montréal

Lifetimes of metastable states may be determined by computing complex eigenvalues of a complex-symmetric matrix. One way to do this is to use the complex-symmetric version of Cullum and Willoughby's non-orthogonalised Lanczos procedure [1]. It is also possible use shift and invert complex-symmetric Lanczos to compute eigenvalues window by window [2]. Both of these approaches require doing complex matrix-vector products. In this talk I shall outline a new method for obtaining complex eigenvalues from real matrix-vector products. The complex-symmetric eigenvalue problem is written as a quadratic real eigenvalue problem which is linearized [3]. We attempted to solve the linearized (asymmetric) problem using the symmetric indefinite Lanczos method [4] but discovered that the symmetric indefinite Lanczos method works poorly for the purpose of computing many eigenpairs. Instead, we use the standard two-sided Lanczos algorithm (exploiting the symmetry of the eigenvalue problem to reduce the number of matrix-vector products by a factor of two) to compute approximate eigenvectors. It is not possible to get very accurate eigenvectors. We use groups of approximate eigenvectors as a basis and compute extremely accurate complex eigenvalues.

[1] Jane Cullum, Wolfgang Kerner, Ralph Willoughby A generalized nonsymmetric Lanczos procedure, Comput. Phys. Comm., 53 (1989)
[2] Bill Poirier and T. Carrington A preconditioned inexact spectral transform method for calculating resonance energies and widths, as applied to HCO, J. Chem. Phys. 116 1215-1227 (2002).
[3] A. Neumaier and V.A. Mandelshtam, Further generalization and numerical implementation of pseudo-time Schroedinger equations for quantum scattering calculations, J. Theor. Comp. Chem. 1, 1 (2002).
[4] Templates for the Solution of Algebraic Eigenvalue Problems Zhaojun Bai, James Demmel, Jack Dongarra, Axel Ruhe, and Henk van der Vorst



On the Computation of Optical Lasing Modes of Axisymmetric VCSEL Devices


Peter Arbenz (email: arbenz@inf.ethz.ch)
Institute of Computational Science
ETH Zurich, Switzerland

The computation of optical modes inside axisymmetric cavity resonators with a general spatial permittivity profile is a formidable computational task. In order to avoid spurious modes the vector Helmholtz equations are discretised by a mixed finite element approach. We formulate the method for first and second order Nedelec edge and Lagrange nodal elements. We discuss how to accurately compute the element matrices and to solve the resulting large sparse complex symmetric eigenvalue problems. We validate our approach by three numerical examples that contain varying material parameters and absorbing boundary conditions (ABC).



Locking issues for finding a large number of eigenvalues of symmetric matrices


Andreas Stathopoulos (email: andreas@cs.wm.edu)
Department of Computer Science
College of William and Mary

Iterative methods that find many eigenvalues of large, sparse, symmetric matrices have to deal with the converged eigenpairs. Typically some form of deflation is used to keep them from reappearing, and to speed up convergence. One of the more stable forms of deflation is known as locking. This usually means that converged eigenvectors Q are removed from the search space V, and all current and subsequent vectors in V are orthogonalized against Q. The alternative, sometimes called soft-locking, is to keep Q in the basis, making sure the orthogonality of V is maintained. Converged eigenvectors must be flagged as converged so that they are not targeted for further improvement by the algorithm. In soft-locking, the converged eigenvectors participate in the Rayleigh-Ritz projection, hence they improve at every step. However, the basis becomes large, and restarting and Rayleigh Ritz projection become too expensive.

We show that, in some cases, locking is susceptible to numerical instability, especially when the eigenvectors are not computed close to full accuracy. In particular, if k eigenvectors Q have converged and are locked with residual tolerance e, any subspace method (like Lanczos or Jacobi-Davidson) may be unable to converge to the k+1 eigenvector to full accuracy e. The reason is that the k+1 eigenvector of A is not the same as the lowest eigenvector of (I-QQT)A(I-QQT). The same reason may cause the correction equation of Jacobi-Davidson type methods to converge without having produced the appropriate correction. We explore possible solutions to these problems.


Minisymposium: Recent Advances in Multilevel Methods I (M2)


Adaptive Algebraic Multigrid


Scott MacLachlan (email: Scott.MacLachlan@colorado.edu)
Department of Applied Mathematics
University of Colorado at Boulder

Numerical simulation of physical processes is often constrained by our ability to solve the complex linear systems at the core of the computation. Classical geometric and algebraic multigrid methods rely on (implicit) assumptions about the character of these matrices in order to develop appropriately complementary coarse-grid correction processes for a given relaxation scheme. The aim of the adaptive multigrid framework is to reduce the restrictions imposed by such assumptions, thus allowing for efficient black-box multigrid solution of a wider class of problems.

There are, however, many challenges in altogether removing the reliance on assumptions about the errors left after relaxation. In this talk, we discuss work to date on an adaptive AMG algorithm that chooses the components of the coarse-grid correction cycle based on automated analysis of the performance of relaxation. Fundamental measures of the need for and quality of a coarse-grid correction will be discussed, along with related techniques for choosing coarse grids and interpolation operators. We will also discuss the (lack of) computability of these ideal measures, and suggest cost-efficient alternatives.

 



Self-adaptative Multigrid via Subcycling


Tim Chartier (email: tichartier@davidson.edu)
Department of Mathematics
Davidson College
Co-authors: Edmond Chow

Considerable efforts in recent multigrid research have concentrated on algebraic multigrid schemes. A vital aspect of this work is uncovering algebraically smooth error modes in order to construct effective multigrid components. Many existing algebraic multigrid algorithms rely on assumptions regarding the nature of algebraic smoothness. For example, a common assumption is that smooth error is essentially constant along `strong connections'. Performance can degrade as smooth error for a problem differs from this assumption. Through tests on the homogeneous problem (Ax = 0) adaptive multigrid methods expose algebraically smooth error.

The method presented in this talk uses relaxation and subcycling on complementary grids as an evaluative tool in correcting multigrid cycling. Each complementary grid is constructed with the intent of dampening a subset of the basis of algebraically smooth error. The particular implementation of this framework manages smooth error in a manner analogous to spectral AMGe. Numerical results will be included.



Adaptive Algebraic Multigrid in Quantum Chromodynamics


James Brannick (email: brannick@newton.colorado.edu)
Co-authors: Marian Brezina, Scott MacLachlan, Tom Manteuffel,
Steve McCormick, John Ruge
Department of Applied Mathematics
University of Colorado at Boulder

Standard algebraic multigrid methods assume explicit knowledge of so-called algebraically-smooth or near-kernel components, which loosely speaking are errors that give relatively small residuals. Tyically, these methods automatically generate a sequence of coarse problems under the assumption that the near-kernel is locally constant. The difficulty in applying algebraic multigrid to lattice QCD is that the near-kernel components can be far from constant, often exhibiting little or no apparent smoothness. In fact, the local character of these components appears to be random, depending on the randomness of the so-called "gauge" group. Hence, no apriori knowledge of the local character of the near-kernel is readily available.

This talk develops an adaptive algebraic multigrid (AMG) preconditioner suitable for the linear systems arising in lattice QCD. The method is a recently developed extension of smoothed aggregation, aSA [1], in which good convergence properties are achieved in situations where explicit knowledge of the near-kernel components may not be available. This extension is accomplished using the method itself to determine near-kernel components automatically, by applying it carefully to the homogeneous matrix equation, Ax=0. The coarsening process is modified to use and improve the computed components. Preliminary results with model 2D QCD problems suggest that this approach may yield optimal multigrid-like performance that is uniform in matrix dimension and gauge-group randomness.

[1] Brezina, M., Falgout, R., MacLachlan, S., Manteuffel, T., S. McCormick, S., Ruge, J.: Adaptive Smoothed Aggregation (aSA). SIAM J. Sci. Comp., 25 (2004), pp. 1806-1920.



Algebraic Multigrid (AMG) Preconditioning for Higher-Order Finite Elements


Luke Olson (email: lolson@dam.brown.edu)
Division of Applied Mathematics
Brown University

In this talk we consider two related approaches for solving linear systems that arise from a higher-order finite element discretization of elliptic partial differential equations. The first approach explores direct application of an algebraic-based multigrid method (AMG) to iteratively solve the linear systems that result from higher-order discretizations. While the choice of basis used on the discretization has a significant impact on the performance of the solver, results indicate that AMG is capable of solving operators from both Poisson's equation and a first-order system least-squares (FOSLS) formulation of Stoke's equation in a scalable manner, nearly independent of basis order, p. The second approach incorporates preconditioning based on a bilinear finite element mesh overlaying the entire set of degrees of freedom in the higher-order scheme. AMG is applied to the operator based on the bilinear finite element discretization and is used as a preconditioner in a conjugate gradient (CG) iteration to solve the algebraic system derived from the high-order discretization. This approach is also nearly independent of p. We present several numerical examples that support each method and discuss the computational advantages of the preconditioning implementation.


Minisymposium: Recent Advances in Multilevel Methods II (M3)


Reducing Complexity in Algebraic Multigrid


Hans De Sterck (email: hdesterck@uwaterloo.ca)
Department of Applied Mathematics, University of Waterloo (Canada)

Algebraic multigrid (AMG) is a very efficient iterative solver and preconditioner for large unstructured linear systems. Traditional coarsening schemes for AMG can, however, lead to computational complexity growth as problem size increases, resulting in increased memory use and execution time, and diminished scalability. We propose new parallel AMG coarsening schemes, that are based on solely enforcing a maximum independent set property, resulting in sparser coarse grids. The new coarsening techniques remedy memory and execution time complexity growth for various large three-dimensional (3D) problems. If used within AMG as a preconditioner for Krylov subspace methods, the resulting iterative methods tend to converge fast. For some difficult problems, however, these methods don't converge well enough, and improved interpolation is necessary in order to obtain scalable methods.

*This is joint work with Ulrike Yang. This work was performed under the auspices of the U.S. Department of Energy by the University of California Lawrence Livermore National Laboratory under contract number W-7405-Eng-48 and subcontract number B545391.




On parallel algebraic multigrid preconditioners for systems of PDEs


Ulrike Meier Yang (email: umyang@llnl.gov)
Center for Applied Scientific Computing
Lawrence Livermore National Laboratory (USA)

Algebraic multigrid (AMG) is a very efficient, scalable algorithm for solving large linear systems on unstructured grids. When solving linear systems derived from systems of partial differential equations (PDEs) often a different approach is required than for those derived from a scalar PDE. There are mainly two approaches, the function approach (also known as the "unknown" approach), and the nodal or "point" approach. The function approach defines coarsening and interpolation separately for each function. The nodal approach uses AMG in a block manner, where all variables that correspond to the same grid node are coarsened, interpolated and relaxed together. While the function approach is much easier to implement and often more efficient, there are problems for which this approach is not sufficient and the more expensive nodal approach is needed.

In this paper we present several parallel implementations of both approaches using various coarsening schemes and interpolation operators. Advantages and disadvantages of both approaches are discussed, and numerical results are presented.

*This work was performed under the auspices of the U.S. Department of Energy by University of California Lawrence Livermore National Laboratory under contract number W-7405-Eng-48.



Scalability advances in algebraic multigrid for Maxwell's Equations


Jonathan Hu (email: jhu@ca.sandia.gov)
Department of Computational Mathematics and Algorithms
Sandia National Laboratories (USA)

In this talk, we present algorithmic and parallel performance improvements to Sandia's algebraic multigrid solver for the eddy current equations. The motivation for these improvements is a 3600 processor, multiphysics Zpinch simulation on Sandia's parallel machine, Janus. We begin with a very brief overview of our AMG solver for the eddy current equations. We then discuss the various challenges that large-scale Zpinch simulation present to our particular AMG solver, such as operator load-balancing and memory constraints, as well as our solutions to these issues. We finish with numerical results that show vastly improved scalability in the AMG solver.



A Multilevel Method for Image Registration


Eldad Haber (email: haber@mathcs.emory.edu)
Department of Mathematics and Computer Science
Emory University (USA)

In this paper we introduce a new framework for image registration. Our formulation is based on consistent discretization of the optimization problem coupled with a multigrid solution of the linear system which evolve in a Gauss-Newton iteration. We show that our discretization is h-elliptic independent of parameter choice and therefore a simple multigrid implementation can be used. To overcome potential large nonlinearities and to further speed up computation, we use a multilevel continuation technique. We demonstrate the efficiency of our method on a realistic highly nonlinear registration problem.


Minisymposium: Preconditioning linear and nonlinear iterations I (M4)

An analysis of partitioned nonlinear systems


Dhavide Aruliah (email: dhavide.aruliah@uoit.ca)
University of Ontario Institute of Technology (Canada)

I will present some recent work looking at iterative methods for a fairly general family of partitioned systems of nonlinear equations that seems to occur in numerous applications (e.g., discretisation of elliptic partial differential equations). In such systems, various analytic problem formulations exist based on a natural partitioning of the unknowns or introduction of auxiliary unknowns. The existence of distinct analytic formulations prompts the question of which analytic system should be solved. While the problems are mathematically equivalent, once fed into appropriate variants of Newton's method, distinct discrete dynamical systems result. These distinct formulations suggest different preconditioning strategies for inner iterations of inexact Newton methods in turn. I will present some examples of such systems and some suggestions of how to choose between alternative equivalent formulations of nonlinear systems of equations that arise in practical applications.



Multi-Preconditioned Conjugate Gradient


Robert Bridson (email: rbridson@cs.ubc.ca)
University of British Columbia (Canada)

This talk introduces a variation on standard preconditioned Krylov methods for solving linear systems, where we extend our search space using several preconditioners instead of just one. Our multi-preconditioned conjugate gradient (MPCG) solver can prove useful for problems where it is not clear which preconditioning approach is the most effective. At the very least, MPCG can robustly identify which to select for a faster PCG iteration, but in some cases can automatically exploit synergy between the preconditioners for asymptotically faster convergence. Numerical experiments for bending problems and domain decomposition highlight the capabilities of MPCG. Unfortunately, similar to flexible methods, the recurrence relation is not short in general. However, in at least one case (where the two preconditioners sum to the original matrix, as in 2D ADI) I will prove the recurrence is short. I will also explore the possibilities for truncation in other cases.



Preconditioned Newton-Krylov iterations for large-scale continuation


Homer Walker (email: walker@wpi.edu)
Worcester Polytechnic Institute (USA)

We consider the application of preconditioned Krylov subspace methods to continuation problems, in which the object is to track a solution curve as a continuation parameter varies. Continuation methods are typically of predictor-corrector form, in which, at each step along the curve, Newton-like corrector iterations are used to return to the curve from an initial "predicted" point. The linear system that characterizes corrector steps is underdetermined, and we discuss the adaptation of preconditioned Krylov subspace methods to solving this system. The focus is on a particular approach that preserves problem conditioning, allows the use of preconditioners constructed for the fixed-parameter case, and has certain other advantages. The effectiveness of the approach is shown in numerical experiments on large-scale problems.



Constraint Preconditioning for saddle-point systems


Andy Wathen (email: wathen@comlab.ox.ac.uk)
Oxford University (Britain)

Saddle-point systems arise in any problem with constraints: the incompressibility condition is a constraint on the momentum equations for incompressible flows, coupling in multiphysics problems is solving one problem subject to the other being satisfied (and conversly by duality) and throughout Optimization constrained problems are ubiqitous. In this talk we will describe preconditioners for linear saddle-point systems which preserve the constraints. Such `Constraint Preconditioners' have attractive theoretical properties which we will describe.


Minisymposium: Preconditioning linear and nonlinear iterations II (M5)

Approximate factorisation constraint preconditioners


Sue Dollar (email: hsd@comlab.ox.ac.uk)
Oxford University (Britain)

The problem of finding good preconditioners for the numerical solution of a certain important class of indefinite systems is considered. These systems are of a saddle point structure $$ \left[\begin{array}{cc} A & B^T \\ B & 0 \\ \end{array}\right] \left[\begin{array}{c} x \\ y \\ \end{array}\right] = \left[\begin{array}{c} c \\ d \\ \end{array}\right], $$ where $A$ is a $n$ by $n$ real, symmetric matrix and $B$ is a $m$ by $n$ real matrix $(m\leq n)$. The idea of using ``constraint preconditioners'' that have a specific 2 by 2 block structure was introduced by Keller, Gould and Wathen (SIMAX vol 21, 2000). They showed that the preconditioned system have some favourable spectral properties. However, even with simple choices for these constraint preconditioners, the factoring of the preconditioner, $K$, may become infeasible for large problems. Suppose that we use a constraint preconditioner of the form $$K = P B P^T,$$ where solutions with each of the matrices $P$, $B$ and $P^T$ are easily obtained. In particular, rather than obtaining $P$ and $B$ from a given $K$, $K$ is derived implicitly from specially chosen $P$ and $B$. We shall investigate how these ``approximate factorization constraint preconditioners'' can be successfully applied to saddle-point problems of the form considered by Keller, Gould and Wathen.



A block diagonal preconditioner for saddle point linear systems
arising from mixed finite element formulation of
time-harmonic Maxwell's equations


Chen Greif (email: greif@cs.ubc.ca)
University of British Columbia (Canada)
Co-authors: Dominik Schoetzau

The mixed finite element formulation of the time-harmonic Maxwell's equation motivates the preconditioners considered in this talk. The saddle point linear system that represents the discretization has a (1, 1) block with high nullity. A discrete divergence-free property and spectral equivalence results that we derive are exploited, and a special block diagonal preconditioner is proposed for the problem. We also present a class of algebraic preconditioners that do not rely on computation of Schur complements or their approximations and may work well for a large class of saddle point systems.



All-at-once inversion of time domain electromagnetic data


Eldad Haber (email: haber@mathcs.emory.edu)
Department of Mathematics and Computer Science, Emory University (USA)

We develop an inversion methodology for 3D electromagnetic data when the forward model consists of Maxwell's equations in which the permeability is constant but electrical conductivity can be highly discontinuous. The goal of the inversion is to recover the conductivity given measurements of the electric and/or magnetic fields. A standard Tikhonov regularization is incorporated and we use an inexact, all-at-once methodology, solving the forward problem and the inverse problem simultaneously in one iterative process. This approach allows development of highly efficient algorithms. Here we present the basic methodology for the all-at-once approach. We then briefly review the time domain forward modelling equations, develop the linearized Gauss-Newton system of equations to be solved at each iteration, and show how these equations can be solved using a preconditioned conjugate gradient solution. We invert a synthetic data set generated from a surface loop transmitter and receivers in four boreholes.


Minisymposium: Combinatorial and Computational Aspects
of the Monomer-Dimer Problem (M6)


Combinatorial and Computational Aspects of the Monomer-Dimer Problem


Shmuel Friedland (e-mail: friedlan@uic.edu)
Co-author: Uri N. Peled (email: uripeled@uic.edu)
Department of Mathematics, Statistics, and Computer Science,
University of Illinois at Chicago, Chicago, Illinois 60607-7045, USA

The exponential growth rate h of the number of configurations on a multi-dimensional grid arises in many fields. It is called topological entropy in mathematics, entropy in physics and multi-dimensional capacity in engineering. In the 1-dimensional case the entropy equal to the spectral radius of a certain digraph. There are very few 2-dimensional models where the value of h is known in a closed form. In this talk we first survey the theory of the computation of h by using lower and upper bounds that converge to h. As a demonstration of these techniques, we compute the topological entropy of the monomer-dimer covers of the 2-dimensional grid to 8 decimal digits and of the 3-dimensional grid with an error smaller than 1.35

Second we discuss a more delicate problem to compute the exponential growth rate h(p), with prescribed densities of configurations. For that purpose we introduce the notion of pressure P, which is a convex function in d-dimensional space. We can compute its values using lower and upper values which converge to P. The values of partial derivatives of P give the value of h(p), which can be estimated from below and above.

As an application we compute the h(p) for the monomer-dimer configurations in 2-dimensional grid. First we let p be the density of the dimers. Then our lower and upper bounds confirm the heuristic computations of Baxter 1968. Second we show that we can compute the h(x,y), where x and y are the densities of the dimers in the direction X and Y in the plane.



Lower Bounds for Partial Matching in Regular Bipartite Graphs with
Application to the Monomer-Dimer Problem


Elliot Krop (email: ekrop1@math.uic.edu)
Co-author: Shmuel Friedland (email: friedlan@uic.edu)
Department of Mathematics, Statistics, and Computer Science,
University of Illinois at Chicago, Chicago, Illinois 60607-7045, USA

The exponential growth h(p) of the monomer dimer configurations in the d-dimensional lattice, with prescribed density p of dimers, is equivalent to the exponential growth of the number of partial matching in the bipartite graphs corresponding to d-dimensional lattice tori with even length sides. The number of this partial matching is equal to a sum of subpermanents of certain orders in the corresponding integer value matrix.

In the work of Friedland and Peled it was shown that h(p) can be bounded below by a(p), using the sharp lower bound on the subpermanents of doubly stochastic matrices. As in the work of Schrijver on permanents of 0-1 matrices, corresponding to regular bipartite graphs, it is plausible that one can improve the lower bound a(p).

In this talk we present a conjectured lower bound b(p) of h(p), which indeed is better than a(p). We give a probabilistic motivation for the form of b(p).

We give numerical examples which test the conjectured function b(p). This is done by using the notion of pressure introduced by Friedland and Peled. Our case is relatively simple case of this theory, since it is reduced to the one dimensional subshift of finite type. Our computation are based on calculations of spectral radii of certain transfer matrices and their derivatives.


Minisyposium: Multigrid and Optimized Schwarz Preconditioners
for High-Order Finite-Elements (M7)


Optimized Multiplicative, Additive and Restricted Additive Schwarz Preconditioning


Amik St-Cyr (email: amik@ucar.edu)
National Center for Atmospheric Research
1850 Table Mesa Drive, Boulder, Colorado, 80305
Co-authors: Martin J. Gander (University of Geneva, Switzerland) and
Stephen J. Thomas (National Center for Atmospheric Research)

We demonstrate that a small modification of the multiplicative, additive and restricted additive Schwarz preconditioner at the algebraic level, motivated by optimized Schwarz methods defined at the continuous level, leads to a significant reduction in the iteration count of the iterative solver. In the context of optimized Schwarz methods this approach permits the use of low order subdomain preconditioners for high order spectral element discretizations, while retaining spectral accuracy in the solution. Numerical experiments using finite difference and spectral element discretizations of the modified Helmholtz problem illustrate the effectiveness of the new approach.



A Non-Overlapping Schwarz Method with
Higher-Order Transmission Conditions for
Time-Harmonic Maxwell Problems


Marinos N. Vouvakis (email:vouvakis.1@osu.edu)
The Ohio State University
Co-authors: Seung-Cheol Lee and Jin-Fa Lee

During the years Domain Decomposition (DD) methods have been proven very efficient and effective numerical techniques for the analysis of Partial Differential Equations. Among other advantages, it suffices to stress their parallelization ability, scalability, ability to systematically couple different numerical methods into hybrid schemes, efficient exploitation of geometrical redundancies and symmetries, and relaxing meshing and adaptive meshing strategies.

In the hart of every DD algorithm is the transmission condition (TC). This is the boundary condition imposed on the interfaces of subdomains to facilitate the information passing between domains. As it is natural, the choice of TC can lead to a successful or a divergent DD algorithm. Appropriately designed TCs have been used to optimize and accelerate convergence of the DD iteration. Such accelerations have been reported for elliptic and Helmholtz problems (M. Gander F. Magoules and F. and Nataf, SIAM SISC., 24, 38-50, 2002), by adopting Robin-type TCs with complex Robin constants.

Recently the authors (S.-C. Lee, M. N. Vouvakis, and J.-F. Lee, J. of Comput. Phy., 203, 1-21, 2005) have proven in the continuous level, that the situation for EM problems, that involve the vector curl-curl operator, is more complicated. With this insight developed in previous work and the powerful Fourier analysis machinery, new higher-order transmission conditions have been developed. Such TCs grand superior convergence properties since they damp both propagative and evanescent error fields.

We intend to present a variant of the Despres-Schwarz Algorithm for time-harmonic electromagnetic problems. In this talk new results for a non-conforming DD method based on higher-order TCs will be given. Convergence curves and computational times will be presented in order to demonstrate the power of the method. Results on real-world challenging radiation and scattering problems such as very large antenna arrays, antenna arrays mounted on large platforms, and antenna arrays in the presence of hybrid radomes and EBGs will be presented.



A Higher-Order Discontinuous Galerkin Discretization of the Unsteady Stokes Problem


Khosro Shahbazi (email: shahbazi@mie.utoronto.ca)
Mechanical Engineering, University of Toronto
Co-authors: Paul Fischer and C. Ross Ethier

Discontinuous Galerkin methods for the Stokes problem with equal and mixed orders of polynomial approximation for the velocity and pressure have been recently developed (see e.g., [P. Hansbo and M. Larson, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 1895-1905] and [D. Schotzau, C. Schwab, and A. Toselli, Math. Models Methods Appl. Sci., 13 (2003), pp. 1413-1436]). However, corresponding efficient numerical solution procedures for the unsteady Stokes problem have not yet been proposed. Our goal is to propose and test such a scheme as the foundation of an unsteady incompressible Navier-Stokes solver.

We have adopted the approximate algebraic splitting scheme (which has been successfully used in the context of the continuous Galerkin method) for the nodal high-order discontinuous Galerkin discretization of the unsteady Stokes problem. Among the available discontinuous Galerkin methods for treating the second order Laplacian, we use the interior penalty method introduced in [D. Arnold, SIAM J. Numer. Anal., 19 (1982), pp. 742-762]. The high-order nodal bases, defined over a standard tetrahedron, consist of Lagrange polynomials calculated based on the nodal set reported in [J. Hesthaven and C. Teng, SIAM J. Sci. Comput., 21 (2000), pp. 2352-2380].

In the proposed splitting scheme, a d-dimensional positive definite Hemholtz problem is solved for the velocity vector and an approximate consistent Poisson equation is solved for the pressure at each time step. For a typical case of a small viscosity and a small time step size, the Helmholtz problem is effectively preconditioned with a block diagonal preconditioner. The pressure system, on the other hand, requires a more sophisticated preconditioner. For this system, we devise a two-level additive Schwarz preconditioner with a local fine part and a global coarse part. The local part corresponds to the summation of the restrictions of the discretized consistent Poisson operator on each element of the triangulation. The global coarse part corresponds to the discretization of the Poisson operator (which is spectrally equivalent to the consistent Poisson operator, but has higher sparsity) on the same mesh but with lower approximation order k=1. We will present data showing the good performance of this preconditioner with respect to the different approximation orders.


Student papers' talks

Performance Optimization and Modeling of Blocked Sparse Kernels


Alfredo Buttari (email: buttari@cs.utk.edu),
Victor Eijkhout, Julien Langou and Salvatore Filippone
Tor Vergata University and University of Tennessee

We present a method for automatically selecting optimal implementations of sparse matrix-vector operations. Our software `AcCELS' (Accelerated Compress-storage Elements for Linear Solvers) involves a setup phase that probes machine characteristics, and a run-time phase where stored characteristics are combined with a measure of the actual sparse matrix to find the optimal kernel implementation. We present a performance model that is shown to be accurate over a large range of matrices.



Recovery Patterns for Iterative Methods in a Parallel Unstable Environment


Zizhong Chen (email: zchen@cs.utk.edu)
G. Bosilca, Z. Chen, J. Dongarra and J. Langou
University of Tennessee

A simple checkpoint-free fault-tolerant scheme for parallel iterative methods is given. Assuming that when one processor fails, all its data is lost and the system is recovered with a new processor, this scheme computes a new approximate solution from the data of the non-failed system. The iterative method is then restarted from this new vector. The main advantage of this technique over standard checkpoint is that there is no extra computation added in the iterative solver. In particular, if no failure occurs, the fault-tolerant application is the same as the original application. The main drawback is that the convergence after failure of the method is no longer the same as the original method. In this paper, we present this recovery technique as well as some implementations of checkpoints in iterative methods. Finally, experiments are presented to compare the two techniques. The fault tolerant MPI library is the FT-MPI library. Iterative linear solvers are considered. (Iterative eigensolvers are considered in [*]).

[*] George Bosilca, Zizhong Chen, Jack Dongarra, and Julien Langou. Recovery patterns for iterative methods in a parallel unstable environment. Technical Report UT-CS-04-538, Department of Computer Science, the University of Tennessee, December 2004.



Public-key Cryptosystems Based on Extended Chebyshev Polynomials


Wang Dahu (email: wangdahu-69@126.com) and Wei Xueye
Beijing jiaotong University, China

Chebyshev polynomials have been recently proposed for designing public-key cryptosystems, which is shown insecure due to the inherent periodicity trait of trigonometric function. An attack can easily recover the corresponding plaintext from a given ciphertext. In the paper, semi-group property based on extended Chebyshev polynomials is used, and two schemes based on such polynomials can avoid this attack and improve security. The first algorithm works on real number and is quite efficient. The second is defined over finite fields and is as secure as RSA. Moreover key agreement, digital signature and entity authentication schemes based on the algorithms are given respectively. In the end, a single consolidate formula used in public-key cryptosystems is concluded.

Keywords: Extended Chebyshev polynomials, Security, Public-key encryptosystem, Entity authentication, Key exchange, Digital signature



Extending constraint preconditioners for saddle point problems


Sue Dollar (email: sue.dollar@comlab.ox.ac.uk)
Oxford University

The problem of finding good preconditioners for the numerical solution of a certain important class of indefinite linear systems is considered. These systems are of a saddle point structure

 -     -  - -   - -
| A B^T || x | | c |
|       ||   |=|   |
| B -C  || y | | d |
 -     -  - -   - -
where A \in R^{n x n}, C \in R^{m x m} are symmetric and B \in R^{m x n}.

In Constraint preconditioning for indefinite linear systems, SIAM J. Matrix Anal. Appl., 21 (2000), Keller, Gould and Wathen introduced the idea of using constraint preconditioners that have a specific 2 by 2 block structure for the case of C being zero. We shall extend this idea by allowing the (2,2) block to be non-zero. Results concerning the spectrum and form of the eigenvectors are presented, as are numerical results to validate our conclusions.

Key words. preconditioning, indefinite linear systems, Krylov subspace methods



Optimization of a Dynamic Separation Problem Using MINLP Techniques


Stefan Emet (email: semet@abo.fi) and Tapio Westerlund
Abo Akademi University

In the present paper a dynamic separation problem is modeled and solved using Mixed Integer Nonlinear Programming (MINLP) techniques. The objective is to maximize the profit for continuous cyclic operation, and at the same time, to find the optimal configuration for the separation column system. The dynamics of the chromatographic separation process is modeled as a boundary value problem which is solved, within the optimization, using an iterative finite difference method. The separation of a fructose-glucose mixture is solved using the Extended Cutting Plane (ECP) method. It is shown that the production planning can be done efficiently for different purity requirements, such that all the output of a system can be utilized. Using a process design that is optimized it is thus possible to use existing complex systems, or to design new systems, more effciently and also to reduce the energy costs or the costs in general.



Multigrid Preconditioning for Anisotropic BTTB Systems


Rainer Fischer (email: rainer.fischer@mytum.de) and Thomas Huckle
Technical University of Munich

Multigrid methods are highly efficient solution techniques for large sparse multilevel Toeplitz systems which are positive semidefinite and ill-conditioned. In this paper, we develop multigrid methods which are especially designed for anisotropic two-level Toeplitz (BTTB) matrices. First, a method is described for systems with anisotropy along coordinate axes as a suitable combination of semicoarsening and full coarsening steps. Although the basic idea is known from the solution of PDEs, we present it here in a more formal way using generating functions and their level curves. This enables us not only to give a formal proof of convergence, but also to carry over the results to systems with anisotropy in other directions. After introducing new coordinates to describe these more complicated systems in terms of generating functions, we can solve them with the same efficiency. We hope that this will be useful in future applications such as the design of new multilevel preconditioners and solvers for the 2D-Helmholtz equation.



Dispatching buses in a depot using block patterns


Mohamed Hamdouni (email: mohamed.hamdouni@polymtl.ca),
Guy Desaulniers, Odile Marcotte, François Soumis, Marianne van Putten
Département de Mathématiques et génie Industriel, École polytechnique & GERAD

In this article we consider the problem of assigning parking slots to buses of different types so that the required buses can be dispatched easily in the morning. More precisely, if a bus of a certain type is needed at a given time, the buses that precede it in the lane must have departed already. Thus care must be taken to ensure that the buses arriving in the evening are parked in an order compatible with the types required for the morning departures. Maneuvers (i.e., rearrangements of buses within lanes) might be necessary to achieve this goal. Since the transit authorities need robust solutions to this problem (known as the dispatching problem in the literature), we formulate a model in which the depot lanes are filled according to specific patterns, called one-block or two-block patterns. We present two versions of this model, study their properties and show that some real-life instances can be solved within reasonable times by a commercial MIP solver.



Cache Efficient Data Structures and Algorithms for d-Dimensional Problems


Judith Hartman and Andreas Krahnke (email: akrahnke@hotmail.com)
TU München

Modern algorithms in numerical simulation need to combine efficient mathematical methods with concepts from computer science for nowadays computer architectures. In high-performance computing, memory access is a crucial factor and more important than pure cpu power. This memory boundedness can be reduced by utilizing modern hierarchical memory structures. However, many classical numerical methods lack the necessary data locality. In this paper we focus on the numerical solution of PDEs. First, we use a space-filling curve, the Peano curve, for the discretization and linearization of the domain. Then, we employ a system of stacks for processing the grid points linearly in a cache efficient way. In addition, this technique requires less memory for the organizational overhead, than for example a hierarchical approach using tree-structures.



A multigrid method for anisotropic PDE's in Elastic Image Registration


Lars Hoemke (email: hoemke@am.uni-duesseldorf.de)
Research Center Juelich

This paper deals with the solution of a nonlinear ill-conditioned inverse problem arising in digital image registration. In the first part of the paper, we define the problem as the minimization of a regularized non-linear least-squares functional, which measures the image difference and smoothness of the transformation. We study inexact Newton methods for solving this problem, i.e. we linearize the functional around a current approximation and replace the Hessian by a suitable operator in order to obtain well-posed subproblems in each step of the iteration.
These anisotropic subproblems are solved using a multigrid solver. Due to the anisotropy in the coefficients of the underlying equations, standard multigrid solvers suffer from poor convergence rates. We discuss modifications to the multigrid components, specifically to the smoothing procedure, the interpolation and the coarse grid correction. Numerical results that demonstrate the improvements obtained with these new components are given.



Parameter Tuning of Adaptive LQR-Repetitive controllers
Based on Genetic Algorithm:
Application to Uninterruptible Power Supply Systems


Mahdi Jalili-Kharaajoo (email: mahdijalili@ece.ut.ac.ir),
Mohammadreza Sadri and Farzad Habibipour Roudsari
Azad University and Iran Telecommunication Research Center, Tehran, Iran

Uninterruptible Power Supplies (UPS) systems, which supply emergency power in case of utility power failures, are used to interface critical loads such as computers and communication equipment to a utility power grid. Many discrete time controllers have been designed to control a single-phase inverter for use in UPS. In order to apply advanced control techniques in low coat micro controllers, a proper choice of the controller gains can be made off-line, which simplifies the control algorithm. Evolutionary algorithms are optimisation and search procedures inspired by genetics and the process of natural selection. This form of search evolves throughout generations improving the features of potential solutions by means of biologically inspired operations. On the ground of the structures undergoing optimisation the reproduction strategies, the genetic operators' adopted, evolutionary algorithms can be grouped in: evolutionary programming, evolution strategies, classifier systems, genetic algorithms and genetic programming. The genetic algorithms behave much like biological genetics. The genetic algorithms are an attractive class of computational models that mimic natural evaluation to solve problems in a wide variety of domains. A genetic algorithm comprises a set of individual elements (the population size) and a set of biologically inspired operators defined over the population itself etc. a genetic algorithms manipulate a population of potential solutions to an optimisation (or search) problem and use probabilistic transition rules. According to evolutionary theories, only the most suited elements in a population are likely to survive and generate offspring thus transmitting their biological heredity to new generations. In this paper, an adaptive LQR with repetitive controller (LQR-RP) for single-phase UPS application is designed. In the proposed controller, a RLS estimator identifies the plant parameters, which are used to compute the LQR gains periodically. The quadratic cost function parameters are chosen in order to reduce the energy of the control signal. The repetitive control action improves the response of the controller mainly in the presence of cyclic loads. It can be shown that using the repetitive controller with forgetting exponential coefficients, the control results can be improved significantly, which leads to an output voltage with low Total Harmonic Distortion (THD). There is not any efficient analytic mechanism to obtain the optimal values of the parameters of theses forgetting coefficients. So, we have to use intelligent optimization methods. Genetic Algorithms (GAs) is one the most efficient methods for this application. In this paper, using GAs the optimal tuning of repetitive controller's parameters is performed. Simulation results on the closed-loop system with obtained parameters conform the ability of the proposed method to achieve an output signal to low and admissible THD.



A family of root-finding algorithms converging faster to multiple roots


Yi Jin (email: yjin@paul.rutgers.edu)
Rutgers University

In this paper, we derive recursive and closed formulae for a family of root-finding algorithms that, contrary to typical algorithms, have higher order of convergence for multiple roots than for simple roots.

Key words: root-finding, order of convergence, symmetric functions, generating functions



Multilevel two-dimensional phase unwrapping


Iddit Shalem (email: shalemi@cs.technion.ac.il) and Irad Yavneh
Technion Israel Institute of Technology

Two-dimensional phase unwrapping is the problem of deducing unambiguous "phase" from values known only modulo 2 \pi. Many authors agree that the objective of phase unwrapping should be to find a weighted minimum of the number of places where adjacent (discrete) phase values differ by more than \pi. This NP-hard problem is of considerable practical interest, largely due to its importance in interpreting data acquired with synthetic aperture radar (SAR) inter-ferometry. Consequently, many heuristic algorithms have been proposed. Here we present an efficient multi-level graph algorithm for the approximate solution of an equivalent problem -- minimal residue cut in the dual graph.



Performance Study and Analysis of Parallel Multilevel Preconditioners


Chi Shen (email: cshen@crs.uky.edu) and Jun Zhang
Kentucky State University

The significant gap between peak and realized performance of parallel systems motivates the need for performance analysis. In order to predict the performance of a class of parallel multilevel ILU preconditioners (PBILUM), we build two performance prediction models for both the preconditioner construction phase and the solution phase. These models combine theoretical features of the preconditioners with estimates on computational cost, communications overhead, etc. Experimental simulations show that our model prediction based on certain reasonable assumptions is close to the simulation results. The models may be used to predict the performance of this class of parallel preconditioners.



Path-following and augmented Lagrangian methods for
contact problems in linear elasticity


Georg Stadler (email: ge.stadler@uni-graz.at)
University of Graz

A certain regularization technique for contact problems leads to a family of problems that can be solved efficiently using infinite-dimensional semismooth Newton methods, or in this case equivalently, primal-dual active set strategies. We present two procedures that use a sequence of regularized problems to obtain the solution of the original contact problem: first-order augmented Lagrangian, and path-following methods. The first strategy is based on a multiplier-update, while pathfollowing with respect to the regularization parameter uses theoretical results about the path value function to increase the regularization parameter appropriately. Comprehensive numerical tests investigate the performance of the proposed strategies for both a 2D as well as a 3D contact problem.

Back to Home page

Top