[1] "Solving the linear least squares problem on a linear array of processors" by Ahmed Sameh (Sept. 1982) [no abstract] [2] "A fast poisson solver for multiprocessors" by Ahmed Sameh (Jan. 1983) [no abstract] [3] "An overview of parallel algorithms in numerical linear algebra" by Ahmed Sameh (March 1983) [no abstract] [4] "The computation and communication complexity of a parallel banded system solver" by Ahmed Sameh (March 1983) [no abstract] Communicated by: Ahmed Sameh Obtainable from: Ahmed Sameh (Request in writing) ----- [5] Alan George and Esmond Ng,"Solution of sparse underdetermined systems of linear equations", Research report CS-82-39, Dept. of Computer Science, University of Waterloo, September 1982. (Presented in Sparse Matrix Symposium 1982.) In this paper we consider the problem of computing the minimal l2-solution to a consistent underdetermined linear system Ax=b, where A is m by n with m<=n. The method of solution is to reduce A to lower trapezoidal form [L O] using orthogonal transformations, where L is m by m and lower triangular. The method can be implemented efficiently if the matrix AA' is sparse. However if A contains some dense columns, AA' may be unacceptably dense. We present a method for handling these dense columns. The problem of solving a rank-deficient underdetermined system is also considered. [6] Esmond Ng,"Row elimination in sparse matrices using rotations", Research report CS-83-01, Dept. of Computer Science, University of Waterloo, January 1983. One way of solving a system of linear equations Ax=b, where A is an m by n matrix, is to use a QR-decomposition of A (or A' if m=n, and let Pr and Pc be permutation matrices of order m and n respectively. Suppose PrAPc is reduced to upper trapezoidal form [R' O]' using Givens rotations, where R is n by n and upper triangular. The sparsity structure of R depends only on Pc. For a fixed Pc, the number of arithmetic operations required to compute R depends on Pr. In this paper, we consider row ordering strategies which are appropriate when Pc is obtained from nested dissection orderings of A'. Recently, it was shown that so-called "width-2" nested dissection orderings of A'A could be used to simultaneously obtain good row and column orderings for A. In this paper, we show that the conventional (width-1) nested dissection orderings can also be used to induce good row orderings. In part I of this paper, we analyze the application of Givens rotations to a sparse matrix A using a bipartite graph model. [9] Alan George and Esmond Ng,"A partial pivoting implementation of Gaussian elimination for sparse systems", Research report CS-83-15, Dept. of Computer Science, University of Waterloo, May 1983. (Submitted to SIAM J. Sci. Stat. Comput.) In this paper, we consider the problem of solving a sparse nonsingular system of linear equations. We show that the structures of the triangular matrices obtained in the LU-decomposition of a sparse nonsingular matrix A using Gaussian elimination with partial pivoting are contained in those of the Cholesky factors of A', provided that the diagonal elements of A are nonzero. Based on this result, a method for solving sparse linear systems is then described. The main advantage of this method is that the numerical computation can be done using a static data structure. Numerical experiments comparing this method with other implementations of Gaussian elimination for solving sparse linear systems are presented and the results indicate that the method proposed in this paper is quite competitive with other approaches. Communicated by: Esmond Ng Obtainable from: Department of Computer Science University of Waterloo Waterloo, Ontario Canada N2L 3G1 ----- [10] Richard H. Bartels and Nezam Mahdavi-Amiri, "On Generating Test Problems for Nonlinear Programming Algorithms" We present an approach to test-problem generation which pro- vides a way of constructing example nonlinear programming prob- lems from a wide variety of given functions. The approach per- mits one to specify an arbitrary number of points at each one of which the problem should satisfy some "interesting" conditions (i.e. be optimal or stationary or degenerate) and to determine the characteristics of functions and derivatives at these points (i.e. choose predetermined values, gradients, Hessians and Lagrange multipliers). We give a sample of results obtained by using our approach to generate test problems for two algorithms to minimize nonlinear- least-squares objectives subject to nonlinear constraints. Communicated by: Richard Bartels Obtainable from: Department of Computer Science University of Waterloo Waterloo, Ontario ----- [11] V. S. Manoranjan, A. R. Mitchell, and J. Ll. Morris, "Numerical Solutions of the 'Good' Boussinesq Equation", Report NA/59, October, 1982. The 'good' Boussinesq equation is studied numerically using Galerkin methods and soliton solutions are shown to exist. An analytic formula for the 2-soliton interaction is presented and verified by numerical experiment. [12] D. F. Griffiths, A. R. Mitchell, "A Numerical Study of the Nonlinear Schroedinger Equation", Report NA/52, May, 1982. This paper presents a numerical study of the Nonlinear (cubic) Schroedinger (NLS) equation. Of particular interest are soliton solutions. The paper presents a one-parameter family of numerical schemes for solving the equation and advances a number of arguments, based on energy conservation, which isolates an optimum of the parameter. The numerical algorithm employs piecewise linear basis functions for space discretizations, and this is compared with a class of finite difference methods based upon a Crank-Nicolson scheme. Numerical results for the evolution of a single soliton and for the interaction of two solitons are provided. Communicated by: J. Ll. Morris, University of Waterloo. Obtainable from: D. F. Griffiths Department of Mathematics University of Dundee Dundee DD1 4HN Scotland, U.K. ----- [13] Thomas F. Coleman and Andrew R. Conn, "On the Local Convergence of a Quasi-Newton Method for the Nonlinear Programming Problem". In this paper we propose a new local quasi-Newton method to solve the equality constrained nonlinear programming prob- lem. The pivotal feature of the algorithm is that a projec- tion of the Hessian of the Lagrangian is approximated by a sequence of symmetric positive definite matrices. The matrix approximation is updated at every iteration by a projected version of the DFP or BFGS formula. We establish that the method is locally convergent and the sequence of x-values converges to the solution at a 2-step Q-superlinear rate. [14] Andrew R. Conn and Nicholas I.M. Gould, "ON THE LOCATION OF DIRECTIONS OF INFINITE DESCENT FOR NON-LINEAR PROGRAMMING ALGORITHMS". There is much current interest in general equality constrained programming problems, both for their own sake and for their applicability to active set methods for nonlinear programming. In the former case, typically, the issues are existence of solutions and their determination. In the latter instance, non-existence of solutions gives rise to directions of infinite descent. Such direction may subsequently be used to deter- mine a more desirable active set. The generalised Cholesky decomposition of relevant matrices enables us to answer the ques- tion of existence and to determine directions of infinite descent (when applicable) in an efficient and stable manner. The methods examined are related to implemen- tations that are suitable for null, range and Lagrangian methods. Communicated by: Andrew R. Conn Obtainable from: Andrew R. Conn Computer Science Department University of Waterloo Waterloo, Ontario CANADA N2L 3G1 ----- Title: "Bibliography of Stanford Computer Science Reports, 1962-1983", Author: "Kathryn A. Berg", Number: "STAN-CS-83-962", HCprice: "$5.00", MFprice: "$2.00", Note: "62 pages", Date: "March 1983" A chronological listing of all technical reports and theses published by the department. The availability and cost of each report is noted. Also included are the AI Memos and a number of the Heuristic Programming Project publications. Communicated by: Gene Golub Obtainable from: Computer Science Department Stanford University Stanford, California USA 94305 ----- [15] Performance of Various Computers Using Standard Linear Algebra Software in a Fortran Environment Jack J. Dongarra Mathematics and Computer Science Division Argonne National Laboratory Argonne, Illinois 60439 Tele: 312-972-7246 Arpa: Dongarra@anl-mcs In this note we compare a number of different computer systems for solving dense systems of linear equations using the LINPACK software in a Fortran environment. There are about 50 computers compared, ranging from a Cray X-MP to the 68000 based systems such as Apollo and Masscomp. [16] Squeezing the Most out of an Algorithm in Cray Fortran Jack J. Dongarra Mathematics and Computer Science Division Argonne National Laboratory and Stanley C. Eisenstat Department of Computer Science and Research Center for Scientific Computation Yale University Abstract - This paper discusses a technique for achieving super-vector performance on a Cray in a purely Fortran environment (i.e., without resorting to assembly language). The technique can be applied to a wide variety of algorithms in linear algebra, and is beneficial in other architectural settings. Communicated by: Jack Dongarra Obtainable from: Jack Dongarra Mathematics and Computer Science Division Argonne National Laboratory Argonne, Illinois 60439 ----- [17] "Block Preconditioning For The Conjugate Gradient Method" P. Concus, G.H. Golub and G. Meurant Block preconditionings for the conjugate gradient method are investigated for solving positive definite block tridiagonal systems of linear equations arising from discretization of boundary value problems for elliptic partial differentialequations. The preconditionings rest on the use of sparse approximate matrix inverses to generate incomplete block cholesky factorizations. Carrying out of the factorizations can be guaranteed under suitable conditions. Numerical experiments on test problems for two dimensions indicate that a particularly attractive preconditioning, which uses special properties of tridiagonal matrix inverses, can be computationally more efficient for the same computer storage than other preconditionings, including the popular point incomplete Cholesky factorization. Key words: conjugate gradient method, elliptic partial differential equations, incomplete factorization, iterative methods, preconditioning, sparse matrices. Communicated by: Gene Golub Computer Science Department Stanford University Stanford, California USA 94305 Available from: The authors ------- [18] Inaccuracy in Quasi-Newton Methods: Local Improvement Theorems (Rice University MASC Technical Report 83-11, March 1983.) by J.E. Dennis, Jr. Mathematical Sciences Department, Rice University, Houston, Texas 77251. and Homer F. Walker Department of Mathematics, University of Houston, Houston, Texas 77004. In this paper, we consider the use of bounded-deterioration quasi-Newton methods implemented in floating-point arithmetic to find solutions to F(x) = 0 where only inaccurate F-values are available. Our analysis is for the case where the relative error in F is less than one. We obtain theorems specifying local rates of improvement and limiting accuracies depending on the nearness to Newton's method of the basic algorithm, the accuracy of its implementation, the relative errors in the function values, the accuracy of the solutions of the linear systems for the Newton steps, and the unit-rounding errors in the addition of the Newton steps. Communicated by: John Dennis Available from: The authors --------------- [19] "Sharp Estimates for Multigrid Rates of Convergence" by Randolph E. Bank[$] and Craig C. Douglas[%] [$] Department of Mathematics, University of California at San Diego [%] Department of Computer Science, Yale University Yale University, Department of Computer Science Research Report #277 Available on request from: (U.S. Mail) (Phone) Craig C. Douglas 203/436-3761 Yale University Department of Computer Science (ARPAnet) 10 Hillhouse Ave. DOUGLAS-CRAIG@YALE P.O. Box 2158 Yale Station New Haven, CT 06520 Summary: In this paper, we prove the convergence of the multilevel iterative method for solving linear equations that arise from elliptic partial differential equations. Our theory is presented entirely in terms of the generalized condition number K of the matrix A and the smoothing matrix B. This leads to a completely algebraic analysis of the method as an iterative technique for solving linear equations; the properties of the elliptic equation and discretization procedure enter only when we seek to estimate K, just as in the case of most standard iterative methods. Here we consider the fundamental two level iteration, and the V and W cycles of the j-level iteration (j > 2). We prove that the V and W cyles converge even when only one smoothing iteration is used. We present several examples of the computation of K using both Fourier analysis and standard finite element techniques. We compare the predictions of our theorems with the actual rate of convergence. Our analysis also shows that adaptive iterative methods, both fixed (Chebyshev) and adaptive (conjugate gradients and conjugate residuals), are effective as smoothing procedures. ------- ================================================== NUMERICAL ANALYSIS TITLES -- Volume 1, Number 2 September, 1983 ================================================== [1] J.J. Dongarra and C.B. Moler, EISPACK - A Package for Solving Matrix Eigenvalue Problems, Argonne National Laboratory, Mathematics and Computer Science Division, ANL/MCS-TM-12, August 1983. This is a two part paper on the software package EISPACK. The first part describes the development and usage of EISPACK and is intended as an introduction to the package. The second part describes changes made to the package for the current release and how to install the package. EISPACK is a systemized collection of subroutines that compute the eigenvalues and eigenvectors of different classes of matrices. Submitted by: Jack Dongarra Obtainable from: Jack Dongarra Mathematics and Computer Science Division Argonne National Laboratory Argonne, Illinois 60439 --------------- [2] ESTIMATION OF SPARSE HESSIAN MATRICES AND GRAPH COLORING PROBLEMS Thomas F. Coleman and Jorge J. More' Large scale optimization problems often require an approximation to the Hessian matrix. If the Hessian matrix is sparse then estimation by differences of gradients is attractive because the number of required differences is usually small compared to the dimension of the problem. The problem of estimating Hessian matrices by differences can be phrased as follows: Given the sparsity structure of a symmetric matrix A, obtain p vectors d(1), ..., d(p) such that Ad(1), ..., Ad(p) determine A uniquely with p as small as possible. We approach this problem from a graph theoretic point of view and show that both direct and indirect approaches to this problem have a natural graph coloring interpretation. The complexity of the problem is analyzed and efficient practical heuristic procedures are developed. Numerical results illustrate the differences between the various approaches. Submitted by: Jorge J. More' Obtainable from: Jorge More' Mathematics and Computer Science Division Argonne National Laboratory Argonne, Illinois 60439 --------------- [3] SOFTWARE FOR ESTIMATION OF SPARSE JACOBIAN MATRICES Thomas F. Coleman, Burton S. Garbow, and Jorge J. More' In many nonlinear problems it is necessary to estimate the Jacobian matrix of a nonlinear mapping F. In large scale problems the Jacobian of F is usually sparse, and then estimation by differences is attractive because the number of differences can be small compared to the dimension of the problem. For example, if the Jacobian matrix is banded then the number of differences needed to estimate the Jacobian matrix is, at most, the width of the band. In this paper we describe a set of Fortran 77 subroutines whose purpose is to estimate the Jacobian matrix of a mapping F with the least possible number of function evaluations. Submitted by: Jorge J. More' Obtainable from: Jorge More' (address above) --------------- [4] THE MINPACK PROJECT Jorge J. More', Danny C. Sorensen, Burton S. Garbow, and Kenneth E. Hillstrom The MINPACK project is a research effort whose goal is the development of a systematized collection of quality optimization software. In this paper we provide an overview of the MINPACK project. We describe the circumstances that led to the project and the various changes that occurred during the next four years. We also describe the techniques we have used to produce reliable, easy-to-use, transportable software. We discuss those topics that are relevant to the development of software optimization libraries: scale invariance, robustness, convergence testing, user interface and documentation standards. Submitted by: Jorge J. More' Obtainable from: Jorge J. More' (address above) --------------- [5] A Nonlinear Version of Gauss's Minimum Variance Theorem with Applications to an Errors-in-the-Variables Model G. W. Stewart TR-1263, April 1983 This paper consists of two parts. In the first, a theorem of Gauss, commonly known as the Gauss- Markov theorem, is generalized and extended to apply to nonlinear least squares estimation. The result justifies the use of nonlinear least squares under the classical assumptions that the errors have mean zero and equal variances, with the additional proviso that the variance is suffi- ciently small. In the second part of the paper the theory developed in the first is applied to an errors-in-the-variables model to justify the use of a form of latent root regression to estimate regression coefficients. However, it is further shown that in terms of the theory ordinary least squares is equally good, which suggests that the latter is to be prefered on grounds of computa- tional simplicity. Submitted by: G. W. Stewart Obtainable from: G. W. Stewart University of Maryland Department of Computer Science College Park, MD 20742 U.S.A. --------------- [6] On the Asymptotic Behavior of Scaled Singular Value and QR Decompositions. G. W. Stewart TR-1307, July 1983 Asymptotic expressions are derived for the singu- lar value decomposition of a matrix some of whose columns approach zero. Expressions are also derived for the QR factorization of a matrix some of whose rows approach zero. The expressions give insight into the method of weights for approximat- ing the solutions of constrained least squares problems. Submitted by: G. W. Stewart Obtainable from: G. W. Stewart (address above) --------------- [7] A Jacobi-like Algorithm for Computing the Schur Decomposition of a Non-Hermitian Matrix G. W. Stewart TR-1321, August 1983 (Dedicated to Peter Henrici on his sixtieth birthday.) This paper describes an iterative method for reducing a general matrix to upper triangular form by unitary similarity transformations. The method is similar to Jacobi's method for the symmetric eigenvalue problem in that it uses plane rotations to annihilate off-diagonal elements, and when the matrix is Hermitian it reduces to a variant of Jacobi's method. Although the method cannot com- pete with the QR algorithm in serial implementa- tion, it admits of a parallel implementation in which a double sweep of the matrix can be done in time proportional to the order of the matrix. Submitted by: G. W. Stewart Obtainable from: G. W. Stewart (address above) --------------- [8] On Lipschitzian Proximity Maps J. R. Respess and E. W. Cheney CNA-174, July 1981 [No abstract] Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis R. L. Moore Hall, 13.150 University of Texas Austin, Texas 78712 (At the time the report is sent, an invoice in the amount of $1 per copy to cover duplication and postage will be included.) --------------- [9] The Characterization of Best Approximations in Tensor Product Spaces W. A. Light and E. W. Cheney CNA-175, September 1981 Given compact Hausdorff spaces S and T, a subspace W of C(SxT) is defined by W = G*C(T) + C(S)*H where G and H are finite dimensional subspaces of C(S) and C(T) respectively, and "*" represents the tensor product. Necessary and sufficient conditions are developed for a function in C(SxT) to have a best approximation in W. Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis (address above, $1 charge) --------------- [10] Generalized Conjugate Gradient Acceleration of Iterative Methods Kang Chang Jea CNA-176, February 1982 Acceleration of basic iterative methods is considered for solving large sparse systems of linear equations which are ei- ther not symmetrizable or for which no symmetrization matrix is conveniently available. Work on developing a method which will work as well as Chebyshev acceleration but requires no eigenvalue estimates is described. Algorithms are analyzed and numerical re- sults given. Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis (address above, $1 charge) --------------- [11] Adapting ITPACK Routines for Use on a Vector Computer David R. Kincaid, Thomas Oppe, and David M. Young CNA-177, August 1982 A description of recent work involving vectorization of the ITPACK package of iterative algorithms for the solution of large sparse systems of linear equations is given along with nu- merical results from runs on the CDC CYBER 203 and 205 computers. A more complete vectorization involving, in some cases, total re- design of algorithms is projected. Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis (address above, $1 charge) --------------- [12] ITPACK on Supercomputers David R. Kincaid and Thomas Oppe CNA-178, September 1982 Results of preliminary testing of a vector version of ITPACK are presented. Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis (address above, $1 charge) --------------- [13] The Best Approximation of Bivariate Functions by Separable Functions M. von Golitschek and E. W. Cheney CNA-179, December 1982 [No abstract] Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis (address above, $1 charge) --------------- [14] The ITPACK Project: Past, Present, and Future David R. Kincaid and David M. Young CNA-180, March 1983 The objectives of the ITPACK project are reviewed, pro- gress to date is summarized, and plans for future work are out- lined. The present packages for use on scalar and vector compu- ters are described. An overview is given of developmental work on other software involving preconditioners and nonsymmetric proce- dures. Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis (address above, $1 charge) --------------- [15] Accelerating Nonsymmetrizable Iterative Methods David M. Young, Kang Chang Jea, and David R. Kincaid CNA-181, March 1983 We discuss recent work done at the University of Texas and elsewhere on accelerating basic iterative methods for the solution of systems of linear equations which have a nonsym- metrizable coefficient matrix. Approaches discussed include that of Concus and Golub and of Widlund as well as the use of normal equations, Chebyshev acceleration, and generalized conjugate gra- dient procedures. Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis (address above, $1 charge) --------------- [16] The ITPACK Software Package David M. Young and David R. Kincaid CNA-182, June 1983 We discuss iterative methods for the solution of sys- tems of linear equations in general and their implementation in ITPACK. Future modifications of the package including vectoriza- tion and modularization are considered. Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis (address above, $1 charge) --------------- [17] The Embedding of Proximinal Sets Carlo Franchetti and E. W. Cheney CNA-183, June 1983 Given Banach spaces Y < X < Z (where "<" denotes the subset relation) such that Y is proximinal in X, we discuss con- ditions under which Y remains proximinal when considered as a subspace of Z. Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis (address above, $1 charge) --------------- [18] Minimal Projections in Tensor Product Spaces Carlo Franchetti and E. W. Cheney CNA-184, July 1983 [No abstract] Submitted by: David Kincaid Obtainable from: Center for Numerical Analysis (address above, $1 charge) --------------- [19] Splines in Interactive Computer Graphics Richard H. Bartels (Invited presentation, Dundee Conference, 1983) Computer graphics, particularly interactive computer graphics, is not, as the name might imply, concerned with drawing graphs, but rather with the broadest issues of manipulating, transforming, and displaying information in visual format. It is interactive in so far as operations can be carried out in real time -- which requires algorithms of high computational efficiency and low com- plexity. Splines are a valuable tool in graphics, but they are often applied in a way not used by the mathematician. This difference raises computational issues which the numerical analyst might otherwise never see. This talk will provide a brief introduction to such issues and follow with a study of two current developments. We begin with a review of the graphics en- vironment, mentioning the modelling and display process and pointing out some of the costly issues. The novel use of splines in interactive graphics comes through the construction of sur- faces as weighted averages of selected points, called control vertices in which B-splines are taken as the weighting functions. Some examples will illustrate the characteristics of this use of B-splines. With this background we consider two recent develop- ments. The first is the control-vertex recurrence of Riesenfeld, Cohen, and Lyche; the second is Barsky's work on geometric vs. mathematical continuity, and his introduction of Beta-splines. We will close with some results on current research concerned with a synthesis of these two developments. Submitted by: Richard Bartels Obtainable from: Richard Bartels University of Waterloo Department of Computer Science Waterloo, Ontario Canada N2L 3G1 --------------- [20] NUMERICAL DECOMPOSITION OF A CONVEX FUNCTION Mike Lamoureux Henry Wolkowicz Department of Mathematics Department of Mathematics Stanford University University of Alberta Given the n by p orthogonal matrix A and the convex function f: R(n) to R, we find two orthogonal matrices P and Q such that f is 'almost' constant on the convex hull of the columns of P and -P, f is 'sufficiently' nonconstant on the column space of Q, and the column spaces of P and Q provide an orthogonal direct sum decomposition of the column space of A. This provides a numerically stable algorithm for calculating the cone of directions of constancy, at a point x, of a convex function. Applications to convex programming are discussed. Submitted by: Henry Wolkowicz Obtainable from: Henry Wolkowicz Department of Mathematics University of Alberta Edmonton, Canada T6G 2G1 --------------- ================================================== NUMERICAL ANALYSIS TITLES -- Volume 2, Number 1 March, 1984 ================================================== [1] Daniel Szyld: A Two-level Iterative Method for Large Sparse Generalized Eigenvalue Calculations We present a method for the solution of the generalized symmetric eigenvalue problem. The matrices in the pencil are assumed to be too expensive to factor. The outer level of iteration is a step of either Inverse or Rayleigh Quotient iteration. The choice is made at each step to guarantee superlinear convergence to an eigenvalue in (or close to) a prescribed interval. The indefinite linear system is solved at each step with the conjugate gradient type method SYMMLQ, for which a variational formulation and error bounds are presented. Submitted by: Daniel Szyld Obtainable from: Daniel B. Szyld Institute for Economic Analysis New York University 269 Mercer Street, Second Floor New York, NY 10003 --------------- [2] Daniel Szyld: Preliminary Results in Implementing a Model of the World Economy on the Cyber 205: A Case of Large Sparse Nonsymmetric Linear Equations presented at the Cyber 200 Applications Seminar, Oct.10-12, 1983, NASA Godard Ctr. Md. A brief description of the Institutes' Model of the Wrold Economy is presented, together with our experience in converting the software to vector code.Includes running times. Submitted by: Daniel Szyld Obtainable from: Daniel B. Szyld (address above) --------------- [3] G. W. Stewart A Generalization of the Eckart-Young Approximation Theorem TR-1325 Department of Computer Science University of Maryland at College Park September 1983 The Eckart-Young theorem solves the problem of approximating a matrix by one of lower rank. How- ever, the approximation generally differs from the original in all its elements. In this paper it is shown how to obtain a best approximation of lower rank in which a specified set of columns of the matrix remains fixed. The paper concludes with some applications of the generalization. Submitted by: G. W. Stewart Obtainable from: G. W. Stewart Department of Computer Science and Institute for Physical Science and Technology The University of Maryland at College Park College Park, Maryland 20742 --------------- [4] Dianne P. O'Leary and G. W. Stewart DATA-FLOW ALGORITHMS FOR PARALLEL MATRIX COMPUTATIONS In this work we develop some algorithms and tools for solv- ing matrix problems on parallel processing computers. Operations are synchronized through data-flow alone, which makes global synchronization unnecessary and enables the algorithms to be implemented on machines with very simple operating systems and communications protocols. As exam- ples, we present algorithms that form the main modules for solving Liaponuv matrix equations. We compare this approach to wavefront array processors and systolic arrays, and note its advantages in handling missized problems, in evaluating variations of algorithms or architectures, in moving algo- rithms from system to system, and in debugging parallel algorithms on sequential machines. Submitted by: G. W. Stewart Obtainable from: G. W. Stewart (address above) --------------- [5] GENE H. GOLUB AND DAVID MAYERS THE USE OF PRE-CONDITIONING OVER IRREGULAR REGIONS Some ideas and techniques for solving elliptic pde.'s over irregular regions is discussed. The basic idea is to break up the domain into subdomains and then to use the pre-conditioned conjugate gradient method for obtaining the solution over the entire domain. The solution of Poisson's equation over a T-shaped region is described in some detail and a numerical example is given. Submitted by: Gene H. Golub Obtainable from: Gene H. Golub Computer Science Department Stanford University Stanford, California 94305 --------------- [6] Robert Schreiber Computing Genrealized Inverses and Eigenvalues of Symmetric Matrices Using Systolic Arrays. Manuscript NA-83-03. Stanford NA Project New systolic arrays methods are presented for tridiagonalization of a dense, symmetric n x n matrix in O(n**3/2) time, tridiagonalization of a matrix of bandwidth b in O(bn) time, and finding the eigenvalues of a symmetric tridiagonal matrix in O(n) time. Systolic arrays for implementing an iterative method for the generalized inverse are given. They allow computation of the generalized inverse of an m x n matrix in O(m) time. Submitted by: Gene Golub Obtainable from: Gene Golub (address above) --------------- [7] Rodrigo Fontecilla Technical Report 1370 Department of Computer Science University of Maryland at College Park January 1984 In solving nonlinear equality constrained optimization prob- lems using secant methods we often have to solve two or three linear systems of equations at each iteration. We use the theory developed by Dembo, Eisenstat and Steihaug in inexact Newton methods to solve at least one of the linear system inexactly. We analyze the case when the number of constraints is much smaller than the number of variables. In this case we solve the smallest linear system (the one used to find the multipliers) exactly while solving the bigger system (the one giving the step on the variables) inexactly at each iteration. Algorithms using the DFP and the BFGS secant updates are proposed. Conditions on the residual are given in order to allow the user the possibility of a trade- off between the rate of convergence (q-superlinear) and the amount of work per iteration. Submitted by: Rodrigo Fontecilla Obtainable from: Rodrigo Fontecilla <rod%umcp-cs@csnet-cic> Department of Computer Science and Institute for Physical Science and Technology University of Maryland College Park, Maryland 20742 --------------- [8] Richard H. Bartels and John J. Jezioranski On Least Squares Fitting Using Orthogonal Multinomials Forsythe (1) has given a method for generating basis polynomials in a single variable which are orthogonal with respect to a given inner product. Weisfeld (2) later demonstrated that Forsythe's approach could be extended to polynomials in an arbitrary number of variables. In this paper we sharpen Weisfeld's results and present a program for computing weighted, multinomial, least squares approximations to discrete data. ----- (1) G. E. Forsythe (1957), Generation and Use of Orthogonal Polynomials for Data-Fitting with a Digital Computer, J. SIAM 5, 74-87. (2) M. Weisfeld (1959), Orthogonal polynomials in Several Variables, Numer. Math. 1, 38-40. Submitted by: Richard H. Bartels Obtainable from: Richard H. Bartels University of Waterloo Department of Computer Science Waterloo, Ontario Canada N2L 3G1 --------------- [9] Ulrich Tulowitzki Inexact Successive Linearization Methods for Variational Inequality Problems We extend the inexact Newton rate of convergence characterization of Dembo, Eisenstat and Steihaug to variational inequality problems. Our characterization theorem is expressed in terms of the relative error in the solution of a linearized subproblem and does not require the set of active inequalities to remain constant. An immediate consequence of our theorem is a practical termination criterion for truncating an iterative procedure applied to the subproblems in successive quadratic programming algorithms for nonlinear optimization. Submitted by: Ulrich Tulowitzki Obtainable from: Ulrich Tulowitzki <Tulowitzki@YALE.ARPA> Yale University School of Organization and Management Box 1-A New Haven, CT 06520 --------------- [10] Ron S. Dembo and Siddhartha Sahi A Convergent Framework for Constrained Nonlinear Optimization Working Paper, Series B #69 School of Organization and Management Yale University September, 1983 We describe a simple, practical algorithmic framework for constrained nonlinear optimization. Any algorithm that may be expressed in this framework, and most existing algorithms can, will converge to a local minimum from an arbitrary feasible starting point. The framework is particularly suitable for the analysis and development of algorithms for large-scale optimization, since it permits radical changes in the active set to occur at each step and can be implemented in terms of quantities that are easily computed. Submitted by: Ron Dembo Obtainable from: Ron Dembo Yale University School of Organization and Management Box 1-A New Haven, CT 06520 --------------- ================================================== NUMERICAL ANALYSIS TITLES -- Volume 2, Number 2 June, 1984 ================================================== [1] Ron S. Dembo and Ulrich Tulowitzki On the Minimization of Quadratic Functions Subject to Box Constraints SOM Working Paper Series B #71 Yale University School of Organization and Management 1983 We study efficient algorithms for very large quadratic problems subject to box constraints an propose some new convergent methods that permit the addition and deletion of many constraints each time a search direction is calculated. Numerical experiments, on problems of up to 10,000 variables, indicate that the perfor- mance of some of the proposed algorithms appears to be insensitive to the number of binding constraints at the optimal solution or to the starting point. Also it appears as if solving the box constrained problem is at worst marginally more expensive than solving the same problem without constraints. Submitted by: Ron Dembo Obtainable from: Ron Dembo Yale University School of Organization and Management P.O. Box 1A New Haven, Conn. 06520 --------------- [2] Ron S. Dembo A Primal Truncated Newton Algorithm with Application to Large-Scale Nonlinear Network Optimization SOM Working Paper Series B #72 Yale University School of Organization and Management March 1983 We describe a new, convergent, primal-feasible algorithm for linearly constrained optimization. It is capable of rapid asymptotic behavior and has relatively low storage requirements. Its application to large scale nonlinear network optimization is discussed and computational results on problems of over 2,000 variables and 1,000 constraints are presented. Indications are that it could prove to be significantly better than known methods for this class of problem. Submitted by: Ron Dembo Obtainable from: Ron Dembo (address above) --------------- ================================================= NUMERICAL ANALYSIS TITLES -- Volume 2, Number 3 September, 1984 ================================================= [1] R. C. Y. Chin, G. W. Hedstrom, and F. A. Howes Considerations on Solving Problems with Multiple Scales Lawrence Livermore National Laboratory Report UCRL-90791 (No abstract given.) Submitted by: Gerald Hedstrom, L-387 Obtainable from: Gerald Hedstrom, L-387 P. O. Box 808 LLNL Livermore, California 94550 (na.hedstrom@su-score.arpa) --------------- [2] R. C. Y. Chin, G. W. Hedstrom, and Lewis Thigpen Numerical methods in seismology J. Comput. Phys. 54 (1984), 18-56. (No abstract given.) Submitted by: Gerald Hedstrom, L-387 Obtainable from: Gerald Hedstrom, L-387 P. O. Box 808 LLNL Livermore, California 94550 (na.hedstrom@su-score.arpa) --------------- [3] R. C. Y. Chin, G. W. Hedstrom, and Lewis Thigpen Matrix methods in synthetic seismograms Geophys. J. R. astron. Soc. 77 (1984), 483-502. Submitted by: Gerald Hedstrom, L-387 Obtainable from: Gerald Hedstrom, L-387 P. O. Box 808 LLNL Livermore, California 94550 (na.hedstrom@su-score.arpa) --------------- [4] W. H. Enright and J. D. Pryce Two Fortran Packages for Assessing Initial Value Methods University of Toronto Computer Science Department TR 167/83 (No abstract given.) Submitted by: Kenneth R. Jackson Obtainable from: USENET {decvax linus ihnp4 allegra floyd decwrl garfield qucis cornell mcgill-vision sask watmath uw-beaver ubc-vision}!utcsrgv!enright or CSNET enright@toronto --------------- [5] Allan D. Jepson and Alastair Spence The Numerical Solution of Nonlinear Equations having Several Parameters University of Toronto Computer Science Department TR 168/83 (No abstract submitted.) Submitted by: Kenneth R. Jackson Obtainable from: USENET { decvax linus ihnp4 allegra floyd decwrl garfield qucis cornell mcgill-vision sask watmath uw-beaver ubc-vision }!utcsrgv!jepson or CSNET jepson@toronto --------------- [6] Tony F. Chan and Kenneth R. Jackson The use of Iterative Linear-Equation Solvers in Codes for Large Systems of Stiff IVPs for ODEs University of Toronto Computer Science Department TR 170/84. (No abstract given.) Submitted by: Kenneth R. Jackson Obtainable from: USENET { decvax linus ihnp4 allegra floyd decwrl garfield qucis cornell mcgill-vision sask watmath uw-beaver ubc-vision }!utcsrgv!krj CSNET krj@toronto --------------- [7] Robert Schreiber Computing Generalized Inverses and Eigenvalues of Symmetric Matrices Using Systolic Arrays. Stanford NA Manuscript NA-83-03. New systolic array methods are presented for tridiagonalization of a dense, symmetric n x n matrix in O(n**3/2) time, tridiagonalization of a matrix of bandwidth b in O(bn) time, and finding the eigenvalues of a symmetric, tridiagonal matrix in O(n) time. Systolic arrays for implementing an iterative method for the generalized inverse are given. They allow computation of the generalized inverse of an m x n matrix in O(m) time. Submitted by: Robert Schreiber Obtainable from: Robert Schreiber ESL Inc. 495 Java Drive, Sunnyvale, CA 94086. --------------- [8] Robert Schreiber Implementation of Adaptive Array Algorithms Stanford NA Manuscript NA-84-28. The computational procedures used to implement some standard and some recently proposed adaptive methods for direction finding and beamforming by sensor arrays are considered. Among these methods are the minimum variance beamforming method and several high-resolution methods that are based on the covariance matrix of the signal. We discuss recursive implementation of these procedures: the covariance matrix is periodically updated by the addition of a matrix of rank one. We propose some efficient, numerically stable algorithms for updating the solution. Submitted by: Robert Schreiber Obtainable from: Robert Schreiber ESL Inc. 495 Java Drive, Sunnyvale, CA 94086. --------------- [9] S. S. Chen, J. J. Dongarra, and C. C. Hsiung Multiprocessing Linear Algebra Algorithms on the CRAY X-MP-2: Experiences with Small Granularity Argonne National Laboratory Report MCS-TM-24 (February 1984) This paper gives an overview of the CRAY X-MP-2 general-purpose multiprocessor system and discusses how it can be used effective- ly to solve problems that have small granularity. An implementa- tion is described for linear algebra algorithms that solve sys- tems of linear equations when the matrix is general and when the matrix is symmetric and positive definite. Submitted by: Jack J. Dongarra Obtainable from: Jack J. Dongarra Mathematics and Computer Science Division Argonne National Laboratory Argonne, Illinois 60439 --------------- [10] J. J. Dongarra Performance of Various Computers Using Standard Linear Equations Software in a Fortran Environment Argonne National Laboratory Report MCS-TM-23 (updated August 1984) This note compares the performance of different computer systems while solving dense systems of linear equations using the LINPACK software in a Fortran environment. About 50 computers, ranging from a CRAY X-MP to the 68000-based systems such as the Apollo and SUN workstations, are compared. Submitted by: Jack J. Dongarra Obtainable from: Jack J. Dongarra Mathematics and Computer Science Division Argonne National Laboratory Argonne, Illinois 60439 --------------- [11] J. J. Dongarra and A. Sameh On Some Parallel Banded System Solvers Argonne National Laboratory Report MCS-TM-27 (March 1984) This paper describes algorithms for solving narrow banded systems and the Helmholtz difference equations that are suitable for mul- tiprocessing systems. The organization of the algorithms highlight the large grain parallelism inherent in the problems. Submitted by: Jack J. Dongarra Obtainable from: Jack J. Dongarra Mathematics and Computer Science Division Argonne National Laboratory Argonne, Illinois 60439 --------------- [12] J.J. Dongarra, E.L. Lusk, R.A. Overbeek, B.T. Smith, and D.C. Sorensen New Directions in Software for Advanced Computer Architectures Argonne National Laboratory Report MCS/TM 32 (August, 1984) This report contains a collection of short position papers that were prepared for the Workshop on Taxonomy of Parallel Algorithms sponsored by Los Alamos National Laboratory and held at Santa Fe, New Mexico Nov. 30 - Dec. 2, 1983. These papers represent the authors views on directions of research that should be undertaken with regards to programming methodology, algorithm design, and software implementation for advanced computer architectures. Several projects have been initiated by the authors within the general framework outlined in these articles since they were written. The concepts formulated here are expected to serve as a foundation for future work in multiprocessing. Submitted by: Jack J. Dongarra Obtainable from: Jack J. Dongarra Mathematics and Computer Science Division Argonne National Laboratory Argonne, Illinois 60439 --------------- [13] D. C. Sorensen Analysis of Pairwise Pivoting in Gaussian Elimination Argonne National Laboratory Report MCS-TM-26 (February 1984) The method of Gaussian Elimination using triangularization by elementary stabilized matrices constructed by pairwise pivoting is analyzed. It is shown that a variant of this scheme which is suitable for implementation on a parallel computer is numerically stable although the bound is larger than the one for the standard partial pivoting algorithm. Submitted by: Jack J. Dongarra Obtainable from: Jack J. Dongarra Mathematics and Computer Science Division Argonne National Laboratory Argonne, Illinois 60439 --------------- [14] D. C. Sorensen Buffering for Vector Performance on a Pipelined MIMD Machine Argonne National Laboratory Report MCS-TM-29 (April 1984) A technique is presented for obtaining vector performance from a pipelined MIMD computer that does not have hardwired vector in- structions. The specific computer in mind is the Denelcor HEP, but the technique might influence the use and possibly even the design of future machines with this type of architecture. This preliminary report presents the basic idea and demonstrates that it can be implemented. Buffering blocks of data to registers is used in conjunction with pipelined floating point operations to achieve vector performance. Empirical evidence is presented to show that up to 5.8 megaflop performance is possible from the Denelcor HEP on very regular tasks such as matrix vector pro- ducts. While this rate is not in the "super-computer" range, it is certainly respectable given the hardware capabilities of the HEP (this machine is rated at 10 MIPS peak). This performance indicates that an apparently minor refinement to the architectur- al design could provide very efficient vector operations in addi- tion to the parallelism and low-overhead synchronization already offered by the HEP architecture. Submitted by: Jack J. Dongarra Obtainable from: Jack J. Dongarra Mathematics and Computer Science Division Argonne National Laboratory Argonne, Illinois 60439 --------------- [15] Thomas F.Coleman and Alex Pothen The Sparse Null Space Problem Tecnical Report:TR 84-598 Computer Science Cornell University The sparse null space basis problem is the following. Let A be a given t*n matrix, t < n, of rank t. Find a matrix N, with the fewest nonzeroes in it, whose columns span the null space of A. This problem arises in the design of practical algorithms for large scale numerical optimization problems. In this paper we analyze the complexity of this problem, using combinatorial tools, and we propose approximation algorithms. Submitted by: Thomas F. Coleman Obtainable from: Thomas F.Coleman Computer Science Department Cornell University Ithace NY 14853 (coleman@cornell.csnet) --------------- [16] Richard H. Byrd An Example of Irregular Convergence in Some Constrained Optimization Methods That Use the Projected Hessian In this paper we give examples illustrating the behavior of the Coleman- Conn horizontal vertical method and of successive quadratic programming with a Hessian approximation exact on the tangent space of the constraints. One example shows that these methods in general are not one-step superlinearly convergent. Submitted by: Richard H. Byrd Obtainable from: Richard H. Byrd Department of Computer Science University of Colorado Boulder, Colorado 80309 --------------- [17] Richard H. Byrd On the Convergence of Constrained Optimization Methods with Accurate Hessian Information on a Subspace In this paper we analyze methods of the type proposed by Coleman and Conn for nonlinearly constrained optimization. We show that if the reduced Hessian approximation is sufficiently accurate, then the method generates a sequence of iterates which converges one-step superlinearly. This result applies to a quasi-Newton implementation. If the reduced exact Hessian is used the method has an R-order equal to that of the secant method. We also prove a similar result for a modified version of successive quadratic programming. Finally some parallels between convergence results for methods that approximate the reduced Hessian, and multiplier methods that use the reduced Hessian inverse are pointed out. Submitted by: Richard H. Byrd Obtainable from: Richard H. Byrd Department of Computer Science University of Colorado Boulder, Colorado 80309 --------------- [18] G. H. Golub, Pierre Manneback, and Phillippe Toint A COMPARISON BETWEEN SOME DIRECT AND ITERATIVE METHODS FOR LARGE SCALE GEODETIC LEAST SQUARES PROBLEMS The purpose of this paper is to describe and compare some numerical methods for solving large dimensional linear least squares problems that arise in geodesy and, more specially, from Doppler posi- tioning. The methods that are considered are the direct orthogonal decompo- sition, and the combination of conjugate gradient type algorithms with projects as well as the exploitatioin of "Property A". Numerical results are given and the respective advantage of thef methods are discussed with respect to such parameters as CPLU time, input/output and storage require- ments. Extensions of the results to more general problems are also discussed. Submitted by: G. H. Golub Obtainable from: G. H. Golub Stanford University Computer Science Department Stanford, CA 94305 (golub@navajo.arpa) --------------- [19] Franklin T. Luk A Jacobi-like Method for Computing the QR-decomposition Technical Report 84-612 Department of Computer Science Cornell University A parallel Jacobi-like method for computing the QR-decomposition of an nxn matrix is proposed. It requires O(n**2) processors and O(n) units of time. The method can be extended to handle an mxn matrix (m>=n). The requirements become O(n**2) processors and O(m) time. Submitted by: Franklin T. Luk Obtainable from: Librarian Department of Computer Science Cornell University Ithaca, NY 14853 --------------- [20] Franklin T. Luk A Triangular Processor Array for Computing the Singular Value Decomposition. Technical Report 84-625 Department of Computer Science Cornell University A triangular processor array for computing a singular value decomposition (SVD) of an mxn (m>=n) matrix is proposed. A Jacobi-type method is used to first triangularize the given matrix and then diagonalize the resultant triangular form. The requirements are O(m) time and n**2/4 + O(n) processors. Submitted by: Franklin T. Luk Obtainable from: Librarian Department of Computer Science Cornell University Ithaca, NY 14853 --------------- [21] Richard P. Brent and Franklin T. Luk The Solution of Singular Value Problems Using Systolic Arrays. Technical Report 84-626 Department of Computer Science Cornell University This paper concerns the computation of the singular value decomposition using systolic arrays. Two different linear time algorithms are presented. Submitted by: Franklin T. Luk Obtainable from: Librarian Department of Computer Science Cornell University Ithaca, NY 14853 --------------- ================================================== NUMERICAL ANALYSIS TITLES -- Volume 2, Number 4 November 23, 1984 ================================================== [1] Finbarr O'Sullivan and Grace Wahba A Cross Validated Bayesian Retrieval Algorithm For Non-Linear Remote Sensing Experiments To appear, Journal of Computational Physics We present a fully automated retrieval algorithm for estimating the solution of a mildly non linear Fredholm Integral Equation of the First Kind, when the data are discrete, noisy (experimental) measurements. This algorithm combines a simple Gauss-Newton interation with an extended form of Generalized Cross Validation , applied in the context of Tihonov Regularization with moment discretization. The performance of the algorithm is illustrated in the context of a remote sensing problem arising in satellite meteorology. Submitted by: Grace Wahba Obtainable from: Grace Wahba Department of Statistics University of Wisconsin, Madison (na.wahba@su-score.arpa) --------------- [2] ANDREW R. CONN NONLINEAR PROGRAMMING, EXACT PENALTY FUNCTIONS AND PROJECTION TECHNIQUES FOR NON-SMOOTH FUNCTIONS We present a personal overview of various approaches to solving nonlinear programs with nonlinear constraints that make use of the L-1 exact penalty function. The advantages, disadvantages and related remaining difficulties of these approaches will be considered. Finally some recent research and extensions are given. Submitted by: Andrew R. Conn Obtainable from: Andrew R. Conn Computer Science Department University of Waterloo Waterloo, Ontario N2L 3G1 CANADA (arconn@waterloo.csnet) --------------- [3] Richard H. Byrd and Robert B. Schnabel Continuity of the Null Space Basis and Constrained Optimization Univ. of Colorado Comp. Sci. Report CU-CS-271-84 Many constrained optimization algorithms use a basis for the null space of the matrix of constraint gradients. Recently, methods have been proposed that enable this null space basis to vary continuously as a function of the iterates in a neighborhood of the solution. This paper reports results from topology showing that, in general, there is no continuous function that generates the null space basis of all full rank rectangular matrices of a fixed size. We also give some indication of where these discontinuities must occur. We then propose an alternative implementation of a class of constrained optimization algorithms that uses approximations to the reduced Hessian of the Lagrangian but is independent of the choice of null space basis. Submitted by: Richard Byrd Obtainable from: Richard Byrd Department of Computer Science University of Colorado Boulder, Colorado 80309 (richard@boulder.csnet) --------------- [4] Y. Saad, A. Sameh, and P. Saylor Solving Elliptic Difference Equations on a Linear Array of Processors. In this paper we consider the organization of three iterative methods for solving self adjoint elliptic difference on a set of linearly con- nected processors. These algorithms are the cyclic Chebyshev semi-iterative scheme, a preconditioned conjugate gradient method, and a generalization of the Chebyshev method. We also consider their performance on this multi- processor as a function of the cost of interprocessor communication. Submitted by: Paul Saylor Obtainable from: Paul Saylor Computer Science Dept. University of Illinois Urbana-Champaine, Illinois (saylor@uiuc.csnet) --------------- [5] Y. Saad, A. Sameh, and P. Saylor An Optimum Semi-iterative Method for Solving Any Linear Set With a Square Matrix. An algorithm is presented to solve Ax=b by computing optimum iteration parameters for Richardson's iteration. The algorithm re- quires some information on the location of the eigenvalues of A. The algorithm yields parameters well-suited for matrices for which Chebyshev parameters are not appropriate. It therefore supplements the Manteuffel algorithm, developed for the Chebyshev case. Numerical examples are described. Submitted by: Paul Saylor Obtainable from: Paul Saylor (as above) --------------- [6] Y. Saad, A. Sameh, and P. Saylor A Generalization of the Golub-Welsch Algorithm to Complex Orthogonal and Kernel Polynomials. The Golub-Welsch algorithm is a stable algorithm to compute the roots of real orthogonal polynomials and kernel polynomials. It is generalized to the complex cases. Submitted by: Paul Saylor Obtainable from: Paul Saylor (as above) --------------- [7] R.T. Gregory A finite number system for digital computers which gives exact results with rational operands. University of Tennessee Computer Science Report CS-84-57 A mapping is established between those rational numbers a/b with both a and b bounded in absolute value by a constant N, and a finite set of integers. When N is chosen to be the largest integer satisfying the inequality 2N*2 + 1 <= p where p is a very large prime, then the mapping is one-to-one and onto the set (0,1,2,...,p-1). Computation involving the rational numbers can be replaced by computation on the integers modulo p. In practice more than one modulus is used and this report deals with the multi- modulus case. Submitted by: Robert T. Gregory (shortly before his death) Obtainable from: Computer Science Department University of Tennessee Knoxville, Tenn. --------------- [8] Daniel Boley, University of Minnesota A Parallel Method for the Generalized Eigenvalue Problem University of Minnesota Computer Science Tech Report 84-21 Sept 1984. We present a parallel method to solve the generalized eigenvalue problem on a linear array of processors, each connected to their nearest neighbors and operating synchronously. We also include a wraparound connection from end to end. Our method is based on the well-known QZ algorithm of Moler and Stewart, which simultaneously reduces two nxn matrices to upper triangular form by orthogonal or unitary tranformations. We show how this algorithm may be partitioned and distributed over n+1 processors, achieving a speed-up over the serial algorithm of O(n). We use the concept of windows to describe the action of each processor at each step. We show how to incorporate single shifts, and how to apply orthogonal plane rota- tions on either side of a matrix without the need to transpose the matrix itself. Submitted by: Daniel Boley Obtainable from: Daniel Boley University of Minnesota (boley@umn-cs.csnet) --------------- [9] Ito, Kazufumi An iterative method for indefinite systems of linear equations ICASE Report No. 84-13, April 2, 1984, 21 pages Submitted to SIAM J. Numer. Anal. An iterative method for solving nonsymmetric indefinite linear systems is proposed. The method involves the successive use of a modified version of the conjugate residual method. A numerical example is given to illustrate the method. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft ICASE, Mail Stop l32C NASA Langley Research Center Hampton, VA 23665 (804) 865-3992. --------------- [10] Adams, Loyce M. and Harry F. Jordan Is SOR color-blind? ICASE Report No. 84-14, May 21, 1984, 39 pages Submitted to SIAM J. Sci. Statist. Comput. The work of Young [1971] showed that the Red/Black ordering and the natural rowwise ordering of matrices with Property A, such as those arising from the 5-point discretization of Poisson's equation, lead to SOR iteration matrices with identical eigenvalues. With the advent of parallel computers, multi-color point SOR schemes have been proposed for more complicated stencils on 2-dimensional rectangular grids, see Adams and Ortega [1982] for example, but to our knowledge, no theory has been provided for the rate of convergence of these methods relative to that of the natural rowwise scheme. New results show that certain matrices may be reordered so the resulting multi-color SOR matrix has the same eigenvalues as that for the original ordering. In addition, for a wide range of stencils, we show how to choose multi-color orderings so the multi-color SOR matrices have the same eigenvalues as the natural rowwise SOR matrix. The strategy for obtaining these results is based on "data flow" concepts and can be used to reach Young's conclusions above for the 5-point stencil. The importance of these results is threefold. Firstly, a constructive and easy means of finding these multi-colorings is a direct consequence of the theory; secondly, multi-color SOR methods can be found that have the same rate of convergence as the natural rowwise SOR method for a wide range of stencils used to discretize partial differential equations; and thirdly, these multi- color SOR methods can be efficiently implemented on a wide class of parallel computers. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [11] Gatski, Thomas B. and Chester E. Grosch A numerical study of the two- and three-dimensional unsteady Navier-Stokes equations in velocity-vorticity variables using compact difference schemes ICASE Report No. 84-15, May 22, 1984, 20 pages Proc. of 9th ICNMED, Saclay, France, June 23-29, 1984 LECTURE NOTES IN PHYSICS Springer-Verlag. A compact finite-difference approximation to the unsteady Navier-Stokes equations in velocity-vorticity variables is used to numerically simulate a number of flows. These include two-dimensional laminar flow of a vortex evolving over a flat plate with an embedded cavity, the unsteady flow over an elliptic cylinder, and aspects of the transient dynamics of the flow over a rearward facing step. The methodology required to extend the two-dimensional formulation to three-dimensions is presented. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [12] Hafez, Mohamed, M. and Phillips, Timothy N. A modified least squares formulation for a system of first-order equations ICASE Report No. 84-16, May 29, 1984, 21 pages Submitted to IMACS J. Numer. Math. Second-order equations in terms of auxiliary variables similar to potential and stream functions are obtained by applying a weighted least squares formulation to a first-order system. The additional boundary conditions which are necessary to solve the higher order equations are determined and numerical results are presented for the Cauchy-Riemann equations. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [13] Hall, Phillip The Gortler vortex instability mechanism in three-dimensional boundary layers ICASE Report No. 84-17, June 6, 1984, 37 pages. Submitted to the Royal Society of London. It is well known that the two-dimensional boundary layer on a concave wall is centrifugally unstable with respect to vortices aligned with the basic flow for sufficiently high values of the Gortler number. However, in most situations of practical interest the basic flow is three-dimensional and previous theoretical investigations do not apply. In this paper the linear stability of the flow over an infinitely long swept wall of variable curvature is considered. If there is no pressure gradient in the boundary layer it is shown that the instability problem can always be related to an equivalent two- dimensional calculation. However, in general, this is not the case and even for small values of the crossflow velocity field dramatic differences between the two and three-dimensional problems emerge. In particular, it is shown that when the relative size of the crossflow and chordwise flow is (R**(-l/2)) where R is the Reynolds number of the flow, the most unstable mode is time- dependent. When the size of the crossflow is further increased, the vortices in the neutral location have their axes locally perpendicular to the vortex lines of the basic flow. In this regime the eigenfunctions associated with the instability become essentially 'centre modes' of the Orr-Sommerfeld equation destabilized by centrifugal effects. The critical Gortler number for such modes can be predicted by a large wavenumber asymptotic analysis; the results suggest that for order unity values of the ratio of the crossflow and chordwise velocity fields, the Gortler instability mechanism is almost certainly not operational. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [14] Gozani, J., Nachshon, A., and Turkel, E. Conjugate gradient coupled with multigrid for an indefinite problem ICASE Report No. 84-18, June 7, 1984, 11 pages Submitted to Proc. of 3rd IMACS International Symposium on Computer Methods for Partial Differential Equations Lehigh University, Bethlehem, PA, June 20-22, 1984 An iterative algorithm for the Helmholtz equation is presented. This scheme is based on the preconditioned conjugate gradient method for the normal equations. The preconditioning is one cycle of a multigrid method for the discrete Laplacian. The smoothing algorithm is red-black Gauss-Seidel and is constructed so it is a symmetric operator. The total number of iterations needed by the algorithm is independent of h. By varying the number of grids, the number of iterations depends only weakly on k when k**3 h**2 is constant. Comparisons with a SSOR preconditioner are presented. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [15] Malik, M. R., Zang, T. A., and Hussaini, M. Y. A spectral collocation method for the Navier-Stokes equations ICASE Report No. 84-19, June 8, 1984, 48 pages Submitted to J. Comput. Phys. A Fourier-Chebyshev spectral method for the incompressible Navier-Stokes equations is described. It is applicable to a variety of problems including some with fluid properties which vary strongly both in the normal direction and in time. In this fully spectral algorithm, a preconditioned iterative technique is used for solving the implicit equations arising from semi- implicit treatment of pressure, mean advection and vertical diffusion terms. The algorithm is tested by applying it to hydrodynamic stability problems in channel flow and in external boundary layers with both constant and variable viscosity. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [16] Davis, Stephen F. TVD finite difference schemes and artificial viscosity. ICASE Report No. 84-20, June 11, 1984, 35 pages Submitted to SIAM J. Sci. Statis. Comput. In this paper we show that the total variation diminishing (TVD) finite difference scheme which was analysed by Sweby can be interpreted as a Lax- Wendroff scheme plus an upwind weighted artificial dissipation term. We then show that if we choose a particular flux limiter and remove the requirement for upwind weighting, we obtain an artificial dissipation term which is based on the theory of TVD schemes, which does not contain any problem dependent parameters and which can be added to existing MacCormack method codes. Finally, we conduct numerical experiments to examine the performance of this new method. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [17] Canuto, C., Hariharan, S. I., and Lustman, L. Spectral methods for exterior elliptic problems ICASE Report No. 84-21, June 12, 1984, 33 pages. Submitted to Numer. Math. This paper deals with spectral approximations for exterior elliptic problems in two dimensions. As in the conventional finite difference or finite element methods, it is found that the accuracy of the numerical solutions is limited by the order of the numerical farfield conditions. We introduce a spectral boundary treatment at infinity, which is compatible with the "infinite order" interior spectral scheme. Computational results are presented to demonstrate the spectral accuracy attainable. Although we deal with a simple Laplace problem throughout the paper, our analysis covers more complex and general cases. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [18] Graif, E. and K. Kunisch Parameter estimation for the Euler-Bernoulli-beam ICASE Report No. 84-22, June 20, 1984, 43 pages. Submitted to IEEE Trans. Auto. Control. An approximation involving cubic spline functions for parameter estimation problems in the Euler-Bernoulli-beam equation (phrased as an optimization problem with respect to the parameters) is described and convergence is proved. The resulting algorithm was implemented and several of the test examples are documented. It is observed that the use of penalty terms in the cost functional can improve the rate of convergence. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [19] Rosen, I. Gary A numerical scheme for the identification of hybrid systems describing the vibration of flexible beams with tip bodies ICASE Report No. 84-23, June 21, 1984, 36 pages. Submitted to J. Math. Anal. Appl. A cubic spline based Galerkin-like method is developed for the identification of a class of hybrid systems which describe the transverse vibration of flexible beams with attached tip bodies. The identification problem is formulated as a least squares fit to data subject to the system dynamics given by a coupled system of ordinary and partial differential equations recast as an abstract evolution equation (AEE) in an appropriate infinite dimensional Hilbert space. Projecting the AEE into spline-based subspaces leads naturally to a sequence of approximating finite diemnsional identification problems. The solutions to these problems are shown to exist, are relatively easily computed, and are shown to, in some sense, converge to solutions to the original identification problem. Numerical results for a variety of examples are discussed. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [20] Banks, H. Thomas and Katherine A. Murphy Estimation of coefficients and boundary parameters in hyperbolic systems ICASE Report No. 84-24, June 21, 1984, 52 pages. Submitted to SIAM J. Control Optim. We consider semi-discrete Galerkin approximation schemes in connection with inverse problems for the estimation of spatially varying coefficients and boundary condition parameters in second order hyperbolic systems typical of those arising in 1-D surface seismic problems. Spline based algorithms are proposed for which theoretical convergence results along with a representative sample of numerical findings are given. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [21] Vardi, A. A new minmax problem ICASE Report No. 84-25, June 26, 1984, 35 pages. Submitted to J. Optim., Theory and Appl. For the classical minimax problem we design a new active set strategy in which there are three types of functions: active, semi-active, and non-active. This technique will help in preventing zigzagging which often occurs when an active set strategy is used. Some of the inequality constraints are handled with slack variables. Also a trust region strategy is used in which at each iteration there is a sphere around the current point in which the local approximation of the function is trusted. The algorithm suggested in the paper was implemented into a successful computer program. Numerical results are provided. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [22] Banks, H. T. and I. G. Rosen Approximation techniques for parameter estimation and feedback control for distributed models of large flexible structures ICASE Report No. 84-26, June 22, 1984, 19 pages. Submitted to Workshop on Identification and Control of Flexible Space Structures, San Diego, CA, June 4-6, 1984, G. Rodgriguez (ed.). We discuss approximation ideas that can be used in parameter estimation and feedback control for Euler-Bernoulli models of elastic systems. Focusing on parameter estimation problems, we outline how one can obtain convergence results for cubic spline-based schemes for hybrid models involving an elastic cantilevered beam with tip mass and base acceleration. Sample numerical findings are also presented. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [23] Turkel, E. and Bram van Leer Flux vector splitting and Runge-Kutta methods for the Euler equations ICASE Report No. 84-27, June 23, 2984, 9 pages. To appear in Proc. of the 9th Intl. Conference on Numerical Methods in Fluid Dynamics, Saclay, France, June 25-29, 1984. Runge-Kutta schemes have been used as a method of solving the Euler equations exterior to an airfoil. In the past this has been coupled with central differences and an artificial viscosity in space. In this study we couple the Runge-Kutta time stepping scheme with an upwinded space approxima tion based on flux-vector splitting. Several acceleration techniques are also considered including a local time step, residual smoothing and multigrid. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [24] Turkel, Eli Fast solutions to the steady state compressible and incompressible fluid dynamics equations ICASE Report No. 84-28, June 23, 1984, 8 pages. To appear in Proc. of the 9th Intl. Conference on Numerical Methods in Fluid Dynamics Saclay, France, June 25-29, 1984. It is well known that for low speed flows the use of the compressible fluid dynamic equations is inefficient. The use of an explicit scheme requires delta-t to be bounded by 1/c. However, the physical parameters change over time scales of order 1/u which is much larger. Hence, it is not appropriate to use explicit schemes for very subsonic flows. Implicit schemes are hard to vectorize and frequently do not converge quickly for very subsonic flows. We shall demonstrate that if one is only interested in the steady state then a minor change to an existing code can greatly increase the efficiency of an explicit method. Even when using an implicit method the proposed changes increase the efficiency of the scheme. We shall first consider the Euler equations for low speed flows and then incompressible flows. We then indicate how to generalize the method to include viscous effects. We also show how to accelerate supersonic flow by essentially decoupling the equations. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [25] Gottlieb, D. Spectral methods for compressible flow problems ICASE Report No. 84-29, June 21, 1984, 20 pages. Proc. 9th Intl. Conference on Numerical Methods in Fluid Dynamics Saclay, France, June 25-29, 1984. In this article we review recent results concerning numerical simulation of shock waves using spectral methods. We discuss shock fitting techniques as well as shock capturing techniques with finite difference artificial viscosity. We also discuss the notion of the information contained in the numerical results obtained by spectral methods and show how this information can be recovered. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [26] Ipsen, Ilse C. F. Singular value decomposition with systolic arrays ICASE Report No. 84-30, July 9, 1984, 22 pages. Proc. of the Society of Photo-Optical Engineers San Diego, CA, August 19-24, 1984. Systolic arrays for determining the Singular Value Decomposition of a matrix A of bandwidth w are presented. After A has been reduced to bidiagonal form B by means of Givens plane rotations, the differential equations. Part II: The linear quadratic control problem. ICASE Report No. 84-31, July 6, 1984, 65 pages. Submitted to SIAM J. Control Optim. The numerical scheme based on the Legendre-tau approximation is proposed to approximate the feedback solution to the linear quadratic optimal control problem for hereditary differential systems. The convergence property is established using Trotter ideas. The method yields very good approximations at low orders and provides an approximation technique for computing closed- loop eigenvalues of the feedback system. A comparison with existing methods (based on "averaging" and "spline" approximations) is made. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [27] Turkel, Eli Acceleration to a steady state for the Euler equations ICASE Report No. 84-32, July 7, 1984, 45 pages. Proc. INRIA Euler Workshop (SIAM) Rocquencourt, France, December 7-9, 1983. A multi-stage Runge-Kutta method is analyzed for solving the Euler equations exterior to an airfoil. Highly subsonic, transonic and supersonic flows are evaluated. Various techniques for accelerating the convergence to a steady state are introduced and analyzed. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [28] Bokhari, Shahid H. and A. D. Raza Augmenting computer networks ICASE Report No. 84-33, July 8, 1984, 16 pages. Proc. 1984 Intl. Conference on Parallel Processing Bellaire, MI, August 21-24, 1984. Three methods of augmenting computer networks by adding at most one link per processor are discussed. 1. A tree of N nodes may be augmented such that the resulting graph has diameter no greater than 4 log2((N+2)/3)-2. This 0(N**3) algorithm can be applied to any spanning tree of a connected graph to reduce the diameter of that graph to 0(log N). 2. Given a binary tree T and a chain C of N nodes each, C may be augmented to produce C' so that T is a subgraph of C'. This algorithm is 0(N) and may be used to produce augmented chains or rings that have diameter no greater than 2 log2 (N+2)/3) and are planar. 3. Any rectangular two-dimensional 4 (8) nearest neighbor array of size N = 2**k may be augmented so that it can emulate a single step shuffle- exchange network of size N/2 in 3(t) time steps. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [29] Reed, Daniel A. and Merrell L. Patrick A model of asynchronous iterative algorithms for solving large, sparse, linear systems ICASE Report No. 84-34, July 31, 1984, 27 pages. Proceedings of the 1984 International Conference on Parallel Processing August 21-24, 1984, Bellaire, MI. Solving large, sparse, linear systems of equations is one of the fundamental problems in large scale scientific and engineering computation. A model of a general class of asynchronous, iterative solution methods for linear systems is developed. In the model, the system is solved by creating several cooperating tasks that each compute a portion of the solution vector. This model is then analyzed to determine the expected intertask data transfer and task computational complexity as functions of the number of tasks. Based on the analysis, recommendations for task partitioning are made. These recommendations are a function of the sparseness of the linear system, its structure (i.e., randomly sparse or banded), and dimension. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [30] Reed, Daniel A. and Merrell L. Patrick Parallel, iterative solution of sparse linear systems: Models and architectures ICASE Report No. 84-35, August 1, 1984, 45 pages. Submitted to Parallel Computing. Solving large, sparse, linear systems of equations is a fundamental problem in large scale scientific and engineering computation. A model of a general class of asychronous, iterative solution methods for linear systems is developed. In the model, the system is solved by creating several cooperating tasks that each compute a portion of the solution vector. A data transfer model predicting both the probability that data must be transferred between two tasks and the amount of data to be transferred is presented. This model is used to derive an execution time model for predicting parallel execution time and an optimal number of tasks given the dimension and sparsity of the coefficient matrix and the costs of computation, synchronization, and communication. The suitability of different parallel architectures for solving randomly sparse linear systems is discussed. Based on the complexity of task scheduling, one parallel architecture, based on a broadcast bus, is presented and analyzed. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [31] Tadmor, Eitan The well-posedness of the Kuramoto-Sivashinsky equation. ICASE Report No. 84-36, August 7, 1984, 22 pages. Submitted to SIAM J. Appl. Math. The Kuramoto-Sivashinsky equation arises in a variety of applications, among which are modeling reaction-diffusion systems, flame-propagation and viscous flow problems. It is considered here, as a prototype to the larger class of generalized Burgers equations: those consist of quadratic nonlinearity and arbitrary linear parabolic part. We show that such equations are well-posed, thus admitting a unique smooth solution, continuously dependent on its initial data. As an attractive alternative to standard energymethods, existence and stability are derived in this case, by "patching" in the large short time solutions without "loss of derivatives". Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [32] Lustman, Liviu The time evolution of spectral discretizations of hyperbolic systems ICASE Report No. 84-37, August 7, 1984, 12 pages. Submitted to SIAM J. Numer. Anal. A Chebyshev collocation spectral method, applied to hyperbolic systems is considered, particularly for those initial boundary value problems which possess only solutions tending to zero at large times. It is shown that the numerical solutions of the system will also vanish at infinity, if and only if, the numerical solution of a scalar equation of the same type does. This result is then generalized for other spectral approximations. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [33] Bayliss, A., C. I. Goldstein, and E. Turkel On accuracy conditions for the numerical computation of waves ICASE Report No. 84-38, August 9, 1984, 16 pages. Submitted to J. Comput. Phys. The Helmholtz equation with a variable index of refraction, and a suitable radiation condition at infinity serves as a model for a wide variety of wave propagation problems. Such problems can be solved numerically by first truncating the given unbounded domain and imposing a suitable outgoing radiation condition on an artificial boundary and then solving the resulting problem on the bounded domain by direct discretization (for example, using a finite element method). In practical applications, the mesh size h and the wave number K, are not independent but are constrained by the accuracy of the desired computation. It will be shown that the number of points per wavelength, measured by l/Kh, is not sufficient to determine the accuracy of a given discretization. For example, the quantity K**3 h**2 is shown to determine the accuracy in the L-two norm for a second-order discretization method applied to several propagation models. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [34] Hariharan, S. I. Numerical solutions of acoustic wave propagation problems using Euler computations ICASE Report No. 84-39, August 10, 1984, 29 pages. Proc. of the AIAA 9th Aeroacoustics Conference Williamsburg, VA, August 15-17, 1984. This paper reports solution procedures for problems arising from the study of engine inlet wave propagation. The first problem is the study of sound waves radiated from cylindrical inlets. The second one is a quasi-one- dimensional problem to study the effect of nonlinearities and the third one is the study of nonlinearities in two dimensions. In all three problems Euler computations are done with a fourth-order explicit scheme. For the first problem results are shown in agreement with experimental data and for the second problem comparisons are made with an existing asymptotic theory. The third problem is part of an ongoing work and preliminary results are presented for this case. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [35] Tadmor, Eitan The exponential accuracy of Fourier and Tchebyshev differencing methods ICASE Report No. 84-40, August 20, 1984, 24 pages. Submitted to SIAM J. Numer. Anal. It is shown that when differencing analytic functions using the pseudospectral Fourier or Tchebyshev methods, the error committed decays to zero at an exponential rate. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [36] Gannon, Dennis and John Van Rosendale On the impact of communication complexity in the design of parallel numerical algorithms ICASE Report No. 84-41, August 22, 1984, 37 pages. To appear in IEEE Trans. on Computers. This paper describes two models of the cost of data movement in parallel numerical algorithms. One model is a generalization of an approach due to Hockney, and is suitable for shared memory multiprocessors where each processor has vector capabilities. The other model is applicable to highly parallel nonshared memory MIMD systems. In the second model, algorithm performance is characterized in terms of the communication network design. Techniques used in VLSI complexity theory are also brought in, and algorithm independent upper bounds on system performance are derived for several problems that are important to scientific computation. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [37] Ito, Kazufumi Legendre-Tau approximation for functional differential equations, Part III: Eigenvalue approximations and uniform stability. ICASE Report No 84-42, August 24, 1984, 34 pages. Proc. Second International Conference on Control Theory for Distributed Parameter Systems and Applications Vorau, Austria, 1984. The stability and convergence properties of the Legendre-tau approximation for hereditary differential systems are analyzed. We derive a characteristic equation for the eigenvalues of the resulting approximate system. As a result of this derivation we are able to establish that the uniform exponential stability of the solution semigroup is preserved under approximation. It is the key to obtaining the convergence of approximate solutions of the algebraic Riccati equation in trace norm. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [38] Berger, M. J. On conservation at grid interfaces ICASE Report No. 84-43, September 7, 1984, 25 pages. Submitted to J. Comput. Phys. This paper considers the solution of hyperbolic systems of conservation laws on discontinuous grids. In particular, we consider what happens to conservation at grid interfaces. A procedure is presented to derive conservative difference approximations at the grid interfaces for two- dimensional grids which overlap in an arbitrary configuration. The same procedures are applied to compute interface formulas for grids which are refined in space and/or time, and for continuous grids where a switch in the scheme causes the discontinuity. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [39] Osher, S. and S. Chakravarthy Very high order accurate TVD schemes ICASE Report No. 84-44, September 10, 1984, 61 pages. Submitted to SIAM J. Numer. Anal. A systematic procedure for constructing semi-discrete families of 2m-1 order accurate, 2m order dissipative, variation diminishing, 2m+1 point bandwidth, conservation form approximations to scalar conservation laws is presented. Here m is any integer between 2 and 8. Simple first-order forward time discretization, used together with any of these approximations to the space derivatives, also results in a fully discrete, variation diminishing algorithm. These schemes all use simple flux limiters, without which each of these fully discrete algorithms is even linearly unstable. Extensions to systems, using a nonlinear field-by-field decomposition are presented, and shown to have many of the same properties as in the scalar case. For linear systems, these nonlinear approximations are variation diminishing, and hence convergent. A new and general criterion for approximations to be variation diminishing is also given. Finally, numerical experiments using some of these algorithms are presented. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [40] Elcrat, A. R. and L. N. Trefethen Classical free-streamline flow over a polygonal obstacle ICASE Report No. 84-45, September 10, 1984, 61 pages. Submitted to J. Comput. Appl. Math. In classical Kirchhoff flow, an ideal incompressible fluid flows past an obstacle and around a motionless wake bounded by free streamlines. Since 1869 it has been known that in principle, the two-dimensional Kirchhoff flow over a polygonal obstacle can be determined by constructing a conformal map onto a polygon in the log-hodograph plane. In practice, however, this idea has rarely been put to use except for very simple obstacles, because the conformal mapping problem has been too difficult. This paper presents a practical method for computing flows over arbitrary polygonal obstacles to high accuracy in a few seconds of computer time. We achieve this high speed and flexibility by working with a modified Schwarz-Christoffel integral that maps onto the flow region directly rather than onto the log-hodograph polygon. This integral and its associated parameter problem are treated numerically by methods developed earlier by Trefethen for standard Schwarz-Christoffel maps. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [41] Krause, Egon Review of some vortex relations ICASE Report No. 84-46, September 12, 1984, 8 pages. Submitted to Computer and Fluids. The evaluation of the circulation from numerical solutions of the momentum and energy equations is discussed for incompressible and compressible flows. It is shown how artificial damping directly influences the time rate of change of the circulation. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [42] Tal-Ezer, Hillel A pseudospectral Legendre method for hyperbolic equations with an improved stability condition ICASE Report No. 84-48, September 14, 1984, 46 pages. Submitted to J. Comput. Phys. A new pseudospectral method is introduced for solving hyperbolic partial differential equations. This method uses different grid points than previously used pseudospectral methods: in fact the grid points are related to the zeroes of the Legendre polynomials. The main advantage of this method is that the allowable time step is proportional to the inverse of the number of grid points 1/N rather than to 1/N**2 (as in the case of other pseudospectral methods applied to mixed initial boundary value problems). A highly accurate time discretization suitable for these spectral methods is discussed. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [43] Bayliss, A., C. I. Goldstein, and E. Turkel The numerical solution of the Helmholtz equation for wave propagation problems in underwater acoustics ICASE Report No. 84-49, September 17, 1983, 29 pages. Submitted to J. Comput. Math. Appl. The Helmholtz Equation with a variable index of refraction, and a suitable radiation condition at infinity serves as a model for a wide variety of wave propagation problems. A numerical algorithm has been developed and a computer code implemented that can effectively solve this equation in the intermediate frequency range. The equation is discretized using the finite element method, thus allowing for the modeling of complicated geometries (including interfaces) and complicated boundary conditions. A global radiation boundary condition is imposed at the far field boundary that is exact for an arbitrary number of propagating modes. The resulting large, non-selfadjoint system of linear equations with indefinite symmetric part is solved using the preconditioned conjugate gradient method applied to the normal equations. A new preconditioner is developed based on the multigrid method. This preconditioner is vectorizable and is extremely effective over a wide range of frequencies provided the number of grid levels is reduced for large frequencies. A heuristic argument is given that indicates the superior convergence properties of this preconditioner. The relevant limit to analyze convergence is for K increasing and a fixed prescribed accuracy level. The efficiency and robustness of the numerical algorithm is confirmed for large acoustic models, including interfaces with strong velocity contrasts. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [44] Burns, John A., Ito, Kazufumi, Powers, Robert K. Chandrasekhar equations and computational algorithms for distributed parameter systems ICASE Report No. 84-50, September 20, 1984, 23 pages. Proc. of 23rd Conference on Decision and Control, December 12-14, 1984, Las Vegas, NV. In this paper we consider the Chandrasekhar equations arising in optimal control problems for linear distributed parameter systems. The equations are derived via approximation theory. This approach is used to obtain existence, uniqueness and strong differentiability of the solutions and provides the basis for a convergent computation scheme for approximating feedback gain operators. A numerical example is presented to illustrate these ideas. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [45] Petter E. Bjorstad and Olof B. Widlund Iterative Methods for the Solution of Elliptic Problems on Regions Partitioned into Substructures Finite element problems can often naturally be divided into subproblems which correspond to subregions into which the region has been partitioned or from which it was originally assembled. A class of iterative methods is discussed in which these subproblems are solved by direct methods,while the interaction across the curves or surfaces which divide the region is handled by a conju- gate gradient method. A mathematical framework for this workis provided by regularity theory for elliptic finite element problems and by block Gauss elimination. A full development of the theory,which shows that certain of these methods are optimal, is given for Lagrangian finite element approxi- mations of second order linear elliptic problems in the plane.Results from numerical experiments are also reported. Submitted by: Olaf Widlund Obtainable from: Olaf Widlund Courant Institute 251 Mercer St. N.Y., N.Y. 10012. (WIDLUND@NYU-ACF1.ARPA) --------------- [46] Ron S. Dembo and Ulrich Tulowitzki Local Convergence Analysis for Successive Inexact Quadratic Programming Methods SOM Working Paper Series B #78 Yale School of Organization and Management October, 1984 Successive quadratic programming (SQP) methods for the solution of constrained nonlinear optimization problems are attractive because they converge rapidly from any sufficiently good initial solution. However, solving the quadratic subproblems at each stage can be expensive, particularly if the number of unknowns is large and may be not justified when far away from an optimal solution. Therefore, we consider a class of successive inexact quadratic programming (SIQP) methods that solve the subproblems only approximately. We characterize the rate of convergence in terms of the relative error incurred in the subproblems in a manner that does not assume that the set of active constraints remains constant. This leads to practical termination criteria for truncating various iterative algorithms applied to the quadratic subproblems. Submitted by: Ron Dembo Obtainable from: Ron Dembo Yale University School of Organization and Management P.O. Box 1A New Haven, Conn. 06520 --------------- [47] Ron S. Dembo and Ulrich Tulowitzki Sequential Truncated Quadratic Programming Methods To appear in the Proceedings of the SIAM Conference on Numerical Optimization, Boulder, Colo., June, 1984 In a sequential quadratic programming (SQP) method for nonlinear optimization it is possible to truncate an iterative procedure, applied to each quadratic subproblem, prior to reaching a solution, without affecting the asymptotic rate of convergence. We describe some practical termination criteria which ensure that an SQP method, executed without exact solution of the QP subproblem, inherits the rate of convergence of its (exact) SQP counterpart. Numerical experiments on linearly-constrained problems of up to 2230 variables indicate that substantial savings can be achieved by inexact solution of the QP subproblems with little or no adverse effect on the convergence rate. Submitted by: Ron Dembo Obtainable from: Ron Dembo Yale University School of Organization and Management P.O. Box 1A New Haven, Conn. 06520 --------------- From rhbartels%watcgl%waterloo.csnet@csnet-relay.arpa Tue Jul 16 13:45:20 1985 Received: from csnet-relay (csnet-relay.arpa.ARPA) by anl-mcs.ARPA ; Tue, 16 Jul 85 13:42:53 cdt Received: from waterloo by csnet-relay.csnet id a001168; 16 Jul 85 14:21 EDT Date: Sun, 14 Jul 85 20:30:55 edt From: Richard Bartels <rhbartels%watcgl%waterloo.csnet@csnet-relay.arpa> Message-Id: <8507150030.AA20616@watcgl.UUCP> To: dongarra@anl-mcs Status: R ================================================== NUMERICAL ANALYSIS TITLES -- Volume 2, Number 4 November 23, 1984 ================================================== [1] Finbarr O'Sullivan and Grace Wahba A Cross Validated Bayesian Retrieval Algorithm For Non-Linear Remote Sensing Experiments To appear, Journal of Computational Physics We present a fully automated retrieval algorithm for estimating the solution of a mildly non linear Fredholm Integral Equation of the First Kind, when the data are discrete, noisy (experimental) measurements. This algorithm combines a simple Gauss-Newton interation with an extended form of Generalized Cross Validation , applied in the context of Tihonov Regularization with moment discretization. The performance of the algorithm is illustrated in the context of a remote sensing problem arising in satellite meteorology. Submitted by: Grace Wahba Obtainable from: Grace Wahba Department of Statistics University of Wisconsin, Madison (na.wahba@su-score.arpa) --------------- [2] ANDREW R. CONN NONLINEAR PROGRAMMING, EXACT PENALTY FUNCTIONS AND PROJECTION TECHNIQUES FOR NON-SMOOTH FUNCTIONS We present a personal overview of various approaches to solving nonlinear programs with nonlinear constraints that make use of the L-1 exact penalty function. The advantages, disadvantages and related remaining difficulties of these approaches will be considered. Finally some recent research and extensions are given. Submitted by: Andrew R. Conn Obtainable from: Andrew R. Conn Computer Science Department University of Waterloo Waterloo, Ontario N2L 3G1 CANADA (arconn@waterloo.csnet) --------------- [3] Richard H. Byrd and Robert B. Schnabel Continuity of the Null Space Basis and Constrained Optimization Univ. of Colorado Comp. Sci. Report CU-CS-271-84 Many constrained optimization algorithms use a basis for the null space of the matrix of constraint gradients. Recently, methods have been proposed that enable this null space basis to vary continuously as a function of the iterates in a neighborhood of the solution. This paper reports results from topology showing that, in general, there is no continuous function that generates the null space basis of all full rank rectangular matrices of a fixed size. We also give some indication of where these discontinuities must occur. We then propose an alternative implementation of a class of constrained optimization algorithms that uses approximations to the reduced Hessian of the Lagrangian but is independent of the choice of null space basis. Submitted by: Richard Byrd Obtainable from: Richard Byrd Department of Computer Science University of Colorado Boulder, Colorado 80309 (richard@boulder.csnet) --------------- [4] Y. Saad, A. Sameh, and P. Saylor Solving Elliptic Difference Equations on a Linear Array of Processors. In this paper we consider the organization of three iterative methods for solving self adjoint elliptic difference on a set of linearly con- nected processors. These algorithms are the cyclic Chebyshev semi-iterative scheme, a preconditioned conjugate gradient method, and a generalization of the Chebyshev method. We also consider their performance on this multi- processor as a function of the cost of interprocessor communication. Submitted by: Paul Saylor Obtainable from: Paul Saylor Computer Science Dept. University of Illinois Urbana-Champaine, Illinois (saylor@uiuc.csnet) --------------- [5] Y. Saad, A. Sameh, and P. Saylor An Optimum Semi-iterative Method for Solving Any Linear Set With a Square Matrix. An algorithm is presented to solve Ax=b by computing optimum iteration parameters for Richardson's iteration. The algorithm re- quires some information on the location of the eigenvalues of A. The algorithm yields parameters well-suited for matrices for which Chebyshev parameters are not appropriate. It therefore supplements the Manteuffel algorithm, developed for the Chebyshev case. Numerical examples are described. Submitted by: Paul Saylor Obtainable from: Paul Saylor (as above) --------------- [6] Y. Saad, A. Sameh, and P. Saylor A Generalization of the Golub-Welsch Algorithm to Complex Orthogonal and Kernel Polynomials. The Golub-Welsch algorithm is a stable algorithm to compute the roots of real orthogonal polynomials and kernel polynomials. It is generalized to the complex cases. Submitted by: Paul Saylor Obtainable from: Paul Saylor (as above) --------------- [7] R.T. Gregory A finite number system for digital computers which gives exact results with rational operands. University of Tennessee Computer Science Report CS-84-57 A mapping is established between those rational numbers a/b with both a and b bounded in absolute value by a constant N, and a finite set of integers. When N is chosen to be the largest integer satisfying the inequality 2N*2 + 1 <= p where p is a very large prime, then the mapping is one-to-one and onto the set (0,1,2,...,p-1). Computation involving the rational numbers can be replaced by computation on the integers modulo p. In practice more than one modulus is used and this report deals with the multi- modulus case. Submitted by: Robert T. Gregory (shortly before his death) Obtainable from: Computer Science Department University of Tennessee Knoxville, Tenn. --------------- [8] Daniel Boley, University of Minnesota A Parallel Method for the Generalized Eigenvalue Problem University of Minnesota Computer Science Tech Report 84-21 Sept 1984. We present a parallel method to solve the generalized eigenvalue problem on a linear array of processors, each connected to their nearest neighbors and operating synchronously. We also include a wraparound connection from end to end. Our method is based on the well-known QZ algorithm of Moler and Stewart, which simultaneously reduces two nxn matrices to upper triangular form by orthogonal or unitary tranformations. We show how this algorithm may be partitioned and distributed over n+1 processors, achieving a speed-up over the serial algorithm of O(n). We use the concept of windows to describe the action of each processor at each step. We show how to incorporate single shifts, and how to apply orthogonal plane rota- tions on either side of a matrix without the need to transpose the matrix itself. Submitted by: Daniel Boley Obtainable from: Daniel Boley University of Minnesota (boley@umn-cs.csnet) --------------- [9] Ito, Kazufumi An iterative method for indefinite systems of linear equations ICASE Report No. 84-13, April 2, 1984, 21 pages Submitted to SIAM J. Numer. Anal. An iterative method for solving nonsymmetric indefinite linear systems is proposed. The method involves the successive use of a modified version of the conjugate residual method. A numerical example is given to illustrate the method. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft ICASE, Mail Stop l32C NASA Langley Research Center Hampton, VA 23665 (804) 865-3992. --------------- [10] Adams, Loyce M. and Harry F. Jordan Is SOR color-blind? ICASE Report No. 84-14, May 21, 1984, 39 pages Submitted to SIAM J. Sci. Statist. Comput. The work of Young [1971] showed that the Red/Black ordering and the natural rowwise ordering of matrices with Property A, such as those arising from the 5-point discretization of Poisson's equation, lead to SOR iteration matrices with identical eigenvalues. With the advent of parallel computers, multi-color point SOR schemes have been proposed for more complicated stencils on 2-dimensional rectangular grids, see Adams and Ortega [1982] for example, but to our knowledge, no theory has been provided for the rate of convergence of these methods relative to that of the natural rowwise scheme. New results show that certain matrices may be reordered so the resulting multi-color SOR matrix has the same eigenvalues as that for the original ordering. In addition, for a wide range of stencils, we show how to choose multi-color orderings so the multi-color SOR matrices have the same eigenvalues as the natural rowwise SOR matrix. The strategy for obtaining these results is based on "data flow" concepts and can be used to reach Young's conclusions above for the 5-point stencil. The importance of these results is threefold. Firstly, a constructive and easy means of finding these multi-colorings is a direct consequence of the theory; secondly, multi-color SOR methods can be found that have the same rate of convergence as the natural rowwise SOR method for a wide range of stencils used to discretize partial differential equations; and thirdly, these multi- color SOR methods can be efficiently implemented on a wide class of parallel computers. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [11] Gatski, Thomas B. and Chester E. Grosch A numerical study of the two- and three-dimensional unsteady Navier-Stokes equations in velocity-vorticity variables using compact difference schemes ICASE Report No. 84-15, May 22, 1984, 20 pages Proc. of 9th ICNMED, Saclay, France, June 23-29, 1984 LECTURE NOTES IN PHYSICS Springer-Verlag. A compact finite-difference approximation to the unsteady Navier-Stokes equations in velocity-vorticity variables is used to numerically simulate a number of flows. These include two-dimensional laminar flow of a vortex evolving over a flat plate with an embedded cavity, the unsteady flow over an elliptic cylinder, and aspects of the transient dynamics of the flow over a rearward facing step. The methodology required to extend the two-dimensional formulation to three-dimensions is presented. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [12] Hafez, Mohamed, M. and Phillips, Timothy N. A modified least squares formulation for a system of first-order equations ICASE Report No. 84-16, May 29, 1984, 21 pages Submitted to IMACS J. Numer. Math. Second-order equations in terms of auxiliary variables similar to potential and stream functions are obtained by applying a weighted least squares formulation to a first-order system. The additional boundary conditions which are necessary to solve the higher order equations are determined and numerical results are presented for the Cauchy-Riemann equations. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [13] Hall, Phillip The Gortler vortex instability mechanism in three-dimensional boundary layers ICASE Report No. 84-17, June 6, 1984, 37 pages. Submitted to the Royal Society of London. It is well known that the two-dimensional boundary layer on a concave wall is centrifugally unstable with respect to vortices aligned with the basic flow for sufficiently high values of the Gortler number. However, in most situations of practical interest the basic flow is three-dimensional and previous theoretical investigations do not apply. In this paper the linear stability of the flow over an infinitely long swept wall of variable curvature is considered. If there is no pressure gradient in the boundary layer it is shown that the instability problem can always be related to an equivalent two- dimensional calculation. However, in general, this is not the case and even for small values of the crossflow velocity field dramatic differences between the two and three-dimensional problems emerge. In particular, it is shown that when the relative size of the crossflow and chordwise flow is (R**(-l/2)) where R is the Reynolds number of the flow, the most unstable mode is time- dependent. When the size of the crossflow is further increased, the vortices in the neutral location have their axes locally perpendicular to the vortex lines of the basic flow. In this regime the eigenfunctions associated with the instability become essentially 'centre modes' of the Orr-Sommerfeld equation destabilized by centrifugal effects. The critical Gortler number for such modes can be predicted by a large wavenumber asymptotic analysis; the results suggest that for order unity values of the ratio of the crossflow and chordwise velocity fields, the Gortler instability mechanism is almost certainly not operational. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [14] Gozani, J., Nachshon, A., and Turkel, E. Conjugate gradient coupled with multigrid for an indefinite problem ICASE Report No. 84-18, June 7, 1984, 11 pages Submitted to Proc. of 3rd IMACS International Symposium on Computer Methods for Partial Differential Equations Lehigh University, Bethlehem, PA, June 20-22, 1984 An iterative algorithm for the Helmholtz equation is presented. This scheme is based on the preconditioned conjugate gradient method for the normal equations. The preconditioning is one cycle of a multigrid method for the discrete Laplacian. The smoothing algorithm is red-black Gauss-Seidel and is constructed so it is a symmetric operator. The total number of iterations needed by the algorithm is independent of h. By varying the number of grids, the number of iterations depends only weakly on k when k**3 h**2 is constant. Comparisons with a SSOR preconditioner are presented. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [15] Malik, M. R., Zang, T. A., and Hussaini, M. Y. A spectral collocation method for the Navier-Stokes equations ICASE Report No. 84-19, June 8, 1984, 48 pages Submitted to J. Comput. Phys. A Fourier-Chebyshev spectral method for the incompressible Navier-Stokes equations is described. It is applicable to a variety of problems including some with fluid properties which vary strongly both in the normal direction and in time. In this fully spectral algorithm, a preconditioned iterative technique is used for solving the implicit equations arising from semi- implicit treatment of pressure, mean advection and vertical diffusion terms. The algorithm is tested by applying it to hydrodynamic stability problems in channel flow and in external boundary layers with both constant and variable viscosity. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [16] Davis, Stephen F. TVD finite difference schemes and artificial viscosity. ICASE Report No. 84-20, June 11, 1984, 35 pages Submitted to SIAM J. Sci. Statis. Comput. In this paper we show that the total variation diminishing (TVD) finite difference scheme which was analysed by Sweby can be interpreted as a Lax- Wendroff scheme plus an upwind weighted artificial dissipation term. We then show that if we choose a particular flux limiter and remove the requirement for upwind weighting, we obtain an artificial dissipation term which is based on the theory of TVD schemes, which does not contain any problem dependent parameters and which can be added to existing MacCormack method codes. Finally, we conduct numerical experiments to examine the performance of this new method. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [17] Canuto, C., Hariharan, S. I., and Lustman, L. Spectral methods for exterior elliptic problems ICASE Report No. 84-21, June 12, 1984, 33 pages. Submitted to Numer. Math. This paper deals with spectral approximations for exterior elliptic problems in two dimensions. As in the conventional finite difference or finite element methods, it is found that the accuracy of the numerical solutions is limited by the order of the numerical farfield conditions. We introduce a spectral boundary treatment at infinity, which is compatible with the "infinite order" interior spectral scheme. Computational results are presented to demonstrate the spectral accuracy attainable. Although we deal with a simple Laplace problem throughout the paper, our analysis covers more complex and general cases. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [18] Graif, E. and K. Kunisch Parameter estimation for the Euler-Bernoulli-beam ICASE Report No. 84-22, June 20, 1984, 43 pages. Submitted to IEEE Trans. Auto. Control. An approximation involving cubic spline functions for parameter estimation problems in the Euler-Bernoulli-beam equation (phrased as an optimization problem with respect to the parameters) is described and convergence is proved. The resulting algorithm was implemented and several of the test examples are documented. It is observed that the use of penalty terms in the cost functional can improve the rate of convergence. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [19] Rosen, I. Gary A numerical scheme for the identification of hybrid systems describing the vibration of flexible beams with tip bodies ICASE Report No. 84-23, June 21, 1984, 36 pages. Submitted to J. Math. Anal. Appl. A cubic spline based Galerkin-like method is developed for the identification of a class of hybrid systems which describe the transverse vibration of flexible beams with attached tip bodies. The identification problem is formulated as a least squares fit to data subject to the system dynamics given by a coupled system of ordinary and partial differential equations recast as an abstract evolution equation (AEE) in an appropriate infinite dimensional Hilbert space. Projecting the AEE into spline-based subspaces leads naturally to a sequence of approximating finite diemnsional identification problems. The solutions to these problems are shown to exist, are relatively easily computed, and are shown to, in some sense, converge to solutions to the original identification problem. Numerical results for a variety of examples are discussed. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [20] Banks, H. Thomas and Katherine A. Murphy Estimation of coefficients and boundary parameters in hyperbolic systems ICASE Report No. 84-24, June 21, 1984, 52 pages. Submitted to SIAM J. Control Optim. We consider semi-discrete Galerkin approximation schemes in connection with inverse problems for the estimation of spatially varying coefficients and boundary condition parameters in second order hyperbolic systems typical of those arising in 1-D surface seismic problems. Spline based algorithms are proposed for which theoretical convergence results along with a representative sample of numerical findings are given. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [21] Vardi, A. A new minmax problem ICASE Report No. 84-25, June 26, 1984, 35 pages. Submitted to J. Optim., Theory and Appl. For the classical minimax problem we design a new active set strategy in which there are three types of functions: active, semi-active, and non-active. This technique will help in preventing zigzagging which often occurs when an active set strategy is used. Some of the inequality constraints are handled with slack variables. Also a trust region strategy is used in which at each iteration there is a sphere around the current point in which the local approximation of the function is trusted. The algorithm suggested in the paper was implemented into a successful computer program. Numerical results are provided. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [22] Banks, H. T. and I. G. Rosen Approximation techniques for parameter estimation and feedback control for distributed models of large flexible structures ICASE Report No. 84-26, June 22, 1984, 19 pages. Submitted to Workshop on Identification and Control of Flexible Space Structures, San Diego, CA, June 4-6, 1984, G. Rodgriguez (ed.). We discuss approximation ideas that can be used in parameter estimation and feedback control for Euler-Bernoulli models of elastic systems. Focusing on parameter estimation problems, we outline how one can obtain convergence results for cubic spline-based schemes for hybrid models involving an elastic cantilevered beam with tip mass and base acceleration. Sample numerical findings are also presented. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [23] Turkel, E. and Bram van Leer Flux vector splitting and Runge-Kutta methods for the Euler equations ICASE Report No. 84-27, June 23, 2984, 9 pages. To appear in Proc. of the 9th Intl. Conference on Numerical Methods in Fluid Dynamics, Saclay, France, June 25-29, 1984. Runge-Kutta schemes have been used as a method of solving the Euler equations exterior to an airfoil. In the past this has been coupled with central differences and an artificial viscosity in space. In this study we couple the Runge-Kutta time stepping scheme with an upwinded space approxima tion based on flux-vector splitting. Several acceleration techniques are also considered including a local time step, residual smoothing and multigrid. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [24] Turkel, Eli Fast solutions to the steady state compressible and incompressible fluid dynamics equations ICASE Report No. 84-28, June 23, 1984, 8 pages. To appear in Proc. of the 9th Intl. Conference on Numerical Methods in Fluid Dynamics Saclay, France, June 25-29, 1984. It is well known that for low speed flows the use of the compressible fluid dynamic equations is inefficient. The use of an explicit scheme requires delta-t to be bounded by 1/c. However, the physical parameters change over time scales of order 1/u which is much larger. Hence, it is not appropriate to use explicit schemes for very subsonic flows. Implicit schemes are hard to vectorize and frequently do not converge quickly for very subsonic flows. We shall demonstrate that if one is only interested in the steady state then a minor change to an existing code can greatly increase the efficiency of an explicit method. Even when using an implicit method the proposed changes increase the efficiency of the scheme. We shall first consider the Euler equations for low speed flows and then incompressible flows. We then indicate how to generalize the method to include viscous effects. We also show how to accelerate supersonic flow by essentially decoupling the equations. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [25] Gottlieb, D. Spectral methods for compressible flow problems ICASE Report No. 84-29, June 21, 1984, 20 pages. Proc. 9th Intl. Conference on Numerical Methods in Fluid Dynamics Saclay, France, June 25-29, 1984. In this article we review recent results concerning numerical simulation of shock waves using spectral methods. We discuss shock fitting techniques as well as shock capturing techniques with finite difference artificial viscosity. We also discuss the notion of the information contained in the numerical results obtained by spectral methods and show how this information can be recovered. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [26] Ipsen, Ilse C. F. Singular value decomposition with systolic arrays ICASE Report No. 84-30, July 9, 1984, 22 pages. Proc. of the Society of Photo-Optical Engineers San Diego, CA, August 19-24, 1984. Systolic arrays for determining the Singular Value Decomposition of a matrix A of bandwidth w are presented. After A has been reduced to bidiagonal form B by means of Givens plane rotations, the differential equations. Part II: The linear quadratic control problem. ICASE Report No. 84-31, July 6, 1984, 65 pages. Submitted to SIAM J. Control Optim. The numerical scheme based on the Legendre-tau approximation is proposed to approximate the feedback solution to the linear quadratic optimal control problem for hereditary differential systems. The convergence property is established using Trotter ideas. The method yields very good approximations at low orders and provides an approximation technique for computing closed- loop eigenvalues of the feedback system. A comparison with existing methods (based on "averaging" and "spline" approximations) is made. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [27] Turkel, Eli Acceleration to a steady state for the Euler equations ICASE Report No. 84-32, July 7, 1984, 45 pages. Proc. INRIA Euler Workshop (SIAM) Rocquencourt, France, December 7-9, 1983. A multi-stage Runge-Kutta method is analyzed for solving the Euler equations exterior to an airfoil. Highly subsonic, transonic and supersonic flows are evaluated. Various techniques for accelerating the convergence to a steady state are introduced and analyzed. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [28] Bokhari, Shahid H. and A. D. Raza Augmenting computer networks ICASE Report No. 84-33, July 8, 1984, 16 pages. Proc. 1984 Intl. Conference on Parallel Processing Bellaire, MI, August 21-24, 1984. Three methods of augmenting computer networks by adding at most one link per processor are discussed. 1. A tree of N nodes may be augmented such that the resulting graph has diameter no greater than 4 log2((N+2)/3)-2. This 0(N**3) algorithm can be applied to any spanning tree of a connected graph to reduce the diameter of that graph to 0(log N). 2. Given a binary tree T and a chain C of N nodes each, C may be augmented to produce C' so that T is a subgraph of C'. This algorithm is 0(N) and may be used to produce augmented chains or rings that have diameter no greater than 2 log2 (N+2)/3) and are planar. 3. Any rectangular two-dimensional 4 (8) nearest neighbor array of size N = 2**k may be augmented so that it can emulate a single step shuffle- exchange network of size N/2 in 3(t) time steps. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [29] Reed, Daniel A. and Merrell L. Patrick A model of asynchronous iterative algorithms for solving large, sparse, linear systems ICASE Report No. 84-34, July 31, 1984, 27 pages. Proceedings of the 1984 International Conference on Parallel Processing August 21-24, 1984, Bellaire, MI. Solving large, sparse, linear systems of equations is one of the fundamental problems in large scale scientific and engineering computation. A model of a general class of asynchronous, iterative solution methods for linear systems is developed. In the model, the system is solved by creating several cooperating tasks that each compute a portion of the solution vector. This model is then analyzed to determine the expected intertask data transfer and task computational complexity as functions of the number of tasks. Based on the analysis, recommendations for task partitioning are made. These recommendations are a function of the sparseness of the linear system, its structure (i.e., randomly sparse or banded), and dimension. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [30] Reed, Daniel A. and Merrell L. Patrick Parallel, iterative solution of sparse linear systems: Models and architectures ICASE Report No. 84-35, August 1, 1984, 45 pages. Submitted to Parallel Computing. Solving large, sparse, linear systems of equations is a fundamental problem in large scale scientific and engineering computation. A model of a general class of asychronous, iterative solution methods for linear systems is developed. In the model, the system is solved by creating several cooperating tasks that each compute a portion of the solution vector. A data transfer model predicting both the probability that data must be transferred between two tasks and the amount of data to be transferred is presented. This model is used to derive an execution time model for predicting parallel execution time and an optimal number of tasks given the dimension and sparsity of the coefficient matrix and the costs of computation, synchronization, and communication. The suitability of different parallel architectures for solving randomly sparse linear systems is discussed. Based on the complexity of task scheduling, one parallel architecture, based on a broadcast bus, is presented and analyzed. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [31] Tadmor, Eitan The well-posedness of the Kuramoto-Sivashinsky equation. ICASE Report No. 84-36, August 7, 1984, 22 pages. Submitted to SIAM J. Appl. Math. The Kuramoto-Sivashinsky equation arises in a variety of applications, among which are modeling reaction-diffusion systems, flame-propagation and viscous flow problems. It is considered here, as a prototype to the larger class of generalized Burgers equations: those consist of quadratic nonlinearity and arbitrary linear parabolic part. We show that such equations are well-posed, thus admitting a unique smooth solution, continuously dependent on its initial data. As an attractive alternative to standard energymethods, existence and stability are derived in this case, by "patching" in the large short time solutions without "loss of derivatives". Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [32] Lustman, Liviu The time evolution of spectral discretizations of hyperbolic systems ICASE Report No. 84-37, August 7, 1984, 12 pages. Submitted to SIAM J. Numer. Anal. A Chebyshev collocation spectral method, applied to hyperbolic systems is considered, particularly for those initial boundary value problems which possess only solutions tending to zero at large times. It is shown that the numerical solutions of the system will also vanish at infinity, if and only if, the numerical solution of a scalar equation of the same type does. This result is then generalized for other spectral approximations. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [33] Bayliss, A., C. I. Goldstein, and E. Turkel On accuracy conditions for the numerical computation of waves ICASE Report No. 84-38, August 9, 1984, 16 pages. Submitted to J. Comput. Phys. The Helmholtz equation with a variable index of refraction, and a suitable radiation condition at infinity serves as a model for a wide variety of wave propagation problems. Such problems can be solved numerically by first truncating the given unbounded domain and imposing a suitable outgoing radiation condition on an artificial boundary and then solving the resulting problem on the bounded domain by direct discretization (for example, using a finite element method). In practical applications, the mesh size h and the wave number K, are not independent but are constrained by the accuracy of the desired computation. It will be shown that the number of points per wavelength, measured by l/Kh, is not sufficient to determine the accuracy of a given discretization. For example, the quantity K**3 h**2 is shown to determine the accuracy in the L-two norm for a second-order discretization method applied to several propagation models. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [34] Hariharan, S. I. Numerical solutions of acoustic wave propagation problems using Euler computations ICASE Report No. 84-39, August 10, 1984, 29 pages. Proc. of the AIAA 9th Aeroacoustics Conference Williamsburg, VA, August 15-17, 1984. This paper reports solution procedures for problems arising from the study of engine inlet wave propagation. The first problem is the study of sound waves radiated from cylindrical inlets. The second one is a quasi-one- dimensional problem to study the effect of nonlinearities and the third one is the study of nonlinearities in two dimensions. In all three problems Euler computations are done with a fourth-order explicit scheme. For the first problem results are shown in agreement with experimental data and for the second problem comparisons are made with an existing asymptotic theory. The third problem is part of an ongoing work and preliminary results are presented for this case. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [35] Tadmor, Eitan The exponential accuracy of Fourier and Tchebyshev differencing methods ICASE Report No. 84-40, August 20, 1984, 24 pages. Submitted to SIAM J. Numer. Anal. It is shown that when differencing analytic functions using the pseudospectral Fourier or Tchebyshev methods, the error committed decays to zero at an exponential rate. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [36] Gannon, Dennis and John Van Rosendale On the impact of communication complexity in the design of parallel numerical algorithms ICASE Report No. 84-41, August 22, 1984, 37 pages. To appear in IEEE Trans. on Computers. This paper describes two models of the cost of data movement in parallel numerical algorithms. One model is a generalization of an approach due to Hockney, and is suitable for shared memory multiprocessors where each processor has vector capabilities. The other model is applicable to highly parallel nonshared memory MIMD systems. In the second model, algorithm performance is characterized in terms of the communication network design. Techniques used in VLSI complexity theory are also brought in, and algorithm independent upper bounds on system performance are derived for several problems that are important to scientific computation. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [37] Ito, Kazufumi Legendre-Tau approximation for functional differential equations, Part III: Eigenvalue approximations and uniform stability. ICASE Report No 84-42, August 24, 1984, 34 pages. Proc. Second International Conference on Control Theory for Distributed Parameter Systems and Applications Vorau, Austria, 1984. The stability and convergence properties of the Legendre-tau approximation for hereditary differential systems are analyzed. We derive a characteristic equation for the eigenvalues of the resulting approximate system. As a result of this derivation we are able to establish that the uniform exponential stability of the solution semigroup is preserved under approximation. It is the key to obtaining the convergence of approximate solutions of the algebraic Riccati equation in trace norm. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [38] Berger, M. J. On conservation at grid interfaces ICASE Report No. 84-43, September 7, 1984, 25 pages. Submitted to J. Comput. Phys. This paper considers the solution of hyperbolic systems of conservation laws on discontinuous grids. In particular, we consider what happens to conservation at grid interfaces. A procedure is presented to derive conservative difference approximations at the grid interfaces for two- dimensional grids which overlap in an arbitrary configuration. The same procedures are applied to compute interface formulas for grids which are refined in space and/or time, and for continuous grids where a switch in the scheme causes the discontinuity. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [39] Osher, S. and S. Chakravarthy Very high order accurate TVD schemes ICASE Report No. 84-44, September 10, 1984, 61 pages. Submitted to SIAM J. Numer. Anal. A systematic procedure for constructing semi-discrete families of 2m-1 order accurate, 2m order dissipative, variation diminishing, 2m+1 point bandwidth, conservation form approximations to scalar conservation laws is presented. Here m is any integer between 2 and 8. Simple first-order forward time discretization, used together with any of these approximations to the space derivatives, also results in a fully discrete, variation diminishing algorithm. These schemes all use simple flux limiters, without which each of these fully discrete algorithms is even linearly unstable. Extensions to systems, using a nonlinear field-by-field decomposition are presented, and shown to have many of the same properties as in the scalar case. For linear systems, these nonlinear approximations are variation diminishing, and hence convergent. A new and general criterion for approximations to be variation diminishing is also given. Finally, numerical experiments using some of these algorithms are presented. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [40] Elcrat, A. R. and L. N. Trefethen Classical free-streamline flow over a polygonal obstacle ICASE Report No. 84-45, September 10, 1984, 61 pages. Submitted to J. Comput. Appl. Math. In classical Kirchhoff flow, an ideal incompressible fluid flows past an obstacle and around a motionless wake bounded by free streamlines. Since 1869 it has been known that in principle, the two-dimensional Kirchhoff flow over a polygonal obstacle can be determined by constructing a conformal map onto a polygon in the log-hodograph plane. In practice, however, this idea has rarely been put to use except for very simple obstacles, because the conformal mapping problem has been too difficult. This paper presents a practical method for computing flows over arbitrary polygonal obstacles to high accuracy in a few seconds of computer time. We achieve this high speed and flexibility by working with a modified Schwarz-Christoffel integral that maps onto the flow region directly rather than onto the log-hodograph polygon. This integral and its associated parameter problem are treated numerically by methods developed earlier by Trefethen for standard Schwarz-Christoffel maps. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [41] Krause, Egon Review of some vortex relations ICASE Report No. 84-46, September 12, 1984, 8 pages. Submitted to Computer and Fluids. The evaluation of the circulation from numerical solutions of the momentum and energy equations is discussed for incompressible and compressible flows. It is shown how artificial damping directly influences the time rate of change of the circulation. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [42] Tal-Ezer, Hillel A pseudospectral Legendre method for hyperbolic equations with an improved stability condition ICASE Report No. 84-48, September 14, 1984, 46 pages. Submitted to J. Comput. Phys. A new pseudospectral method is introduced for solving hyperbolic partial differential equations. This method uses different grid points than previously used pseudospectral methods: in fact the grid points are related to the zeroes of the Legendre polynomials. The main advantage of this method is that the allowable time step is proportional to the inverse of the number of grid points 1/N rather than to 1/N**2 (as in the case of other pseudospectral methods applied to mixed initial boundary value problems). A highly accurate time discretization suitable for these spectral methods is discussed. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [43] Bayliss, A., C. I. Goldstein, and E. Turkel The numerical solution of the Helmholtz equation for wave propagation problems in underwater acoustics ICASE Report No. 84-49, September 17, 1983, 29 pages. Submitted to J. Comput. Math. Appl. The Helmholtz Equation with a variable index of refraction, and a suitable radiation condition at infinity serves as a model for a wide variety of wave propagation problems. A numerical algorithm has been developed and a computer code implemented that can effectively solve this equation in the intermediate frequency range. The equation is discretized using the finite element method, thus allowing for the modeling of complicated geometries (including interfaces) and complicated boundary conditions. A global radiation boundary condition is imposed at the far field boundary that is exact for an arbitrary number of propagating modes. The resulting large, non-selfadjoint system of linear equations with indefinite symmetric part is solved using the preconditioned conjugate gradient method applied to the normal equations. A new preconditioner is developed based on the multigrid method. This preconditioner is vectorizable and is extremely effective over a wide range of frequencies provided the number of grid levels is reduced for large frequencies. A heuristic argument is given that indicates the superior convergence properties of this preconditioner. The relevant limit to analyze convergence is for K increasing and a fixed prescribed accuracy level. The efficiency and robustness of the numerical algorithm is confirmed for large acoustic models, including interfaces with strong velocity contrasts. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [44] Burns, John A., Ito, Kazufumi, Powers, Robert K. Chandrasekhar equations and computational algorithms for distributed parameter systems ICASE Report No. 84-50, September 20, 1984, 23 pages. Proc. of 23rd Conference on Decision and Control, December 12-14, 1984, Las Vegas, NV. In this paper we consider the Chandrasekhar equations arising in optimal control problems for linear distributed parameter systems. The equations are derived via approximation theory. This approach is used to obtain existence, uniqueness and strong differentiability of the solutions and provides the basis for a convergent computation scheme for approximating feedback gain operators. A numerical example is presented to illustrate these ideas. Submitted by: emily@icase.csnet Obtainable from: Ms. Barbara Kraft (as above) --------------- [45] Petter E. Bjorstad and Olof B. Widlund Iterative Methods for the Solution of Elliptic Problems on Regions Partitioned into Substructures Finite element problems can often naturally be divided into subproblems which correspond to subregions into which the region has been partitioned or from which it was originally assembled. A class of iterative methods is discussed in which these subproblems are solved by direct methods,while the interaction across the curves or surfaces which divide the region is handled by a conju- gate gradient method. A mathematical framework for this workis provided by regularity theory for elliptic finite element problems and by block Gauss elimination. A full development of the theory,which shows that certain of these methods are optimal, is given for Lagrangian finite element approxi- mations of second order linear elliptic problems in the plane.Results from numerical experiments are also reported. Submitted by: Olaf Widlund Obtainable from: Olaf Widlund Courant Institute 251 Mercer St. N.Y., N.Y. 10012. (WIDLUND@NYU-ACF1.ARPA) --------------- [46] Ron S. Dembo and Ulrich Tulowitzki Local Convergence Analysis for Successive Inexact Quadratic Programming Methods SOM Working Paper Series B #78 Yale School of Organization and Management October, 1984 Successive quadratic programming (SQP) methods for the solution of constrained nonlinear optimization problems are attractive because they converge rapidly from any sufficiently good initial solution. However, solving the quadratic subproblems at each stage can be expensive, particularly if the number of unknowns is large and may be not justified when far away from an optimal solution. Therefore, we consider a class of successive inexact quadratic programming (SIQP) methods that solve the subproblems only approximately. We characterize the rate of convergence in terms of the relative error incurred in the subproblems in a manner that does not assume that the set of active constraints remains constant. This leads to practical termination criteria for truncating various iterative algorithms applied to the quadratic subproblems. Submitted by: Ron Dembo Obtainable from: Ron Dembo Yale University School of Organization and Management P.O. Box 1A New Haven, Conn. 06520 --------------- [47] Ron S. Dembo and Ulrich Tulowitzki Sequential Truncated Quadratic Programming Methods To appear in the Proceedings of the SIAM Conference on Numerical Optimization, Boulder, Colo., June, 1984 In a sequential quadratic programming (SQP) method for nonlinear optimization it is possible to truncate an iterative procedure, applied to each quadratic subproblem, prior to reaching a solution, without affecting the asymptotic rate of convergence. We describe some practical termination criteria which ensure that an SQP method, executed without exact solution of the QP subproblem, inherits the rate of convergence of its (exact) SQP counterpart. Numerical experiments on linearly-constrained problems of up to 2230 variables indicate that substantial savings can be achieved by inexact solution of the QP subproblems with little or no adverse effect on the convergence rate. Submitted by: Ron Dembo Obtainable from: Ron Dembo Yale University School of Organization and Management P.O. Box 1A New Haven, Conn. 06520 --------------- From rhbartels%watcgl%waterloo.csnet@csnet-relay.arpa Tue Jul 16 14:50:12 1985 Received: from csnet-relay (csnet-relay.arpa.ARPA) by anl-mcs.ARPA ; Tue, 16 Jul 85 14:46:21 cdt Received: from waterloo by csnet-relay.csnet id aa01168; 16 Jul 85 14:32 EDT Date: Sun, 14 Jul 85 20:31:32 edt From: Richard Bartels <rhbartels%watcgl%waterloo.csnet@csnet-relay.arpa> Message-Id: <8507150031.AA20620@watcgl.UUCP> To: dongarra@anl-mcs Status: R ================================================== NUMERICAL ANALYSIS TITLES -- Volume 3, Number 1 April 1, 1985 (No foolin') ================================================== Note: Responding to a request by John Lewis, and for the benefit of the SIGNUM Newsletter, these titles and abstracts are now ordered alphabetically by insti- tution. ##### ARGONNE ##### All Argonne reports listed here may be obtained by writing to the person listed as submitting the report at the following address: Argonne National Laboratory Mathematics and Computer Science Division 9700 S. Cass Ave. Argonne, IL 60439 ################### [1] Proceedings for the Argonne Workshop on Programming the Next Generation of SuperComputers, October 1984 Report ANL/MCS-TM-34 Argonne National Lab On February 27-28, 1984, an Argonne National Laboratory workshop titled "Programming the Next Generation of Supercomputers" was held at the Amfac Hotel in Albuquerque, New Mexico. The purpose of the workshop was to gather together researchers from the DOE laboratories who were using multiprocessors and to have these investigators relate their experiences with the programming languages supplied with these processors. The goal was to identify those areas where additional language primitives would aid in the development of portable software. Upon the identification of such needs, extensions to the Fortran language would then be proposed to such groups as the DOE Language Working Group and the X3J3 Fortran committee for further study. The topics presented at the workshop included reviews of existing languages that address concurrency and parallel computation, presentations on portable programming methods for parallel processing, users' views of what is needed to program multiprocessors, and the supplier's views of what is needed to drive their processors. This report collects the presentations given at this workshop and records the organized discussion held after each presentation and each session. In most cases, the materials reproduced here are the slides given by each author at the workshop. In some cases, papers were prepared by the authors and included in this report. A summary of the conclusions reached by the workshop is also included. Submitted by: Brian T. Smith <smith@anl-mcs.arpa> --------------- [2] Iain S. Duff Parallel implementation of multifrontal schemes ANL/MCS-TM-49, March, 1985 We consider the direct solution of large sparse sets of linear equations in an MIMD environment. We base our algorithm on a multifrontal approach and show how to distribute tasks among processors according to an elimination tree which can be automat- ically generated from any pivot strategy. We study organiza- tional aspects of the scheme for shared memory multiprocessor configurations and consider implications for multiprocessors with local memories and a communication network. Submitted by: Iain S. Duff <duff@anl-mcs.arpa> --------------- [3] I.S. Duff, A.M. Erisman, C.W. Gear, and J.K. Reid Some remarks on inverses of sparse matrices ANL/MCS-TM-51, March, 1985 We discuss some of the properties of the structure of the inverse of a sparse matrix. In particular, we show that the structural inverse of an irreducible sparse matrix is always full. We first show that Gaussian elimination always yields a full structural inverse and then show the result directly. We believe that our lemmas concerning Gaussian elimination are of interest in their own right. Submitted by: Iain S. Duff <duff@anl-mcs.arpa> --------------- [4] I.S. Duff and C.W. Gear Computing the structural index ANL/MCS-TM-50, March, 1985 The index of many differential/algebraic equations is deter- mined by the structure of the system, that is, by the pattern of nonzero entries in the Jacobians. It is important to be able to determine if the index does not exceed two because such problems can be solved numerically. In this paper we present an algorithm that determines if the index is one, two, or greater, based only on the structure. The algorithm can be exponential in its execu- tion time: we do not know whether it is possible to get an asymp- totically faster algorithm. However, in many practical problems, this algorithm will execute in polynomial time. Submitted by: Iain S. Duff <duff@anl-mcs.arpa> --------------- [5] Jack J. Dongarra, Linda Kaufman and, Sven Hammarling, Squeezing the Most out of Eigenvalue Solvers on High-Performance Computers ANL-MCS/TM 46, January, 1985 This paper describes modifications to many of the standard algorithms used in computing eigenvalues and eigenvectors of matrices. These modifications can dramatically increase the performance of the underlying software on high performance computers without resorting to assembler language, without significantly influencing the floating point operation count, and without affecting the roundoff error properties of the algorithms. The techniques are applied to a wide variety of algorithms and are beneficial in various architectural settings. Submitted by: Jack J. Dongarra <dongarra@anl-mcs.apa> --------------- [6] Jack J. Dongarra, Jeremy Du Croz, Sven Hammarling, and Richard J. Hanson A Proposal for an Extended Set of Basic Linear Algebra Subprograms ANL-MCS/TM 41, December, 1984 This paper describes an extension to the set of Basic Linear Algebra Subprograms. The extensions proposed are targeted at matrix vector operations which should provide for more efficient and portable implementations of algorithms for high performance computers. Submitted by: Jack J. Dongarra <dongarra@anl-mcs.apa> --------------- [7] J.J. Dongarra, A.H. Sameh, and D.C. Sorensen Implementation of Some Concurrent Algorithms for Matrix Factorization ANL-MCS/TM 25, October, 1984 Three parallel algorithms for computing the QR-factorization of a matrix are presented. The discussion is primarily concerned with implementation of these algorithms on a computer that supports tightly coupled parallel processes sharing a large common memory. The particular experiments were conducted on the Denelcor HEP computer. The three algorithms are a Householder method based upon high-level modules, a Windowed Householder method that avoids fork-join synchronization, and a Pipelined Givens method that is a variant of the data-flow type algorithms offering large enough granularity to mask synchronization costs. Computational results indicate that the Pipelined Givens method is preferred and that this is primarily due to the number of array references required by the various algorithms. Submitted by: Jack J. Dongarra <dongarra@anl-mcs.apa> --------------- [8] J.J. Dongarra and D.C. Sorensen A Parallel Linear Algebra Library for the Denelcor HEP ANL-MCS/TM 33, October, 1984 This paper discusses the implementation of a library of algorithms for problems in linear algebra on the Denelcor HEP. The package includes some of the most heavily used subroutines from LINPACK, that is, solution of linear systems based on LU, Cholesky, and QR factorizations as well as the appropriate triangular solvers. The concept followed is to code these routines in terms of high-level modules, which provides a vehicle to achieve both transportability and efficiency across a wide range of architectures. We discuss this concept in the context of a numerical linear algebra software library which is adaptable to highly parallel computing systems. However, the concept is expected to applicable to other libraries as well. Submitted by: Jack J. Dongarra <dongarra@anl-mcs.apa> --------------- [9] M. T. Heath and D.C. Sorensen A Pipelined Givens Method for Computing the QR Factorization of a Sparse Matrix ANL-MCS/TM 47, February, 1985 This paper discusses an extension of the Pipelined Givens method for computing the QR factorization of a real m x n matrix to the case in which the matrix is sparse. When restricted to one process, the algorithm performs the same computation as the serial sparse Givens algorithm of George and Heath. Our implementation is compatible with the data structures used in SPARSPAK. The pipelined algorithm is well suited to parallel computers having globally shared memory and low overhead synchronization primitives, such as the Denelcor HEP, for which computational results are presented. We point out certain synchronization problems that arise in the adaptation to the sparse setting and discuss the effect on parallel speedup of accessing a serial data file. Submitted by: Jack J. Dongarra <dongarra@anl-mcs.apa> --------------- ##### BELL LABS ##### All Bell Laboratories reports were submitted by Eric Grosse <{ucbvax,ihnp4}!research!ehg> and they may be obtained by writing to him at: AT&T Bell Labs 2c471 Murray Hill, NJ 07974 or by requesting them of him through electronic mail. ##################### [10] Eric Grosse AUTOMATIC CHOICE OF COLORS FOR LEVEL PLOTS TR/NAMs 85-1 AT&T Bell Labs Color level plots are becoming popular in scientific computing but often require manually adjusting the color codes. A formula, fit by Friele, MacAdam, and Chickering to ellipses of measured color discrimination, provides the basis for a spectrum of equally spaced colors. The color choice problem is formulated as a nonlinear least squares optimization and solved numerically. A table can be precomputed to adjust the hue coordinate in the hexcone model of Alvy Ray Smith, just as compensation tables are used for gamma correction of CRTs. --------------- [11] William M. Coughran, Jr. and Eric Grosse THE GROVE EDITOR TR/NAMs 83-3 We describe a structure editor, grove, which manipulates groves of trees rather than text files. Various commands for searching and modifying nodes and subtrees are described. We also document some routines for performing basic operations on trees, which facilitate the construction of program- transformation systems. --------------- [12] William M. Coughran, Jr. and Eric Grosse THE PINE PROGRAMMING LANGUAGE TR/NAMs 83-4 Pine is a language intended to bring to scientific computation convenient program transformation, array operations, and simplified calls to libraries. Pine programs can be translated into FORTRAN. The key features added to FORTRAN are array operations, heap-based dynamic memory allocation, more convenient input/output, generic functions, operators of the form +=, expression-language syntax, and control structures. The lexical form of the language, based on a screen editor of trees, is designed to simplify program transformations. This manual describes pine and its processor to those already familiar with the issues involved in compiling into FORTRAN. --------------- [13] William M. Coughran, Jr., Eric Grosse, and Donald J. Rose VARIATION DIMINISHING SPLINES IN SIMULATION TR/NAMs 84-3 Variation diminishing splines provide an effective tool for modeling active elements in circuit simulation. Using quadratic tensor product splines and maintaining uniform sampling at the boundary by linear extension of the data yields an algorithm that is smooth (unlike simple table lookup), shape preserving (unlike simple interpolation), and efficient (30 microseconds to evaluate on a Cray-1A). The rate of convergence to function and derivative values is O(h**2). --------------- [14] Brenda S. Baker, Eric Grosse, and Conor S. Rafferty NON-OBTUSE TRIANGULATION OF A POLYGON TR/NAMs 84-4 We show how to triangulate a polygon without using any obtuse triangles. Such triangulations can be used to discretize partial differential equations in a way that guarantees the resulting matrix is Stieltjes, a desirable property both for computation and for theoretical analysis. A simple divide-and-conquer approach would fail because adjacent subproblems cannot be solved independently, but this can be overcome by careful subdivision. Overlay a square grid on the polygon, preferably with the polygon vertices at grid points. Choose boundary cells so they can be triangulated without propagating irregular points to adjacent cells. The remaining interior is rectangular and easily triangulated. Small angles can also be avoided in these constructions. --------------- [15] Jack J. Dongarra and Eric Grosse DISTRIBUTION OF MATHEMATICAL SOFTWARE VIA ELECTRONIC MAIL TR/NAMs 85-2 A large collection of mathematical software is now available via electronic mail. Messages sent to "netlib@anl-mcs" (on the Arpanet/CSNET) or to "research!netlib" (on the UNIX\(tm network) wake up a server that distributes items from the collection. For example the one-line message, "send index", gets a library catalog by return mail. We describe how to use the service and some of the issues in its implementation. --------------- [15] Linda Kaufman, Jack Dongarra, and Sven Hammarling SQUEEZING THE MOST OUT OF EIGENVALUE SOLVERS ON HIGH-PERFORMANCE COMPUTERS TR/NAMs 84-7 This paper describes modifications to many of the standard algorithms used in computing eigenvalues and eigenvectors of matrices which decrease the number of vector-memory references. These modifications can dramatically increase the performance of the underlying software without resorting to assembly language and without significantly influencing the floating point operation count. --------------- ##### CORNELL ##### [16] Franklin T. Luk A Parallel Method for Computing the Generalized Singular Value Decomposition TR/EE-CEG-85-1 We describe a new parallel algorithm for computing the generalized singular value decomposition of two n-by-n matrices, one of which is nonsingular. Our procedure requires O(n) time and one triangular array of O( n**2 ) processors. Submitted by: Franklin T. Luk <luk%tesla@cornell.arpa> Obtainable from: Franklin T. Luk School of Electrical Engineering Phillips Hall Cornell University Ithaca, New York 14853 --------------- [17] Thomas F. Coleman, Burton S. Garbow, and Jorge J. More' Software for Estimating Sparse Hessian Matrices Argonne National Laboratory Cornell University MCS/TM-43 CSD/TR 85-660 January, 1985 The solution of a nonlinear optimization problem often requires an estimate of the Hessian matrix for a function f. In large scale problems the Hessian matrix is usually sparse, and then estimation by differences of gradients is attractive because the number of differences can be small compared to the dimension of the problem. In this paper we describe a set of subroutines whose purpose is to estimate the Hessian matrix with the least possible number of gradient evaluations. Submitted by: Thomas F. Coleman <coleman%gvax@cornell.arpa> Obtainable from: Thomas F. Coleman Cornell University Computer Science Dept. Ithaca, NY 14853 --------------- ##### NYU: COURANT INSTITUTE ##### [18] James Demmel On the Conditioning of Pole Assignment Computer Science Dept. Technical Report, Jan 85 An algorithm for the pole assignment problem has recently been given by Kautsky,Nichols, and Van Dooren. Given the n by n matrix A, the n by m full rank matrix B (m<n), and the n complex numbers z(i), the problem is to find F such that A+BF has the z(i) as eigenvalues. The quality of the solution is measured by the condition number of X, the matrix of eigenvectors of A+BF. In this note we show that a lower bound on cond(X) is given essentially by the reciprocal of the product of three terms: the condition number of B, the maximum condition number of any A- z(i)I, and the relative distance from the pair (A,B) to the nearest uncontrollable pair. This last term can be interpreted as follows: the lower bound on cond(X) becomes large as the relative distance from the problem to the nearest "infinitely ill-conditioned" problem becomes small. Submitted by: James Demmel <demmal@nyu-csd2.arpa> Obtainable from: James Demmel Courant Institute 251 Mercer Str. New York, NY 10012 --------------- ##### DELFT (NETHERLANDS) ##### [19] Henk A. van der Vorst The Performance of FORTRAN Implementations for Preconditioned Conjugate Gradients on Vector Computers We will consider in detail the performance of FORTRAN implemen- tations for the conjugate gradient algorithm, for the solution of large linear systems, on a number of well-known vectorcomputers: CRAY-1, CRAY X-MP and CYBER 205. Lowerbounds on the CPU-times, required for separate parts of the algorithm, are presented and these are compared to the actually observed CPU-times. It appears that these lowerbounds are reasonably sharp. Submitted by: Henk A. van der Vorst <decvax!mcvax!dutinfd!numan> Obtainable from: Henk A. van der Vorst Department of Mathematics University of Technology at Delft Julianalaan 132 Delft, the Netherlands --------------- ##### EMORY ##### [20] J.M. Borwein and H. Wolkowicz A Simple Constraint Qualification in Infinite Dimensional Programming A new, simple, constraint qualification for infinite dimensional programs with linear programming type constraints is used to derive the dual program. Applications include a proof of the explicit solution of the best interpolation problem of Micchelli, Smith, Swetits, and Ward. Submitted by: Henry Wolkowicz <henry@emory.csnet> Obtainable from: J.M. Borwein Dalhousie University Halifax, Nova Scotia B3H 4H8 Canada (or from Henry, but he didn't include his address). --------------- ##### IBM ##### [21] Jane Cullum and Ralph A. Willoughby A Practical Procedure for Computing Eigenvalues of Large Sparse Nonsymmetric Matrices IBM Research Report, RC 10988, 2/14/85 We propose a Lanczos procedure with no reorthogonalization for computing eigenvalues of very large nonsymmetric matrices. (This procedure can also be used to compute corresponding eigenvectors but that issue will be dealt with in a separate paper.) Such computations are for example, central to transient stability analyses of electrical power systems and for determining parameters for iterative schemes for the numerical solution of partial differential equations. Numerical results for several large matrices are presented to demonstrate the effectiveness of this procedure. Submitted by: Jane Cullum <cullumj%yktvmv@ibm-sj.csnet> Obtainable from: Jane Cullum and Ralph A. Willoughby IBM T.J. Watson Research Center, P.O. Box 218 Yorktown Heights, New York 10598 --------------- ##### ICASE ##### Barbara Kraft ICASE, M/S 132C NASA Langley Research Center Hampton, VA 23665 may be contacted for all ICASE reports, or a message may be sent to: emily@icase.csnet All of these reports were submitted by "Emily". (The list below has been abbreviated to the numerically-oriented reports. Only the original technical-report information is given, although some of these have appeared in proceedings and journals. Sorry. -RHB) ################## [22] Rosen, I. G. Approximation methods for inverse problems involving the vibration of beams with tip bodies. ICASE Report No. 84-51, September 26, 1984, 10 pages. Proc. 23rd IEEE Conference on Decisions and Control, Las Vegas, NV, Decmeber 12-14, 1984. We outline two cubic spline based approximation schemes for the estimation of structural parameters associated with the transverse vibration of flexible beams with tip appendages. The identification problem is formulated as a least squares fit to data subject to the system dynamics which are given by a hybrid system of coupled ordinary and partial differential equations. The first approximation scheme is based upon an abstract semigroup formulation of the state equation while a weak/variational form is the basis for the second. Cubic spline based subspaces together with a Rayleigh-Ritz- Galerkin approach was used to construct sequences of easily solved finite dimensional approximating identification problems. Convergence results are briefly discussed and a numerical example demonstrating the feasibility of the schemes and exhibiting their relative performance for purposes of comparison is provided. --------------- [23] Canuto, C. Boundary Conditions in Chebyshev and Legendre Methods ICASE Report No. 84-52, October 1, 1984, 38 pages We discuss two different ways of treating non-Dirichlet boundary conditions in Chebyshev and Legendre collocation methods for second order differential problems. An error analysis is provided. The effect of preconditioning the corresponding spectral operators by finite difference matrices is also investigated. --------------- [24] Roe, P. L. Generalized formulation of TVD Lax-Wendroff schemes ICASE Report No. 84-53, October 23, 1984, 14 pages The work of Davis which imports the concept of total-variation- diminution (TVD) into non-upwinded, Lax-Wendroff type schemes is reformulated in a way which is easier to analyze. The analysis reveals a class of TVD schemes not observed by Davis. --------------- [25] Bokhari, Shahid H. Shuffle-exchanges on augmented meshes ICASE Report No. 84-54, October 24, 1984, 12 pages Proc. 1985 Intl. Conference on Distributed Computer Systems Denver, CO. Prior research has shown how a mesh connected array of size N = 2**K an integer, can be augmented by adding at most one edge per node such that it can perform a shuffle-exchange of size N/2 in constant time. We now show how to perform a shuffle-exchange of size N on this augmented array in constant time. This is done by combining the available perfect shuffle of size N/2 with the existing nearest neighbor connections of the mesh. By carefully scheduling the different permutations that are composed in order to achieve the shuffle, the time required is reduced to 5 steps, which is shown to be optimal for this network. --------------- [26] Goodman, Jonathan B. and Randall J. LeVeque A geometric approach to high resolution TVD schemes ICASE Report No. 84-55, October 29, 1984, 35 pages We use a geometric approach, similar to Van Leer's MUSCL schemes, to construct a second-order accurate generalization of Godunov's method for solving scalar conservation laws. By making suitable approximations we obtain a scheme which is easy to implement and total variation diminishing. We also investigate the entropy condition from the standpoint of the spreading of rarefaction waves. For Godunov's method we obtain quantitative information on the rate of spreading which explain the kinks in rarefaction waves often observed at the sonic point. --------------- [27] Tin-Chee Wong, C. H. Liu, and James Geer Comparison of uniform perturbation solutions and numerical solutions for some potential flows past slender bodies ICASE Report No. 84-56, October 29, 1984, 32 pages Approximate solutions for potential flow past an axisymmetric slender body and past a thin airfoil are calculated using a uniform perturbation method and then compared with either the exact analytical solution or the solution obtained using a purely numerical method. The perturbation method is based upon a representation of the disturbance flow as the superposition of singularities distributed entirely within the body, while the numerical (panel) method is based upon a distribution of singularities on the surface of the body. It is found that the perturbation method provides very good results for small values of the slenderness ratio and for small angles of attack. Moreover, for comparable accuracy, the perturbation method is simpler to implement, requires less computer memory, and generally uses less computation time than the panel method. In particular, the uniform perturbation method yields good resolution near the regions of the leading and trailing edges where other methods fail or require special attention. --------------- [28] Salas, M. D., S. Abarbanel and D. Gottlieb Multiple steady states for characteristic initial value problems ICASE Report No. 84-57, November 1, 1984, 43 pages The time dependent, isentropic, quasi-one-dimensional equations of gas dynamics and other model equations are considered under the constraint of characteristic boundary conditions. Analysis of the time evolution shows how different initial data may lead to different steady states and how seemingly anamolous behavior of the solution may be resolved. Numerical experimentation using time consistent explicit algorithms verifies the conclusions of the analysis. The use of implicit schemes with very large time steps leads to erroneous results. --------------- [29] Rose, Milton E. A compact finite element method for elastic bodies ICASE Report No. 84-58, November 7, 1984, 38 pages A nonconforming finite element method is described for treating linear equilibrium problems, and a convergence proof showing second order accuracy is given. The close relationship to a related compact finite difference scheme due to Phillips and Rose is examined. A condensation technique is shown to preserve the compactness property and suggests an approach to a certain type of homogenization. --------------- [30] Lustman, Liviu A review of spectral methods ICASE Report No. 84-59, November 19, 1984, 25 pages This paper presents a brief outline of spectral methods for partial differential equations. The basic ideas, together with simple proofs are discussed. An application to potential transonic flow is also reviewed. --------------- [31] Hussaini, M. Y. and W. D. Lakin Existence and non-uniqueness of similarity solutions of a boundary layer problem ICASE Report No. 84-60, November 21, 1984, 17 pages This work considers a Blasius boundary value problem with inhomogeneous lower boundary conditions f(0) = 0 and f'(0) = - L with L strictly positive. The Crocco variable formulation of this problem has a key term which changes sign in the interval of interest. It is shown that solutions of the boundary value problem do not exist for values of L larger than a positive critical value LC. The existence of solutions is proved for 0 < L _ LC by considering an equivalent initial value problem. However, for 0 < L < LC, solutions of the boundary value problem are found to be nonunique. Physically, this non-uniqueness is related to multiple values of the skin friction. --------------- [32] Bayliss, Alvin, Maestrello, Lucio and Turkel, Eli On the interaction of a sound pulse with the shear layer of an axisymmetric jet, III. Nonlinear effects ICASE Report No. 84-61 November 23, 1984, 22 pages The fluctuating field of a jet excited by transient mass injection is simulated numerically. The model is developed by expanding the state vector as a mean state plus a fluctuating state. Nonlinear terms are not neglected, and the effect of nonlinearity is studied. A high order numerical method is used to compute the solution. The results show a significant spectral broadening in the flow field due to the nonlinearity. In addition, large scale structures are broken down into smaller scales. --------------- [33] Swanson, R. C. and Eli Turkel A multistage time-stepping scheme for the Navier-Stokes equations ICASE Report No. 84-62, November 23, 1984, 28 pages A class of explicit multistage time-stepping schemes is used to construct an algorithm for solving the compressible Navier-Stokes equations. Flexibility in treating arbitrary geometries is obtained with a finite-volume formulation. Numerical efficiency is achieved by employing techniques for accelerating convergence to steady state. Computer processing is enhanced through vectorization of the algorithm. The scheme is evaluated by solving laminar and turbulent flows over a flat plate and an NACA 0012 airfoil. Numerical results are compared with theoretical solutions or other numerical solutions and/or experimental data. --------------- [34] Brenier, Yann and Stanley Osher Approximate Riemann solvers and numerical flux functions ICASE Report No. 84-63, December 5, 1984, 33 pages Given a monotone function z(x) which connects two constant states, uL < uR, (uL > uR), we find the unique (up to a constant) convex (concave) flux function, F(u) such that z(x/t) is the physically correct solution to the associated Riemann problem. For z(x/t), an approximate Riemann solver to a given conservation law, we derive simple necessary and sufficient conditions for it to be consistent with any entropy inequality. Associated with any member of a general class of consistent numerical fluxes, hf(uR,uL), we have an approximate Riemann solver defined through z(Z) = (-d/dZ)hfZ(uR,uL), where fZ(u) = f(u) - Zu. We obtain the corresponding F(u) via a Legendre transform and show that it is consistent with all entropy inequalities iff hfZ(uR,uL) is an E flux for each relevant Z. Examples involving commonly used two point numerical fluxes are given, as are comparisons with related work. --------------- [35] Hall, Philip and Mujeeb R. Malik On the instability of a three-dimensional attachment line boundary layer: Weakly nonlinear theory and a numerical approach ICASE Report No. 84-64, December 6, 1984, 65 pages The instability of a three-dimensional attachment line boundary layer is considered in the nonlinear regime. Using weakly nonlinear theory, it is found that, apart from a small interval near the (linear) critical Reynolds number, finite amplitude solutions bifurcate subcritically from the upper branch of the neutral curve. The time-dependent Navier-Stokes equations for the attachment line flow have been solved using a Fourier-Chebyshev spectral method and the subcritical instability is found at wavenumbers that correspond to the upper branch. Both the theory and the numerical calculations show the existence of supercritical finite amplitude (equilibrium) states near the lower branch which explains why the observed flow exhibits a preference for the lower branch modes. The effect of blowing and suction on nonlinear stability of the attachment line boundary layer is also investigated. --------------- [36] Fix, George J. and Manil Suri Three-dimensional mass conserving elements for compressible flows ICASE Report No. 84-65, December 13, 1984, 27 pages A variety of finite element schemes has been used in the numerical approximation of compressible flows particularly in underwater acoustics. In many instances instabilities have been generated due to the lack of mass conservation. In this paper we develop new two- and three-dimensional elements which avoid these problems. --------------- [37] Banks, H. Thomas and James M. Crowley Estimation of material parameters in elastic systems ICASE Report No. 84-66, December 28, 1984, 45 pages In this paper we present theoretical and numerical results for inverse problems involving estimation of spatially varying parameters such as stiffness and damping in distributed models for elastic structures such as Euler-Bernoulli beams. An outline of algorithms we have used and a summary of some of our computational experiences are presented. --------------- [38] Ito, Kazufumi and Robert K. Powers Chansrasekhar equations for infinite dimensional systems ICASE Report No. 84-67, December 31, 1984, 36 pages In this paper we derive the Chandrasekhar equations for linear time invariant systems defined on Hilbert spaces using a functional analytic technique. An important consequence of this is that the solution to the evolutional Riccati equation is strongly differentiable in time and one can define a 'strong' solution of the Riccati differential equation. A detailed discussion on the linear quadratic optimal control problem for hereditary differential systems is also included. --------------- [39] Ortega, James M. and Robert G. Voigt Solution of partial differential equations on vector and parallel computers ICASE Report No. 85-1, January 2, 1985, 173 pages In this paper we review the present status of numerical methods for partial differential equations on vector and parallel computers. A discussion of the relevant aspects of these computers and a brief review of their development is included, with particular attention paid to those characteristics that influence algorithm selection. Both direct and iterative methods are given for elliptic equations as well as explicit and implicit methods for initial-boundary value problems. The intent is to point out attractive methods as well as areas where this class of computer architecture cannot be fully utilized because of either hardware restrictions or the lack of adequate algorithms. A brief discussion of application areas utilizing these computers is included. --------------- [40] Voigt, Robert G. Where are the parallel algorithms? ICASE Report No. 85-2, January 16, 1985, 14 pages Four paradigms that can be useful in developing parallel algorithms are discussed. These include computational complexity analysis, changing the order of computation, asynchronous computation, and divide and conquer. Each is illustrated with an example from scientific computation, and it is shown that computational complexity must be used with great care or an inefficient algorithm may be selected. --------------- [41] Gottlieb, D. and E. Tadmor Recovering pointwise values of discontinuous data within spectral accuracy ICASE Report No. 85-3, January 31, 1985, 20 pages We show how pointwise values of a function, f(x), can be accurately recovered either from its spectral or pseudospectral approximations, so that the accuracy solely depends on the local smoothness of f in the neighborhood of the point x. Most notably, given the equidistant function grid values, its intermediate point values are recovered within spectral accuracy, despite the possible presence of discontinuities scattered in the domain. (Recall that the usual spectral convergence rate decelerates otherwise to first order, throughout). To this end we employ a highly oscillatory smoothing kernel, in contrast to the more standard positive unit-mass mollifiers. In particular, post-processing of a stable Fourier method applied to hyperbolic equations with discontinuous data, recovers the exact solution modulo a spectrally small error. Numerical examples are presented. --------------- [42] Lakin, W. D. Differentiating Matrices for Arbitrarily Spaced Grid Points ICASE Report No. 85-5, January 31, 1985, 23 pages Differentiating matrices allow the numerical differentiation of functions defined at points of a discrete grid. The present work considers a type of differentiating matrix based on local approximation on a sequence of sliding subgrids. Previous derivations of this type of matrix have been restricted to grids with uniformly spaced points, and the resulting derivative approximations have lacked precision, especially at endpoints. The new formulation allows grids which have arbitrarily spaced points. It is shown that high accuracy can be achieved through use of differentiating matrices on non-uniform grids which include "near-boundary" points. Use of the differentiating matrix as an operator to solve eigenvalue problems involving ordinary differential equations is also considered. --------------- [43] Fleming, Peter J. SUBOPT - A CAD program for suboptimal linear regulators ICASE Report No. 85-6, February 4, 1985, 16 pages. This interactive software package provides design solutions for both standard linear quadratic regulator (LQR) and suboptimal linear regulator problems. Intended for time-invariant continuous systems the package is easily modified to include sampled-data systems. LQR designs are obtained by established techniques while the large class of suboptimal problems containing controller and/or performance index options is solved using a robust gradient minimization technique. Numerical examples demonstrate features of the package and recent developments are described. --------------- [44] Banks, H. Thomas and I. Gary Rosen A Galerkin method for the estimation of parameters in hybrid systems governing the vibration of flexible beams with tip bodies ICASE Report No. 85-7, February 5, 1985, 44 pages In this report we develop an approximation scheme for the identification of hybrid systems describing the transverse vibrations of flexible beams with attached tip bodies. In particular, problems involving the estimation of functional parameters (spatially varying stiffness and/or linear mass density, temporally and/or spatially varying loads, etc.) are considered. The identification problem is formulated as a least squares fit to data subject to the coupled system of partial and ordinary differential equations describing the transverse displacement of the beam and the motion of the tip bodies respectively. A cubic spline-based Galerkin method applied to the state equations in weak form and the discretization of the admissible parameter space yield a sequence of approximting finite dimensional identification problems. We demonstrate that each of the approximating problems admits a solution and that from the resulting sequence of optimal solutions a convergent subsequence can be extracted, the limit of which is a solution to the original identification problem. The approximating identification problems can be solved using standard techniques and readily available software. Numerical results for a variety of examples are provided. --------------- [45] Tal-Ezer, Hillel The eigenvalues of the pseudospectral Fourier approximation to the operator sin(2x) (d/dx) ICASE Report No. 85-8, February 7, 1985, In this note we show that the eigenvalues of the pseudospectral Fourier approximation to the operator sin(2x) (d/dx) have real parts either equal to zero or equal to plus or minus one. Whereas this does not prove stability for the Fourier method, applied to the hyperbolic equation dU/dt = sin(2x)(dU/dx) - Pi < x < - Pi; it indicates that the growth in time of the numerical solution is essentially the same as that of the solution to the differential equation. --------------- [46] Tal-Ezer, Hillel Spectral methods in time for parabolic problems ICASE Report No. 85-9, February 8, 1985, 19 pages A pseudospectral explicit scheme for solving linear, periodic, parabolic problems is described. It has infinite accuracy both in time and in space. The high accuracy is achieved while the time resolution parameter M (M = 0(1/(delta t)) for time marching algorithm) and the space resolution parameter N (N = 0(1/(delta t)) have to satisfy M = 0(N**(1 + ep)) ep > 0, compared to the common stability condition M = 0(N**2) which has to be satisfied in any explicit finite order time algorthm. --------------- [47] Tadmor, Eitan Complex symmetric matrices with strongly stable iterates ICASE Report No. 85-10, February 14, 1985, 22 pages We study complex-valued symmetric matrices. A simple expression for the spectral norm of such matrices is obtained, by utilizing a unitarily congruent invariant form. Consequently, we provide a sharp criterion for identifying those symmetric matrices whose spectral norm is not exceeding one: such strongly stable matrices are usually sought in connection with convergent difference approximations to partial differential equations. As an example, we apply the derived criterion to conclude the strong stability of a Lax- Wendroff scheme. --------------- [48] Turkel, Eli Algorithms for the Euler and Navier-Stokes equations for supercomputers ICASE Report No. 85-11, February 14, 1985, 27 pages We consider the steady state Euler and Navier-Stokes equations for both compressible and incompressible flow. Methods are found for accelerating the convergence to a steady state. This acceleration is based on preconditioning the system so that it is no longer time consistent. In order that the acceleration technique be scheme independent this preconditioning is done at the differential equation level. Applications are presented for very slow flows and also for the incompressible equations. --------------- [49] Pratt, Terrence W. PISCES: an environment for parallel scientific computation ICASE Report No. 85-12, February 15, 1985, 30 pages PISCES (Parallel Implementation of Scientific Computing EnvironmentS) is a project to provide high-level programming environments for parallel MIMD computers. Pisces 1, the first of these environments, is a Fortran 77 based environment which runs under the UNIX operating system. The Pisces 1 user programs in Pisces Fortran, an extension of Fortran 77 for parallel processing. The major emphasis in the Pisces 1 design is in providing a carefully specified "virtual machine" that defines the run-time environment within which Pisces Fortran programs are executed. Each implementation then provides the same virtual machine, regardless of differences in the underlying architecture. The design is intended to be portable to a variety of architectures. Currently Pisces 1 is implemented on a network of Apollo workstations and on a DEC VAX uniprocessor via simulation of the task level parallelism. An implementation for the Flexible Computing Corp. FLEX/32 is under construction. This paper provides an introduction to the Pisces 1 virtual computer and the Fortran 77 extensions. An example of an algorithm for the iterative solution of a system of equations is given. The most notable features of the design are the provision for several different "granularities" of parallelism in programs and the provision of a "window" mechanism for distributed access to large arrays of data. --------------- [50] Harabetian, Eduard A convergent series expansion for hyperbolic systems of conservation laws ICASE Report No. 85-13, February 20, 1985, 88 pages We consider the discontinuous piecewise analytic initial value problem for a wide class of conservation laws that includes the full three-dimensional Euler equations. The initial interaction at an arbitrary curved surface is resolved in time by a convergent series. Among other features the solution exhibits shock, contact, and expansion waves as well as sound waves propagating on characteristic surfaces. The expansion waves correspond to the one-dimensional rarefactions but have a more complicated structure. The sound waves are generated in place of zero strength shocks, and they are caused by mismatches in derivatives. --------------- [51] Abarbanel, Saul and David Gottlieb Information content in spectral calculations ICASE Report No. 85-14, February 21, 1985, 13 pages Spectral solutions of hyperbolic partial differential equations induce a Gibbs phenomenon near local discontinuities or strong gradients. A procedure is presented for extracting the piecewise smooth behavior of the solution out of the oscillatory numerical solution data. The procedure is developed from the theory of linear partial differential equations. Its application to a non-linear system (the two-dimensional Euler equations of gas dynamics) is shown to be efficacious for the particular situation. --------------- [52] Cox, Christopher L. and George J. Fix On the accuracy of least squares methods in the presence of corner singularities ICASE Report No. 85-15, February 21, 1985, 24 pages This paper treats elliptic problems with corner singularities. Finite element approximations based on variational principles of the least squares type tend to display poor convergence properties in such contexts. Moreover, mesh refinement or the use of special singular elements do not appreciably improve matters. Here we show that if the least squares formulation is done in appropriately weighted space, then optimal convergence results in unweighted spaces like L two. --------------- [53] Maday, Yvon Analysis of spectral operators in one-dimensional domains ICASE Report No. 85-17, February 28, 1985, 26 pages We prove results concerning certain projection operators on the space of all polynomials of degree less than or equal to N with respect to a class of one-dimensional weighted Sobolev spaces. These results are useful in the theory of the approximation of partial differential equations with spectral methods. --------------- [54] Roe, Philip L. Discrete models for the numerical analysis of time-dependent multidimensional gas dynamics ICASE Report No. 85-18, March 1, 1985, 36 pages This paper explores a possible technique for extending to multidimensional flows some of the upwind-differencing methods that have proved highly successful in the one-dimensional case. Attention here is concentrated on the two-dimensional case, and the flow domain is supposed to be divided into polygonal computational elements. Inside each element the flow is represented by a local superposition of elementary solutions consisting of plane waves not necessarily aligned with the element boundaries. --------------- ##### LIVERMORE ##### [55] R. C. Y. Chin, G. W. Hedstrom, F. A. Howes, and J. R. McGraw Parallel computation of multiple-scale problems Report UCRL-92007 Lawrence Livermore National Laboratory <No Abstract> Submitted by: Gerald Hedstrom <hedstrom@lll-crg.arpa> Obtainable from: Gerald Hedstrom Lawrence Livermore National Laboratory Livermore, CA 94550 --------------- ##### MARYLAND ##### All of the Maryland titles were submitted by: Rodrigo Fontecilla <rod@maryland.arpa> They may be obtained from: Rodrigo Fontecilla University of Maryland Dept. of Computer Science College Park, MD 20742 #################### [56] Rodrigo Fontecilla A globally convergent algorithm for nonlinear equality constrained optimization using a double line search. In this paper, we present a globally convergent algo- rithm based on the 2-step algorithms developed by the author. Two different line searches with Goldstein-Armijo conditions are implemented. We show that the horizontal step is a descent direction for the objective function and that the vertical step is a descent direction for the $l sub 2#- norm of the constraints. Further, we show that the total step (horizontal+vertical) is a descent direction for both, the Lagrangian and the $l sub 2#-norm of the constraints. It is shown that the algorithm generates a sequence {$x sub k#} converging in the weak sense, i.e. the gradient of the Lagrangian and the constraints converge to zero. Fast con- vergence is guaranteed close to the solution with a 2-step q-quadratic rate. The BFGS secant update is implemented, a global convergence result is also obtained, and close to the solution the convergence is 2-step q-superlinear. Numerical results using a finite difference approximation of the Hes- sian and the BFGS secant update on 30 test problems are presented. --------------- [57] Rodrigo Fontecilla SUBSPACE INVARIANCE OF BROYDEN'S $THETA$-CLASS University of Maryland Computer Science Department Technical Report 1430 Broyden's $theta#-class of updating formulae have the pro- perty that if the current approximation matrix is symmetric and positive definite, then the updated matrix maintains those same properties under certain conditions. It is shown that if the current approximation matrix is symmetric and positive definite on a subspace of $Rn#, then the updated matrix is symmetric and positive definite along the same subspace under certain conditions. Applications of these results to the implementation of quasi-Newton methods for solving nonlinear constrained optimization problems are presented. --------------- ##### MINNESOTA ##### All of the reports from Minnesota were submitted without abstracts by: Panos Pardalos They are available from him at: Computer Science Department University of Minnesota Lind Hall 136 Minneapolis MN 55455 ##################### [58] J.B.Rosen Computational Experience with Large-Scale Constrained Global Optimization Tech. Report 84-13 June 1984 --------------- [59] J.B.Rosen and P.M.Pardalos Global Minimization of Large-Scale Constrained Concave Quadratic Problems by Separable Programming Tech. Report 84-29 October 1984 --------------- [60] P.M.Pardalos and J.B.Rosen Reduction of Nonlinear Integer Separable Programming Problems to Zero-One Linear Programming Problems Tech. Report 84-25 October 1984 --------------- [61] P.M.Pardalos and J.B.Rosen Methods for Global Concave Minimization: A Bibliographic Survey Tech. Report 84-31 November 1984 --------------- [62] P.M.Pardalos and J.B.Rosen A Global Optimization Approach to the Linear Complementarity Problem December 1984 --------------- ##### OAK RIDGE ##### [63] M. T. Heath A. J. Laub C. C. Paige R. C. Ward ORNL UC-Santa Barbara McGill Univ. ORNL COMPUTING THE SINGULAR VALUE DECOMPOSITION OF A PRODUCT OF TWO MATRICES Tech. Rept. ORNL-6118, January 1985 An algorithm is developed for computing the singular value decomposition of a product of two general matrices without explicitly forming the product. The algorithm is based on an earlier Jacobi-like method due to Kogbetliantz and uses plane rotations applied to the two matrices separately. A triangular variant of the basic algorithm is developed that reduces the amount of work required. Submitted by: Bob Ward <ward@ornl-msr.arpa> Obtainable from: Mathematical Sciences Oak Ridge National Laboratory P. O. Box Y, Bldg. 9207-A Oak Ridge, Tennessee 37831 ------------------- [64] Alan George Michael T. Heath Joseph Liu University of Waterloo ORNL York University PARALLEL CHOLESKY FACTORIZATION ON A MULTIPROCESSOR Tech. Rept. ORNL-6124, March 1985 A parallel algorithm is developed for Cholesky factorization on a multiprocessor. The algorithm is based on self-scheduling of a pool of tasks. The subtasks in several variants of the basic elimination algorithm are analyzed for potential concurrency in terms of precedence relations, work profiles, and processor utilization. This analysis is supported by simulation results. The most promising variant, which we call column-Cholesky, is identified and implemented for the Denelcor HEP multiprocessor. Experimental results are given for this machine. Submitted by: Bob Ward Obtainable from: Mathematical Sciences Oak Ridge National Laboratory P. O. Box Y, Bldg. 9207-A Oak Ridge, Tennessee 37831 ----------------- [65] R. E. Funderlic C. D. Meyer, Jr. Oak Ridge National Laboratory North Carolina State University SENSITIVITY OF THE STATIONARY DISTRIBUTION VECTOR FOR AN ERGODIC MARKOV CHAIN Technical Report ORNL-6098, January 1985 Stationary distribution vectors p for Markov chains with associated transition matrices T are important in the analysis of many models in the mathematical sciences, such as queuing networks, input-output economic models and compartmental tracer analysis models. The purpose of this paper is to provide insight into the sensitivity of p to perturbations in the transition probabilities of T and to understand some of the difficulties in computing an accurate p. The group inverse A# of I-T is shown to be of fundamental importance in understanding sensitivity or conditioning of p. The main result shows that if there is a state that is accessible from every other state and the corresponding column of T has no small off-diagonal elements, then p cannot be sensitive to small perturbations in T. Ecological examples are given. A new stable algorithm for calculating A# is described. Submitted by: Bob Ward Obtainable from: Mathematical Sciences Oak Ridge National Laboratory P. O. Box Y, Bldg. 9207-A Oak Ridge, Tennessee 37831 --------------- [66] George Ostrouchov Symbolic Givens Reduction in Large Sparse Least Squares ORNL-6102 Oak Ridge National Laboratory December 1984 Orthogonal Givens factorization is a popular method for solving large sparse least squares problems. In order to exploit sparsity and to use a fixed data structure in Givens reduction, a preliminary symbolic factorization needs to be performed. Some results on row-ordering and structure of rows in a partially reduced matrix are obtained using a graph-theoretic representation. These results provide a basis for a symbolic Givens factorization. Column-ordering is also discussed, and an algorithm for symbolic Givens reduction is developed and tested. --------------- ##### UBC ##### [67] Manfred R. Trummer Nystroem's method versus Fourier type methods for the numerical solution of integral equations. Tech. Rep. 84-23 It is shown that Nystroem's method and Fourier type methods produce the same approximation to a solution of an integral equation at the collocation points for Nystroem's method. The quadrature rule for numerical integration must have these collocation points as abscissa. Submitted by: Manfred R. Trummer <trummer@ubc.csnet> Obtainable from: Department of Computer Science The University of British Columbia Vancouver, B.C., Canada V6T 1W5 --------------- [68] Manfred R. Trummer An efficient implementation of a conformal mapping method using the Szego kernel. Tech. Rep. 84-24 An implementation, based on iterative techniques, of a method to compute the Riemann mapping function is presented. The method has been recently introduced by N. Kerzman and the author. It expresses the Szego kernel as the solution of an integral equation of the second kind. It is shown how to treat symmetric regions. The algorithm is tested on five examples. The numerical results show that the method is competitive, with respect to accuracy, stability, and efficiency. Submitted by: Manfred R. Trummer <trummer@ubc.csnet> Obtainable from: Department of Computer Science The University of British Columbia Vancouver, B.C., Canada V6T 1W5 --------------- [69] Uri Ascher and G. Bader STABILITY OF COLLOCATION AT GAUSSIAN POINTS Technical Report 84-13 Revised, February, 1985 Symmetric Runge-Kutta schemes are particularly useful for solving stiff two-point boundary value problems. Such A-stable schemes perform well in many cases, but it is demonstrated that in some instances the stronger property of algebraic stability is required. A characterization of symmetric, algebraically stable Runge-Kutta schemes is given. The class of schemes thus defined turns out not to be very rich: The only collocation schemes in it are those based on Gauss points, and other schemes in the class do not seem to offer any advantage over collocation at Gaussian points. Submitted by: Uri Ascher <ascher@ubc.csnet> Obtainable from: Uri Ascher Department of Computer Science University of British Columbia Vancouver, British Columbia Canada V6T 1W5 --------------- [70] Uri Ascher COLLOCATION FOR TWO-POINT BOUNDARY VALUE PROBLEMS REVISITED Technical Report 84-17 November, 1984 Collocation methods for two-point boundary value problems for higher order differential equations are considered. By using appropriate monomial bases, we relate these methods to corresponding one-step schemes for 1st order systems of differential equations. This allows us to present the theory for nonstiff problems in relatively simple terms, refining at the same time some convergence results and discussing stability. No restriction is placed on the meshes used. Submitted by: Uri Ascher <ascher@ubc.csnet> Obtainable from: Uri Ascher Department of Computer Science University of British Columbia Vancouver, British Columbia Canada V6T 1W5 --------------- ##### WATERLOO ##### [71] P. H. Calamai and A. R. Conn A projected Newton method for $l sub p$ location problems. The University of Waterloo Computer Science Department Tech Rept CS-85-07 This paper is concerned with the numerical solution of con- tinuous minisum multifacility location problems involving the $l sub p$ norm, where 1<p<infinity This class of problems is poten- tially difficult to solve because the objective function is not everywhere differentiable. After developing conditions that characterize the minimum of the problem under consideration, a second-order algorithm is presented. This algorithm is based on the solution of a finite sequence of linearly constrained sub- problems. Descent directions for these subproblems are obtained by projecting the Newton direction onto the corresponding con- straint manifold. Univariate minimization is acheived via a spe- cialized line search which recognizes the possibility of first derivative discontinuity ( and second derivative unboundedness) at point along the descent direction. Th algorithm, motivated by earlier work by Calamai and Calamai and Conn and related to methods recently described by Overton and Dax, is shown to pos- sess both global and quadratic convergence properties. Degeneracy can complicate the numerical solution of the sub- problems. This degeneracy is identified, and a method that cir- cumvents this degeneracy is included. An implementation of the algorithm, that exploits the intrinsic structure of the location problem formulation, is then described along with a discusion of numerical results. Submitted by: Andrew R. Conn <arconn%water@waterloo.csnet> Obtainable from: Andrew R. Conn The University of Waterloo Computer Science Department Waterloo, Ontario N2L 3G1 Canada --------------- [72] R. H. Bartels and J. C. Beatty Beta-Splines with a Difference University of Waterloo Computer Science Department Technical Report CS-83-40 Local control of the shape parameters beta-1 and beta-2 in a beta-spline has previously relied on quintic Hermite interpolation of the distinct beta values associated with the joints in a geometrically-continuous piecewise parametric polynomial curve. Here we introduce an alternative means of obtaining local control of geometrically continuous piecewise polynomial curves. The essential idea is to generalize the truncated power functions suitably to account for the discontinuities allowed by geometric continuity, and to generate the beta-splines from a differencing process. The resulting basis functions are piecewise cubic, partition unity, have local support, and are nonegative for "reasonable" values of beta-1 and beta-2. They accommodate variable knot spacing as well as different beta values at each knot. Submitted by: Richard H. Bartels <rhbartels@waterloo.csnet> Obtainable from: Richard H. Bartels or John C. Beatty The University of Waterloo Computer Science Department Waterloo, Ontario N2L 3G1 Canada --------------- ##### YALE ##### [73] Howard C. Elman A Stability Analysis of Incomplete LU Factorizations Report YALEU/DCS/RR-365 <No Abstract> Submitted by: Howard C. Elman <elman%yale-bulldog@yale.arpa> Obtainable from: Howard C. Elman Yale University New Have, Conn. --------------- From GOLUB@SU-SCORE.ARPA Fri Aug 30 14:45:54 1985 Received: from SU-SCORE.ARPA (su-score.arpa.ARPA) by anl-mcs.ARPA ; Fri, 30 Aug 85 14:45:11 cdt Return-Path: <rhbartels%watcgl%waterloo.csnet@csnet-relay.arpa> Received: from csnet-relay by SU-SCORE.ARPA with TCP; Tue 27 Aug 85 17:48:23-PDT Received: from waterloo by csnet-relay.csnet id ae16675; 27 Aug 85 20:33 EDT Date: Tue, 27 Aug 85 13:12:38 edt From: Richard Bartels <rhbartels%watcgl%waterloo.csnet@csnet-relay.arpa> Message-Id: <8508271712.AA25227@watcgl.UUCP> To: na@SU-SCORE.ARPA Subject: T&A Part 1 of 3 Resent-Date: Fri 30 Aug 85 11:47:55-PDT Resent-From: Gene Golub (415/497-3124) <GOLUB@SU-SCORE.ARPA> Resent-To: :; Resent-Message-Id: <12139307538.36.GOLUB@SU-SCORE.ARPA> Status: R ================================================== NUMERICAL ANALYSIS TITLES -- Volume 3, Number 2, Part 1 of 3 -- August 27, 1985 ================================================== Please note: Due to its length, this list is being distributed in 3 parts, each part is about 7 pages in length. ##### AT&T BELL LABORATORIES ##### For copies of these reports, send mail to {ucbvax,ihnp4}!research!wmc or wmc.btl@csnet-relay or na.coughran@su-score or Bill Coughran AT&T Bell Labs 2C419 Murray Hill, NJ 07974 (201) 582-6619 ################### [1] R. E. Bank, W. M. Coughran, Jr., W. Fichtner, D. J. Rose, and R. K. Smith, COMPUTATIONAL ASPECTS OF SEMICONDUCTOR DEVICE SIMULATION NAMs 85-3 In this chapter, we present an overview of the numerical techniques used to solve the coupled system of nonlinear partial differential equations that model the behavior of semiconductor devices. These methods have been incorporated into our device simulation package which has been success- fully used to model complex device structures in two and three space dimensions for both steady-state and transient conditions. --------------- [2] R. E. Bank, W. M. Coughran, Jr., W. Fichtner, E. H. Grosse, D. J. Rose, R. K. Smith, TRANSIENT SIMULATION OF SILICON DEVICES AND CIRCUITS NAMs 85-8 In this paper, we present an overview of the physical prin- ciples and numerical methods used to solve the coupled sys- tem of nonlinear partial differential equations that model the transient behavior of silicon VLSI device structures. We also describe how the same techniques are applicable to circuit simulation. A composite linear multistep formula is introduced as the time-integration scheme. Newton-iterative methods are exploited to solve the nonlinear equations that arise at each time step. We also present a simple data structure for nonsymmetric matrices with symmetric nonzero structures that facilitates iterative or direct methods with substantial efficiency gains over other storage schemes. Several computational examples, including a CMOS latchup problem, are presented and discussed. --------------- ##### CORNELL ##### [3] Franklin T. Luk and Sanzheng Qiao A fast but unstable orthogonal triangularization technique for Toeplitz matrices EE-CEG-85-5 Address: School of Electrical Engineering Cornell University Ithaca, NY 14863 < No Abstract > Submitted by: Frank Luk (luk%tesla@cornell.csnet) Obtainable from: Frank Luk Cornell University Ithaca, New York 14853 --------------- [4] Franklin T. Luk Algorithm-based Fault Tolerance for Parallel Matrix Equation Solvers EE-CEG-85-2 School of Electrical Engineering, Cornell University Ithaca, New York 14853 June 1985 We examine the checksum schemes of Abraham et al. for the computation of the LU-factorization using a multiprocessor array. Their methods are very efficient for detecting a transient error, but quite expensive for correcting it due to the need for a computation rollback. In this paper, we show how to avoid the rollback and how to implement pivoting. We also introduce a new checksum method for solving triangular sets of linear equations. Obtainable from: Frank Luk, above. --------------- ##### COURANT INSTITUTE ##### [5] O. McBryan Computational Methods for Discontinuities in Fluids Lectures in Applied Mathematics vol. 22, AMS, Providence, 1985. < No Abstract > Submitted by: Oliver McBryan (mcbryan@nyu.arpa) Obtainable from: Oliver McBryan Courant Institute 251 Mercer St New York, NY 10012 --------------- [6] O. McBryan and E. Van de Velde Parallel Algorithms for Elliptic Equations Proceedings of the 1984 ARO Novel Computing Environments Conference Stanford University, SIAM , to appear. < No Abstract > Obtainable from: O. McBryan, above. --------------- [7] O. McBryan, E. Van de Velde, and P. Vianna Parallel Algorithms for Elliptic and Parabolic Equations Proceedings of the Conference on Parallel Computations in Heat Transfer and Fluid Flows University of Maryland, November 1984. < No Abstract > Obtainable from: O. McBryan, above. --------------- [8] O. McBryan Using CRAY super-computers as Attached Processors Courant Institute Preprint, 1985. < No Abstract > Obtainable from: O. McBryan, above. --------------- [9] O. McBryan State of the Art of Multiprocessors in Scientific Computation Proceedings of European Weather Center Conference on Multiprocessors Reading, England, Dec 1984, to appear. < No Abstract > Obtainable from: O. McBryan, above. --------------- [10] O. McBryan and E. Van de Velde Parallel Algorithms for Elliptic Equation Solution on the HEP Computer Proceedings of the First HEP Conference University of Oklahoma, March 1985. < No Abstract > Obtainable from: O. McBryan, above. --------------- [11] O. McBryan and E. Van de Velde Parallel Algorithms for Elliptic Equations to appear in Commun. Pure and Appl. Math., Feb 1985. < No Abstract > Obtainable from: O. McBryan, above. --------------- [12] O. McBryan and E. Van de Velde Elliptic Equation Algorithms on Parallel Computers Proceedings of the Conference on Parallel Computers and Partial Differential Equations, Commun. in Applied Numerical Methods, University of Texas, Austin, March 1985, to appear. < No Abstract > Obtainable from: O. McBryan, above. --------------- [13] J. Glimm, B. Lindquist, O. McBryan, and G. Tryggvason Sharp and Diffuse Fronts in Oil Reservoirs: Front Tracking and Capillarity Proceedings of the Houston SPE/SIAM meeting on Mathematics of Reservoir Simulation, SIAM, to appear, Feb. 1985. < No Abstract > Obtainable from: O. McBryan, above. --------------- [14] J. Glimm and O. McBryan A Computational Model for Interfaces Courant Institute Preprint, 1985. < No Abstract > Obtainable from: O. McBryan, above. --------------- [15] James W. Demmel and Bo Kagstrom Computing Stable Eigendecompositions of Matrix Pencils Technical Report # 163 Computer Science Department Courant Institute May, 1985 We discuss how to compute an eigendecomposition of a matrix pencil A-zB when A and B are only known to within a tolerance epsilon. When A-zB is regular (i.e. det(A-zB) is not identically zero) we show how to partition the spectrum and eigenspaces of A-zB into clusters which vary smoothly as A-zB varies within a ball of radius epsilon. When A-zB is singular (a case of interest in systems theory) so that the structures we wish to compute are nongeneric, we show that certain spaces and eigenvalues of the pencil vary smoothly if A-zB varies along a lower dimensional surface as well as within a ball of radius epsilon. This result implies that the usual algorithms for analyzing singular pencils generally compute accurate eigenvalues and spaces. We apply this result by computing perturbation bounds for the controllable subspace and uncontrollable modes of a system dx/dt = Cx + Du. Submitted by: James Demmel (demmel.csd2@nyu) Obtainable from: James Demmel Courant Institute 251 Mercer Str. New York, NY 10012 --------------- [16] Jonathan Goodman, Robert Kohn, and Luis Reyna A Numerical Study of a Relaxed Variational Problem We present the numerical solution of an optimization problem that arises in two phase flow, and in the design of a beam out of two different materials in given proportion for optimum torsional rigidity. The problem is to minimize a Dirichlet integral over all possible functions and over all possible subdomains of a given area. Minimizing over functions with the subdomain fixed would lead to a second order elliptic equation with discontinuous coefficients. A naive (but very plausible) algorithm based directly on this formulation is shown to have "premature termination" at points that are not even stationary points. It is known that problems like this often do not have "classical" solutions. Rather, there are "weak" or "relaxed" solutions that involve microscopic mixing of the two materials. Application of homogenization theory gives a relaxed minimization problem that can numerically be solved. We used a variational multigrid proceedure to get high resolution (256 x 256) in reasonable time on a VAX-780. The numerical results show that the region of microscopic mixing can occupy about 15% of the region in some cases. Submitted by: Jonathan Goodman (goodnan@nyu-acf1.csnet) Obtainable from: Jonathan Goodman Courant Institute of Mathematical Sciences 251 Mercer St. New York, New York, 10012 --------------- [17] Jonathan Goodman, Robert Kohn, and Luis Reyna A Numerical Study of a Relaxed Variational Problem We present the numerical solution of an optimization problem that arises in two phase flow, and in the design of a beam out of two different materials in given proportion for optimum torsional rigidity. The problem is to minimize a Dirichlet integral over all possible functions and over all possible subdomains of a given area. Minimizing over functions with the subdomain fixed would lead to a second order elliptic equation with discontinuous coefficients. A naive (but very plausible) algorithm based directly on this formulation is shown to have "premature termination" at points that are not even stationary points. It is known that problems like this often do not have "classical" solutions. Rather, there are "weak" or "relaxed" solutions that involve microscopic mixing of the two materials. Application of homogenization theory gives a relaxed minimization problem that can numerically be solved. We used a variational multigrid proceedure to get high resolution (256 x 256) in reasonable time on a VAX-780. The numerical results show that the region of microscopic mixing can occupy about 15% of the region in some cases. Obtainable from J. Goodman, above --------------- [18] Jonathan Goodman Convergence of the Random Vortex Method The random vortex method, introduced by Chorin, was intended to compute high Reynolds number incompressible flows in 2 or 3 dimensions. We prove that the method converges (with probability one) for smooth flows without boundaries in two dimensions. The error bounds are independent of the viscosity. The proof relies on the form of smoothing (replacing vortex points by vortex blobs) introduced by Hald in his convergence proof for the (non-random) vortex method for the incompressible Euler equations in 2 dimensions. Obtainable from J. Goodman, above --------------- ##### EMORY ##### [19] Henry Wolkowicz and Phil Smith A NONLINEAR EQUATION FOR LINEAR PROGRAMMING Research Report 16 Emory University We present a characterization of the 'normal' optimal solution of the linear programming problem. This characterization involves the solution of m piecewise linear equations in m unknowns. Submitted by: Henry Wolkowicz csnet: henry@emory bitnet: henry_wolkowicz@uqv-mts.bitnet uucp: alberta!uqv-mts!sunn emory!henry Obtainable from: Henry Wolkowicz Department of Mathematics and Computer Science Emory University Atlanta, Georgia 30322 --------------- From GOLUB@SU-SCORE.ARPA Fri Aug 30 15:40:44 1985 Received: from SU-SCORE.ARPA (su-score.arpa.ARPA) by anl-mcs.ARPA ; Fri, 30 Aug 85 15:39:53 cdt Return-Path: <rhbartels%watcgl%waterloo.csnet@csnet-relay.arpa> Received: from csnet-relay by SU-SCORE.ARPA with TCP; Tue 27 Aug 85 17:57:31-PDT Received: from waterloo by csnet-relay.csnet id af16675; 27 Aug 85 20:36 EDT Date: Tue, 27 Aug 85 13:13:16 edt From: Richard Bartels <rhbartels%watcgl%waterloo.csnet@csnet-relay.arpa> Message-Id: <8508271713.AA25261@watcgl.UUCP> To: na@SU-SCORE.ARPA Subject: T&A Part 2 of 3 Resent-Date: Fri 30 Aug 85 11:48:00-PDT Resent-From: Gene Golub (415/497-3124) <GOLUB@SU-SCORE.ARPA> Resent-To: :; Resent-Message-Id: <12139307553.36.GOLUB@SU-SCORE.ARPA> Status: RO ================================================== NUMERICAL ANALYSIS TITLES -- Volume 3, Number 2, Part 2 of 3 -- August 27, 1985 ================================================== Please note: Due to its length, this list is being distributed in 3 parts, each part is about 7 pages in length. ##### GMD, WEST GERMANY ##### All the papers of the GMD may be ordered from Dr. Kurt Brand Gesellschaft fuer Mathematik und Datenverarbeitung (GMD) F1/T Postfach 1240 D-5205 St. Augustin 1 F.R.G. or send a request by electronic mail to gmap18%dbngmd21.bitnet@wiscvm.arpa No abstracts were submitted. ################### [20] K. Becker Ein Mehrgitterverfahren zur Berechnung subsonischer Potential- stroemungen um Tragflaechenprofile Berichte der Gesellschaft fuer Mathematik und Datenverarbeitung, Bericht Nr. 152, Oldenbourg-Verlag, 1985 --------------- [21] K. Becker, U. Trottenberg Development of Multigrid Algorithms for Problems from Fluid Dynamics Arbeitspapiere der GMD no.111, Gesellschaft fuer Mathematik und Datenverarbeitung, Bonn, 1984. (DM 6,-) --------------- [22] A. Brandt Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics GMD-Studien no. 85, Gesellschaft fuer Mathematikund Daten- verarbeitung, Bonn, 1984. (DM 39,-) --------------- [23] O. Kolp Parallelisierung eines Mehrgitterverfahrens fuer einen Baumrechner Arbeitspapiere der GMD no. 82, Bonn, 1984. (DM 10,-) --------------- [24] R. Lorenz Some New Periodic Chebycheff Spaces Arbeitspapiere der GMD no. 79, Gesellschaft fuer Mathematik und Datenverarbeitung, Bonn, 1984 --------------- [25] J. Ruge, K. Stueben Efficient Solution of Finite Difference and Finite Element Equations by Algebraic Multigrid (AMG) Arbeitspapiere der GMD no. 89, Gesellschaft fuer Mathematik und Datenverarbeitung, Bonn, 1984. (DM 10,-) --------------- [26] B. Ruttmann, K. Solchenbach A Multigrid Solver for the computation of in-cylinder turbulent flows in engines Efficient Solutions of Elliptic Systems, Proceedings of GAMM-Seminar, Kiel, January 27 to 29, 1984, (W. Hackbusch ed.). Notes on Numerical Fluid Mechanics, 10, pp. 87-108, Vieweg, Braunschweig, 1984. --------------- [27] K. Stueben, U. Trottenberg, K. Witsch Software Development Based on Multigrid Techniques Arbeitspapiere der GMD no. 84, Gesellschaft fuer Mathematik und Datenverarbeitung, Bonn, 1984. (DM 10,-) --------------- [28] C.-A. Thole, U. Trottenberg Basic Smoothing Procedures for the Multigrid Treatment of Elliptic 3D-Operators Arbeitspapiere der GMD no. 141, Bonn, 1985. (DM 10,-) --------------- [29] U. Trottenberg Schnelle Loesung elliptischer Differentialgleichungen nach dem Mehrgitterprinzip GAMM-Mitteilungen, 1/84, February 1984. --------------- [30] U. Trottenberg, P. Wypior (Ed.) Rechnerarchitekturen fuer die numerische Simulation auf der Basis superschneller Loesungsverfahren I Workshop "Rechnerarchitektur", Erlangen, 14.-15. Juni 1984. GMD-Studien no. 88, Gesellschaft fuer Mathematik und Daten- verarbeitung, Bonn, 1984. (DM 48,-) --------------- [31] U. Trottenberg, P. Wypior (Ed.) Rechnerarchitekturen fuer die numerische Simulation auf der Basis superschneller Loesungsverfahren II Workshop "Simulation/Anwendungen und Numerik", GMD-Studien no. 102, Gesellschaft fuer Mathematik und Daten- verarbeitung, Bonn, 1985. (DM 48,-) --------------- [32] U. Trottenberg Mehrgitterprinzip und Rechnerarchitektur GAMM-Mitteilungen, 1/84, February 1984. --------------- ##### ICASE ##### All ICASE reports listed here may be obtained by writing to: Barbara Kraft ICASE, M/S 132C NASA Langley Research Center Hampton, VA 23665 ################### [33] Zang, T. A. and M. Y. Hussaini A three-dimensional spectral algorithm for simulations of transition and turbulence ICASE Report No. 85-19, March 6, 1985, 39 pages. AIAA Paper No. 85-0296 presented at the AIAA 23rd Aerospace Sciences Meeting, January 14-17, 1985, Reno, Nevada. A spectral algorithm for simulating three-dimensional, incompressible, parallel shear flows is described. It applies to the channel, to the parallel boundary layer, and to other shear flows with one wall-bounded and two periodic directions. Representative applications to the channel and to the heated boundary layer are presented. --------------- [34] Drummond, J. P., M. Y. Hussaini, and T. A. Zang Spectral methods for modeling supersonic chemically reacting flow fields ICASE Report No. 85- 20, March 8, 1985, 37 pages. AIAA J. A numerical algorithm has been developed for solving the equations describing chemically reacting supersonic flows. The algorithm employs a two- stage Runge-Kutta method for integrating the equations in time and a Chebyshev spectral method for integrating the equations in space. The accuracy and efficiency of the technique have been assessed by comparison with an existing implicit finite-difference procedure for modeling chemically reacting flows. The comparison showed that the new procedure yielded equivalent accuracy on much coarser grids as compared to the finite-difference procedure with resultant significant gains in computational efficiency. --------------- [35] LeVeque, R. J Intermediate boundary conditions for LOD, ADI, and approximate factorization methods ICASE Report No. 85-21, April 17, 1985, 24 pages. Submitted to J. Comput. Phys. A general approach to determining the correct intermediate boundary conditions for dimensional splitting methods is presented and illustrated. The intermediate solution U* is viewed as a second-order accurate approximation to a modified equation. Deriving the modified equation and using the relationship between this equation and the original equation allows us to determine the correct boundary conditions for U*. To illustrate this technique, we apply it to LOD and ADI methods for the heat equation in two and three space dimensions. The approximate factorization method is considered in slightly more generality. --------------- [36] Rosen, I. G. Spline-based Rayleigh-Ritz methods for the approximation of the natural modes of vibration for flexible beams with tip bodies ICASE Report No. 85-22, March 18, 1985, 31 pages. Submitted to Quarterly of Applied Mathematics. Rayleigh-Ritz methods for the approximation of the natural modes for a class of vibration problems involving flexible beams with tip bodies using subspaces of piecewise polynomial spline functions are developed. An abstract operator theoretic formulation of the eigenvalue problem is derived and spectral properties investigated. The existing theory for spline-based Rayleigh-Ritz methods applied to elliptic differential operators and the approximation properties of interpolatory splines are used to argue convergence and establish rates of convergence. An example and numerical results are discussed. --------------- [37] Bayliss, A., L. Maestrello, P. Parikh, and E. Turkel Numerical simulation of boundary layer excitation by surface heating/cooling ICASE Report No. 85-23, March 25, 1985, 22 pages. AIAA Paper No. 85-0565, AIAA Shear Flow Control Conference, March 12-14, 1985, Boulder, CO. This paper is a numerical study of the concept of active control of growing disturbances in an unstable compressible flow by using time periodic, localized surface heating. The simulations are calculated by a fourth-order accurate solution of the compressible, laminar Navier-Stokes equations. Fourth-order accuracy is particularly important for this problem because the solution must be computed over many wavelengths. The numerical results demonstrate the growth of an initially small fluctuation into the nonlinear regime where a local breakdown into smaller scale disturbances can be observed. It is shown that periodic surface heating over a small strip can reduce the level of the fluctuation provided that the phase of the heating current is properly chosen. --------------- [38] Hossain, M., G. Vahala, and D. Montgomery Forced MHD turbulence in a uniform external magnetic field ICASE Report No. 85-24, March 28, 1985, 39 pages. Submitted to Phys. Fluid. Two-dimensional dissipative MHD turbulence is randomly driven at small spatial scales and is studied by numerical simulation in the presence of a strong uniform external magnetic field. A novel behavior is observed which is apparently distinct from the inverse cascade which prevails in the absence of an external magnetic field. The magnetic spectrum becomes dominated by the three longest-wavelength Alfv'n waves in the system allowed by the boundary conditions: those which, in a box size of edge 2 pi, have wave numbers (kx, ky) = (1, 0), (1, 1), and (1, -1), where the external magnetic field is in the x direction. At any given instant, one of these three modes dominates the vector potential spectrum, but they do not constitute a resonantly coupled triad. Rather, they are apparently coupled by the smaller- scale turbulence. --------------- [39] Davis, S. F. Shock capturing ICASE Report No. 85-25, April 26, 1985, 23 pages. To appear in Numerical Methods for Partial Differential Equations, (S. I. Hariharan and T. H. Moulder, eds.), Pitman Press, 1986. This chapter describes recent developments which have improved our understanding of how finite difference methods resolve discontinuous solutions to hyperbolic partial differential equations. As a result of this understanding improved shock capturing methods are currently being developed and tested. Some of these methods are described and numerical results are presented showing their performance on problems containing shocks in one and two dimensions. We begin this discussion by defining what is meant by a conservative difference scheme and showing that conservation implies that, except in very special circumstances, shocks must be spread over at least two grid intervals. These two interval shocks are actually attained in one dimension if the shock is steady and an upwind scheme is used. By analyzing this case, we determine the reason for this excellent shock resolution and use this result to provide a mechanism for improving the resolution of two-dimensional steady shocks. Unfortunately, this same analysis shows that these results cannot be extended to shocks which move relative to the computing grid. To deal with moving shocks and contact discontinuities we introduce total variation diminishing (TVD) finite difference schemes and flux limiters. We show that TVD schemes are not necessarily upwind, but that upwind TVD schemes perform better because they permit a wider choice of flux limiters. The advantage of non-upwind TVD schemes is that they are easy to implement. Indeed, it is possible to add an appropriately chosen artificial viscosity to a conventional scheme such as MacCormack's method and make it TVD. We conclude by presenting some theoretical results on flux limiters and some numerical computations to illustrate the theory. --------------- [40] Brandt, A. and S. Ta'asan Multigrid solutions to quasi-elliptic schemes. ICASE Report No. 85-26, May 3, 1985, 21 pages. Progress and Supercomputing in Computational Fluid Dynamics, (Earl. S. Murman and Saul Abarbanel, eds.), Birkhauser Boston, Inc., (tentative publication date: August 20, 1985). Quasi-elliptic schemes arise from central differencing or finite element discretization of elliptic systems with odd order derivatives on non-staggered grids. They are somewhat unstable and less accurate then corresponding staggered-grid schemes. When usual multigrid solvers are applied to them, the asymptotic algebraic convergence is necessarily slow. Nevertheless, it is shown by mode analyses and numerical experiments that the usual FMG algorithm is very efficient in solving quasi-elliptic equations to the level of truncation errors. Also, a new type of multigrid algorithm is presented, mode analyzed and tested, for which even the asymptotic algebraic convergence is fast. The essence of that algorithm is applicable to other kinds of problems, including highly indefinite ones. --------------- [41] Zang, T. A. and M. Y. Hussaini On spectral multigrid methods for the time-dependent Navier-Stokes equations ICASE Report No. 85-27, May 13, 1985, 24 pages. Presented at the 2nd Copper Mountain Conference on Multigrid Methods, April 1-3, 1985, Copper Mountain, CO. A new splitting scheme is proposed for the numerical solution of the time-dependent, incompressible Navier-Stokes equations by spectral methods. A staggered grid is used for the pressure, improved intermediate boundary conditions are employed in the split step for the velocity, and spectral multigrid techniques are used for the solution of the implicit equations. --------------- [42] Osher, S. and E. Tadmor On the convergence of difference approximations to scalar conservation laws ICASE Report No. 85-28, May 14, 1985, 70 pages. Submitted to Math. Comp. We present a unified treatment of explicit in time, two level, second order resolution, total variation diminishing, approximations to scalar conservation laws. The schemes are assumed only to have conservation form and incremental form. We introduce a modified flux and a viscosity coefficient and obtain results in terms of the latter. The existence of a cell entropy inequality is discussed and such an equality for all entropies is shown to imply that the scheme is an E scheme on monotone (actually more general) data, hence at most only first order accurate in general. Convergence for TVD-SOR schemes approximating convex or concave conservation laws is shown by enforcing a single discrete entropy inequality. --------------- [43] Mehrotra, P. and J. Van Rosendale The BLAZE language: A parallel language for scientific programming ICASE Report No. 85-29, May 15, 1985, 57 pages. Submitted to Parallel Computing. Programming multiprocessor parallel architectures is a complex task. This paper describes a Pascal-like scientific programming language, Blaze, designed to simplify this task. Blaze contains array arithmetic, "forall" loops, and APL-style accumulation operators, which allow natural expression of fine grained parallelism. It also employs an applicative or functional procedure invocation mechanism, which makes it easy for compilers to extract coarse grained parallelism using machine specific program restructuring. Thus Blaze should allow one to achieve highly parallel execution on multiprocessor architectures, while still providing the user with conceptually sequential control flow. A central goal in the design of Blaze is portability across a broad range of parallel architectures. The multiple levels of parallelism present in Blaze code, in principle, allows a compiler to extract the types of parallelism appropriate for the given architecture, while neglecting the remainder. This paper describes the features of Blaze, and shows how this language would be used in typical scientific programming. --------------- [44] Trefethen, L. N. and L. Halpern Well-Posedness of one-way wave equations and absorbing boundary conditions ICASE Report No. 85-30, June 10, 1985, 23 pages. Submitted to Math. Comp. A one-way wave equation is a partial differential equation which, in some approximate sense, behaves like the wave equation in one direction but permits no propagation in the opposite one. The construction of such equations can be reduced to the approximation of the square root of 1 - s2 on [-1,1] by a rational function r(s) = Pm(s)/qn(s). This paper characterizes those rational functions r for which the corresponding one-way wave equation is well-posed, both as a partial differential equation and as an absorbing boundary condition for the wave equation. We find that if r(s) interpolates the square root of 1 - s2 at sufficiently many points in (-1,1), then well- posedness is assured. It follows that absorbing boundary conditions based on Pade approximation are well-posed if and only if (m,n) lies in one of two distinct diagonals in the Pade table, the two proposed by Engquist and Majda. Analogous results also hold for one-way wave equations derived from Chebyshev or least-squares approximation. --------------- [45] Majda, G. A new theory for multistep discretizations of stiff ordinary differential equations: Stability with large step sizes ICASE Report No. 85-31, 70 pages. Submitted to Math. Comp. In this paper we consider a large set of variable coefficient linear systems of ordinary differential equations which possess two different time scales, a slow one and a fast one. A small parameter epsilon characterizes the stiffness of these systems. We approximate a system of o.d.e.s in this set by a general class of multistep discretizations which includes both one- leg and linear multistep methods. We determine sufficient conditions under which each solution of a multistep method is uniformly bounded, with a bound which is independent of the stiffness of the system of o.d.e.s, when the step size resolves the slow time scale but not the fast one. We call this property stability with large step sizes. The theory presented in this paper lets us compare properties of one-leg methods and linear multistep methods when they approximate variable coefficient systems of stiff o.d.e.s. In particular, we show that one-leg methods have better stability properties with large step sizes than their linear multistep counterparts. This observation is consistent with results obtained by Dahlquist and Lindberg {11}, Nevanlinna and Liniger {32} and van Veldhuizen {41}. Our theory also allows us to relate the concept of D- stability (van Veldhuizen {41}) to the usual notions of stability and stability domains and to the propagation of errors for multistep methods which use large step sizes. --------------- [46] Banks, H. T. On a variational approach to some parameter estimation problems. ICASE Report No. 85-32, 38 pages. Invited lecture, Internationl Conference on Control Theory for Distributed Parameter Systems and Applications, July 9 - 14, 1984, Voran, Austria. We consider examples (1-D seismic, large flexible structures, bioturbation, nonlinear population dispersal) in which a variational setting can provide a convenient framework for convergence and stability arguments in parameter estimation problems. --------------- [47] Hariharan, S. I. Absorbing boundary conditions for exterior problems ICASE Report No. 85-33, July 18, 1985, 33 pages. To appear in Numerical Methods for Partial Differential Equations, (S. I Hariharan and T. H. Moulden, eds.), Pitman Press, 1986. In this paper we consider elliptic and hyperbolic problems in unbounded regions. These problems, when one wants to solve them numerically, have the difficulty of prescribing boundary conditions at infinity. Computationally, one needs a finite region in which to solve these problems. The corresponding conditions at infinity imposed on the finite distance boundaries should dictate the boundary condition at infinity and be accurate with respect to the interior numerical scheme. Such boundary conditions are commonly referred to as absorbing boundary conditions. This paper presents a survey and covers our own treatment on these boundary conditions for wave-like equations. --------------- From GOLUB@SU-SCORE.ARPA Fri Aug 30 16:44:09 1985 Received: from SU-SCORE.ARPA (su-score.arpa.ARPA) by anl-mcs.ARPA ; Fri, 30 Aug 85 16:43:21 cdt Return-Path: <rhbartels%watcgl%waterloo.csnet@csnet-relay.arpa> Received: from csnet-relay by SU-SCORE.ARPA with TCP; Wed 28 Aug 85 08:22:34-PDT Received: from waterloo by csnet-relay.csnet id ac19806; 28 Aug 85 7:31 EDT Date: Tue, 27 Aug 85 13:13:45 edt From: Richard Bartels <rhbartels%watcgl%waterloo.csnet@csnet-relay.arpa> Message-Id: <8508271713.AA25302@watcgl.UUCP> To: na@SU-SCORE.ARPA Subject: T&A Part 3 of 3 Resent-Date: Fri 30 Aug 85 11:48:03-PDT Resent-From: Gene Golub (415/497-3124) <GOLUB@SU-SCORE.ARPA> Resent-To: :; Resent-Message-Id: <12139307563.36.GOLUB@SU-SCORE.ARPA> Status: RO ================================================== NUMERICAL ANALYSIS TITLES -- Volume 3, Number 2, Part 3 of 3 -- August 27, 1985 ================================================== Please note: Due to its length, this list is being distributed in 3 parts, each part is about 7 pages in length. ##### JOHNS HOPKINS ##### [48] Ralph Byers and Stephen Nash On the Singular "Vectors" of the Lyapunov Operator Tech Rep 438 Mathematical Sciences Department The Johns Hopkins University June 1985 For a real matrix A, the separation of A' and A is sep(A',-A) = min ||A'X + XA||/||X||, where ||.|| represents the Frobenius norm, and A' is the transpose of A. We discuss the conjecture that the minimizer X is symmetric. This conjecture is related to the numerical stability of methods for solving the matrix Lyapunov equation. The quotient is minimized by either a symmetric matrix or a skew symmetric matrix and is maximized by a symmetric matrix. The conjecture is true if A is 2-by-2, if A is normal, if the minimum is zero, or if the real parts of the eigenvalues of A are of one sign. In general the conjecture is false, but counter examples suggest that symmetric matrices are nearly optimal. Submitted by: Steve Nash (sgn@jhu.csnet) Obtainable from: Steve Nash Department of Mathematical Sciences The Johns Hopkins University Baltimore, MD 21218 --------------- ##### MARYLAND ##### [49] Dianne P. O'Leary, Systolic Arrays for Matrix Transpose and Other Reorderings Computer Science Department University of Maryland TR-1481, March, 1985. In this note, a systolic array is described for computing the transpose of an n-by-n matrix in time 3n-1 using n squared switching processors and n squared bit buffers. A one-dimensional implementation is also described. Arrays are also given to take a matrix in by rows and put it out by diagonals, and vice versa. Submitted by: Dianne Oleary (oleary@maryland.arpa) Obtainable from: Dianne Oleary Department of Computer Science University of Maryland College Park, MD 20742 --------------- [50] Dianne P. O'Leary, G. W. Stewart Assignment and Scheduling in Parallel Matrix Factorization Computer Science Department University of Maryland TR-1486, April, 1984. We consider in this paper the problem of factoring a dense n-by-n matrix on a network consisting of P MIMD processors when the net- work is smaller than the number of elements in the matrix 2 (P < n ). The specific example analyzed is a computational net- work that arises in computing the LU, QR, or Cholesky factoriza- tions. We prove that if the nodes of the network are evenly dis- tributed among processors and if computations are scheduled by a round-robin or a least-recently-executed scheduling algorithm, then optimal order of speed-up is achieved. However, such speed- up is not necessarily achieved for other scheduling algorithms or if the computation for the nodes is inappropriately split across processors, and we give examples of these phenomena. Lower bounds on execution time for the algorithm are established for two important node-assignment strategies. Obtainable from: D. Oleary, above. --------------- ##### OAK RIDGE ##### All Oak Ridge reports listed below may be obtained from: Ms. Dawn Human Mathematical Sciences Oak Ridge National Laboratory P. O. Box Y, Bldg. 9207-A Oak Ridge, TN 37831 electronically: human@ornl-msr.arpa ##################### [51] George Ostrouchov ANOVA Model Fitting via Sparse Matrix Computations ORNL Tech. Report, August 1985 The least-squares computational approach is used in fitting analysis of variance models when data are unbalanced or incomplete. For large models containing factors with many levels, standard statistical packages require enormous amounts of time and storage due to the size of the crossproducts matrix. The model matrix X (often called the design matrix) consists of dummy variables and is sparse, thus great reduction in time and storage can be achieved by viewing this as a sparse matrix problem. A method, based on orthogonal Givens factorization of a maximal rank set of sparse columns of X, for computing estimates and residual sums of squares is presented. Selection of a maximal rank set of sparse columns is a key part of this method and is based on the linear dependencies between columns of the model matrix. The linear dependencies are described using a notation based on index sets associated with model terms. Time and storage requirements of the new algorithm are compared to the requirements of GLIM, which uses the same basic numerical algorithm that is used in most statistical packages. Both requirements of the new algorithm are less by up to three orders of magnitude for large models and are competitive for small models. Submitted by: George Ostrouchov -------------- [52] V. Alexiades J.B. Drake G.A. Geist G.E. Giles A.D. Solomon R.F. Wood Mathematical Modeling of Laser-Induced Ultrarapid Melting and Solidification ORNL-6129, July 1985 Mathematical and numerical aspects of laser induced phase change modeling are dicussed. Several numerical formulations, including explicit and implicit enthalpy, adaptive grid finite volume and front tracking are derived and results are presented. In addition, a novel hyperbolic formulation is explored analytically and numerically. The explicit enthalpy approach allows the greatest flexibility and efficiency in modeling problems with complicated phase transitions. Submitted by: Vasili Alexiades -------------- [53] Max Morris Vasili Alexiades Sensitivity Analysis of a Numerically Simulated HgCdTe Solidification Process by Statistical Methods ORNL-6210, August 1985 Statistical response surface methodology is applied to the problem of determining simple approximations for outputs as functions of input parameter values in a numerical simulation of the solidification of a mercury-cadmium-telluride alloy. The approximating polynomials are then differentiated with respect to parameter values to determine sensitivities. A table of percent changes in output values as a result of one per cent changes in parameter values is presented. The empirical procedure used constitutes an unorthodox application of the statistical methodology, but it is easy to use, quite generally applicable and very effective. Submitted by: Vasili Alexiades -------------- [54] V. Alexiades G. A. Geist A. D. Solomon Numerical Simulation of a HgCdTe Solidification Process ORNL-6127, August 1985 The solidification of a cylindrical ingot of mercury-cadmium-telluride is modeled taking into account both heat conduction and solute diffusion. Values of the relevant thermophysical parameters of the pseudo-binary HgTe-CdTe are compiled. The model is implemented numerically by a finite-difference discretization and results of the simulation of a freezing experiment are reported. Submitted by: Vasili Alexiades -------------- [55] George A. Geist Michael T. Heath Parallel Cholesky Factorization on a Hypercube Multiprocessor ORNL-6190, July 1985 Two types of message-passing parallel algorithms are developed for solving symmetric systems of linear equations on a hypercube multiprocessor. One type involves broadcast communication among processors, and the other involves communication along a ring of processors. Details are provided in the form of C programs that implement the algorithms on a hypercube simulator and which should run with little modification on real hypercube hardware. Performance of the various algorithms is demonstrated by means of processor utilization graphs and parallel speedup curves. Submitted by: Mike Heath -------------- [56] Michael T. Heath Parallel Cholesky Factorization in Message-Passing Multiprocessor Environments ORNL-6150, May 1985 Parallel algorithms are presented for computing the Cholesky factorization on multiprocessors having only private local memory. Synchronization of multiple processes is based on message passing. Several possible processor interconnection networks are considered. Submitted by: Mike Heath -------------- [57] George A. Geist ORNL Efficient Parallel LU Factorization with Pivoting on a Hypercube Multiprocessor ORNL-6211, August 1985 A message-passing parallel algorithm is developed for forming the LU factors of general non-singular matrices on a hypercube multiprocessor. Partial pivoting is performed to ensure numerical stability, but the scheduling of tasks is such that the pivot search in the parallel algorithm is completely masked. Empirical tests show that the load imbalance produced by random pivoting causes a 5%-14% increase in execution time. A complementary parallel triangular solution algorithm is also given. Comparisons with the non-pivoting case are used to demonstrate the efficiency of this new algorithm. Submitted by: Al Geist -------------- ##### PITT ##### [58] Rami G. Melhem A Study of Data Interlock in VLSI Computational Networks for Sparse Matrix Multiplication TR-CSD-505 Department of Computer Science Purdue University April 1985 The general question addressed in this study is: Are regular VLSI networks suitable for sparse matrix computations?. More specifically, we consider a special purpose self-timed network that is designed for certain specific dense matrix computation. We add to each cell in the network the capability of recognizing and skipping operations that involve zero operands, and then ask how efficient this resulting network is for sparse matrix computation?. In order to answer this question, it is necessary to study the effect of data interlock on the performance of self-timed networks. For this, the class of pseudo systolic networks is introduced as a hybrid class between systolic and self-timed networks. Networks in this class are easy to analyze, and provide a means for the study of the worst case performance of self-timed networks. The well known concept of computation front is also generalized to include irregular flow of data, and a technique based on the propagation of such computation fronts is suggested for the estimation of the processing time and the communication time of pseudo systolic networks. Submitted by: Rami Melhem (melhem@purdue OR melhem@pitt) Obtainable from: Rami Melhem University of Pittsburgh Department of Mathematics Pittsburgh, PA 15260 --------------- [59] Rami G. Melhem A Modified Frontal Technique Suitable for parallel Systems ICMA-85-84 Department of Mathematics University of Pittsburgh July 1985 Frontal techniques offer the potential for processing the assembly and the factorization phases of finite element analysis in parallel. However, the rows of the stiffness matrix are assembled and factored in different order, thus depriving frontal solvers of the uniformity desirable in parallel processing. On the other hand, band solution techniques handle the factorization phase in a very uniform way but do not interleave assembly and factorization. In this paper, we suggest a technique that borrows from both frontal and band solvers those characteristics that are advantageous for parallel processing. Moreover, book-keeping and data manipulation are simpler in the suggested techniques than in the classical frontal method. This makes the suggested technique also attractive for sequential systems. Obtainable from: Rami Melhem, as above. --------------- ##### STANFORD ##### [60] Gene H. Golub and Carl D. Meyer Simultaneous Computation of Stationary Probabilities with Estimates of Their Sensitivity SIAM ms.no.10958 < No Abstract > Submitted by: Gene H. Golub (golub@su-navajo.arpa) Obtainable from: Gene H. Golub (golub@su-navajo.arpa) Computer Science Department Stanford University Stanford, CA 94305 --------------- [61] Philip E. Gill, Walter Murray, Michael A. Saunders, J. A. Tomlin, and Margaret H. Wright ON PROJECTED NEWTON BARRIER METHODS FOR LINEAR PROGRAMMING AND AN EQUIVALENCE TO KARMARKAR'S PROJECTIVE METHOD We discuss interior-point methods for linear programming derived by applying a logarithmic barrier transformation and performing projected Newton steps for a sequence of barrier parameters. Under certain conditions, one of these ``projected Newton barrier'' methods is shown to be equivalent to Karmarkar's (1984) projective method for linear programming. Details are given of a specific barrier algorithm and its practical implementation. Numerical results are given for several nontrivial test problems. Submitted by: Philip E. Gill (or.gill@su-sierra.arpa) Obtainable from: Philip E. Gill (or.gill@su-sierra.arpa) Systems Optimization Laboratory Operations Research Department Stanford University Stanford, CA 94305 --------------- ##### SUNY/BUFFALO ##### [62] P.J.Eberlein On the Schur Decomposition of a Matrix for Parallel Computation TR 85-689, June 1985, Cornell University An algorithm to solve the eigenproblem for non-symmetric matrices on an N x N array of mesh connected processors, isomorphic to the architecture described by Brent and Luk for symmetric matrices, is presented. This algorithm is a generalization of the classical Jacobi method which holds promise for parallel architectures. Only unitary transformations are used. The rotational parameters are carefully analysed, and many examples of a working program, simulating the parallel architecture, are given. Submitted by: P.J.Eberlein (eberlein%gvax@cornell.arpa) Obtainable from: P.J.Eberlein Dept. of Computer Science SUNY/Buffalo, Bell Hall, Buffalo, N.Y. 14260 --------------- ##### TORONTO ##### [63] D.W. Decker and A. Jepson Convergence cones near bifurcation TR 171/84. < No Abstract > Submitted by: Ken R. Jackson (krj@toronto.csnet) Obtainable from: The Scientific Computing Group, Department of Computer Science, University of Toronto, Toronto, Ontario, Canada M5S 1A4 --------------- [64] Paul Muir Implicit Runge-Kutta methods for two-point boundary value problems TR 175/84. < No Abstract > Obtainable from the above. --------------- [65] Kevin Burrage The order properties of implicit multivalue methods for ordinary differential equations TR 176/84. < No Abstract > Obtainable from the above. --------------- [66] Kevin Burrage Order and stability properties of explicit multivalue methods TR 177/84. < No Abstract > Obtainable from the above. --------------- [67] W.H. Enright, K.R. Jackson, S.P. Norsett, and P.G. Thomsen Interpolants for Runge-Kutta formulas TR 180/85. < No Abstract > Obtainable from the above. --------------- [68] J.M. Fine Low order Runge-Kutta-Nystrom methods with interpolants TR 183/85. < No Abstract > Obtainable from the above. --------------- ##### WASHINGTON STATE ##### [69] ALAN GENZ FULL SYMMETRIC INTERPOLATORY FOR MULTIPLE INTEGRALS CS-85-136 WASHINGTON STATE UNIVERSITY JULY, 1984 A method is given for the direct determination of the weights for fully symmetric integration rules for the hypercube, using multivariable Lagrange interpolation polynomials. The formulas for the weights lead to new classes of efficient rules. Submitted by: Alan Genz ARPANET: na.genz@su-score CSNET: acg@wsu Obtainable from: Alan Genz Computer Science Department Washington State University Pullman, WA 99164-1210 --------------- [70] ALAN GENZ MATRIX METHODS FOR THE FORCED DIFFUSION EQUATION CS-85-137 WASHINGTON STATE UNIVERSITY JULY, 1985 Finite difference methods for the forced diffusion equation with time dependent boundary conditions are developed using matrices to represent both the approximate solution and the difference operations. This technique allows boundary conditions to be included in a natural way and eliminates the need for an analysis of the commutativity of the spatial difference operations. The matrix methods are used to develop algorithms that are second order accurate in space and time for a standard square domain and an annular domain. Obtainable from: Alan Genz, as above. --------------- ##### WATERLOO ##### [71] Senad Busovaca Handling Degeneracy in a Nonlinear L-1 Algorithm Technical Report CS-85-34 Department of Computer Science University of Waterloo This paper is concerned with handling degeneracy in a nonlinear optimization algorithm based on an active set strategy. Although the solution to the problem is given in the context of nonlinear L-1 optimization, the approach is more general and can be used to overcome the non-uniqueness of dual variables in methods that use optimality conditions construc- tively. Submitted by: Richard Bartels Obtainable from: Richard Bartels University of Waterloo Computer Science Department Waterloo, Ontario N2L 3G1 CANADA --------------- .