NA Digest Sunday, April 21, 1991 Volume 91 : Issue 16

Today's Editor: Cleve Moler

Today's Topics:


From: Mike Heath <>
Date: Sun, 14 Apr 91 21:23:36 CDT
Subject: Last Call for SIAG/LA Prize

This is a reminder that the deadline is April 30 for nominations
of papers for the SIAG/LA Prize for the best paper in applicable
linear algebra for the years 1988-1990. Please see previous
announcements for details and eligibility rules. Nominations
should consist of a complete bibliographic citation and a brief
statement justifying the nomination, and should be directed to
the Prize Committee Chairman at the following address:

Michael T. Heath
Center for Supercomputing Research and Development
305 Talbot Laboratory
University of Illinois
104 South Wright Street
Urbana, IL 61801-2932

Phone: 217-244-6915
Fax: 217-244-1351


From: Rob MacLeod <>
Date: Sun, 14 Apr 91 21:28:58 MDT
Subject: Sparse Solvers for Overdetemrined Systems

Hello na-netters,

I am looking for a solver that will work on a slightly overdetermined linear
system which is very sparse and reasonably large (ca. 3024 by 2832). It
arises from a minimization problem with which I can perform inteprolation
of electric potential values over a three-dimensional surface (the human
thorax). I matrix is not symmetric, or shouldn't be in theory, and I have
not yet checked diagonal dominance or positivity.

In the past, I have used direct QR-solvers on the full matrix in problems
of this type which were a bit smaller (1112 by 1048) in very tolerable
amounts of time (6 minutes) on an IBM RS/6000 520. Unfortunately, the
scaled-up version employing this approach for a 3024 by 2832 case ran for
12 hours before I gave up and stopped it. I need something considerably
faster that will run as a two-step process so that I can either
backsubstitute numerous right hand sides back into the decomposed matrix,
or do an iterative computation that will run in a few minutes.

Any suggestions out there?



Rob MacLeod, Ph.D.
Nora Eccles Harrison Cardiovascular Research and Training Institute (CVRTI)
Building 500
University of Utah
Salt Lake City
Utah 84112


From: Tilak Ratnanather <>
Date: Mon, 15 Apr 91 12:28:59 +0100
Subject: Block Tridiagonal Matrices

The following problem arising in magnetism was brought to
my attention by colleagues who'd greatly appreciate input
from experienced numerical analysts:

First question: is there a program that will compute the
eigenvalues of a BLOCK TRIDIAGONAL (Hermitian) MATRIX. Attempts
to search for one in netlib failed to yield a program. (Fortran
or even C is the preferred language).

Second question: at a more deeper level, the block tridiagonal
Hermitian matrix M may be written in the shorthand form:
(B*, A, B) where A and B are 5x5 matrices and * indicates the
complex conjugate and M is 100x100. Is there a way of reducing this
"large" problem to a sequence of reduced problems (or 5x5) problems?
One step that has been suggested is to make use of the diagonalised
form of A. The entries of A and B may be changed pending the physical

Any ideas or pointers to references would be greatly appreciated. It
may be possible to summarise the responses at a later date.

Thanks a lot

Tilak Ratnanather
Dept of Mathematics
City University
Northampton Square
London EC1V 0HB



From: Roland England <>
Date: Tue, 16 APR 91 14:19:12 GMT
Subject: Shampine/Gordon Integrator in C

In connection with Kathryn Brenan's enquiry, we were also about to start
rewriting the Shampine/Gordon integration software (Adams numerical
integration code for ordinary differential equations) in C.

We would greatly welcome any available information, if this has already
been done elsewhere, and particularly, of course, if there is
a public domain version in existence which we could use.

Thanks for any replies.

Roland England

Tel. +44-908-65-2329

The Open University
Faculty of Mathematics
Milton Keynes
Great Britain


From: Steve Nichols <>
Date: 20 Apr 91 13:57:35 GMT
Subject: C Source for a Stiff ODE Integrator

I'm looking for C source code for a stiff ODE integrator. If anyone
has converted LSODE (from the netlib) to C, this would be perfect.
Thanks for any help,

Steve Nichols
Georgia Tech Physics Department


From: Nick Higham <>
Date: Sat Apr 13 14:31:37 PDT 1991
Subject: Three Measures of Precision in Floating Point Arithmetic

Three measures of precision in floating point arithmetic
by Nick Higham

This note is about three quantities that relate to the precision of
floating point arithmetic. For t-digit, rounded base b arithmetic the
quantities are

(1) machine epsilon (eps), defined as the distance from 1.0 to
the smallest floating point number bigger than 1.0 (and
given by eps = b**(1-t), which is the spacing of the
floating point numbers between 1.0 and b),

(2) mu = smallest floating point number x such that fl(1 + x) > 1, and

(3) unit roundoff u = b**(1-t)/2 (which is a bound for the relative
error in rounding a real number to floating point form).

The terminology I have used is not an accepted standard; for example,
the name machine epsilon is sometimes given to the quantity in (2).
My definition of unit roundoff is as in Golub and Van Loan's book
`Matrix Computations' [1] and is widely used. I chose the notation eps
in (1) because it conforms with MATLAB, in which the permanent
variable eps is the machine epsilon. [Ed. note: Well, not quite.
See my comments below. --Cleve]

The purpose of this note is to point out that it is not necessarily the
case that mu = eps, or that mu = u, as is sometimes claimed in the
literature, and that, moreover, the precise value of mu is difficult to

It is helpful to consider binary arithmetic with t = 3. Using binary
notation we have

1 + u = 1.00 + .001 = 1.001,

which is exactly half way between the adjacent floating point numbers
1.00 and 1.01. Thus fl(1 + u) = 1.01 if we round away from zero when
there is a tie, while fl(1 + u) = 1.00 if we round to an even last
digit on a tie. It follows that mu <= u with round away from zero (and
it is easy to see that mu = u), whereas mu > u for round to even.
I believe that round away from zero used to be the more common choice
in computer arithmetic, and this may explain why some authors define
or characterize u as in (2). However, the widely used IEEE standard
754 binary arithmetic uses round to even.

So far, then it is clear that the way in which ties are resolved
in rounding affects the value of mu. Let us now try to determine the
value of mu with round to even. A little thought may lead one to
suspect that mu <= u(1+eps). For in the b = 2, t = 3 case we have

x = u*(1+eps) = .001*(1+.01) = .00101

=> fl(1 + x) = fl( 1.00101 ) = 1.01,

assuming ``perfect rounding''. I reasoned this way, and decided to
check this putative value of mu in 386-MATLAB on my PC.
MATLAB uses IEEE standard 754 binary arithmetic, which has t = 53
(taking into account the implicit leading bit of 1). Here is what I

>> format compact; format hex
>> x = 2^(-53)*(1+2^(-52)); y = [1+x 1 x]
y =
3ff0000000000000 3ff0000000000000 3ca0000000000001

>> x = 2^(-53)*(1+2^(-11)); y = [1+x 1 x]
y =
3ff0000000000000 3ff0000000000000 3ca0020000000000

>> x = 2^(-53)*(1+2^(-10)); y = [1+x 1 x]
y =
3ff0000000000001 3ff0000000000000 3ca0040000000000

Thus the guess is wrong, and it appears that mu = u*(1+2^(42)*eps)
in this environment! What is the explanation?

The answer is that we are seeing the effect of ``double-rounding'', a
phenomenon that I learned about from an article by Cleve Moler [2].
The Intel floating-point chips used on PCs implement internally the
optional extended precision arithmetic described in the IEEE standard,
with 64 bits in the mantissa [3]. What appears to be happening in the
example above is that `1+x' is first rounded to 64 bits; if
x = u*(1+2^(-i)) and i > 10 then the least significant bit is lost in
this rounding. The extended precision number is now rounded to 53 bit
precision; but when i > 10 there is a rounding tie (since we have
lost the original least significant bit) which is resolved to 1.0,
which has an even last bit.

The interesting fact, then, is that the value of mu can vary even
between machines that implement IEEE standard arithmetic.

Finally, I'd like to stress an important point that I learned from the
work of Vel Kahan: the relative error in addition and subtraction
is not necessarily bounded by u. Indeed on machines such as Crays
that lack a guard digit this relative error can be as large as 1. For
example, if b = 2 and t = 3, then subtracting from 1.0 the next
smaller floating number we have

Exactly: 1.00-

Computed, without 1.00-
a guard digit: .11 The least significant bit is dropped.

The computed answer is too big by a factor 2 and so has relative
error 1! According to Vel Kahan, the example I have given mimics what
happens on a Cray X-MP or Y-MP, but the Cray 2 behaves differently and
produces the answer zero. Although the relative error in
addition/subtraction is not bounded by the unit roundoff u for
machines without a guard digit, it is nevertheless true that

fl(a + b) = a(1+e) + b(1+f),

where e and f are bounded in magnitude by u.

[1] G. H. Golub and C. F. Van Loan, Matrix Computations, Second Edition,
Johns Hopkins Press, Baltimore, 1989.

[2] C. B. Moler, Technical note: Double-rounding and implications for
numeric computations, The MathWorks Newsletter, Vol 4, No. 1 (1990), p. 6.

[3] R. Startz, 8087/80287/80387 for the IBM PC & Compatibles, Third
Edition, Brady, New York, 1988.


Editor's addendum:

I agree with everything Nick has to say, and have a few more comments.

MATLAB on a PC has IEEE floating point with extended precision implemented
in an Intel chip. The C compiler generates code with double rounding.
MATLAB on a Sun Sparc also has IEEE floating point with extended precision,
but it is implemented in a Sparc chip. The C compiler generates code
which avoids double rounding.

On both the PC and the Sparc

eps = 2^(-52) = 3cb0000000000000 = 2.220446049250313e-16

However, on the PC

mu = 2^(-53)*(1+2^(-10)) = 3ca0040000000000 = 1.111307226797642e-16

While on the Sparc

mu = 2^(-53)*(1+2^(-52)) = 3ca0000000000001 = 1.110223024625157e-16

Note that mu is not 2 raised to a negative integer power.

MATLAB on a VAX usually uses "D" floating point (there is also a "G"
version under VMS). Compared to IEEE floating point, the D format has
3 more bits in the fraction and 3 less bits in the exponent. So
eps should be 2^(-55), but MATLAB says eps is 2^(-56). It is actually
using the 1+x > 1 trick to compute what we're now calling mu. There
is no extended precision or double rounding and ties between two floating
point values are chopped, so we can find mu by just trying powers of 2.

On the VAX with D float

eps = 2^(-55) = 2.775557561562891e-17

mu = 2^(-56) = 1.387778780781446e-17

The definition of "eps" as the distance from 1.0 to the next floating
point number is a purely "geometric" quantity depending only on the
structure of the floating point numbers. The point Nick is making is
that the more common definition of what we here call mu involves a
comparison between 1.0 + x and 1.0 and subtle rounding properties of
floating point addition. I now much prefer the simple geometric
definition, even I've been as responsible as anybody for the popularity
of the definition involving addition.

-- Cleve


From: Milo Dorr <>
Date: Sun, 14 Apr 91 15:34:53 PDT
Subject: Meeting on Mathematics, Computations, and Reactor Physics

The American Nuclear Society (Mathematics & Computation Division and
Reactor Physics Division) International Topical Meeting:


April 28 - May 2, 1991
Green Tree Marriott
Pittsburgh, PA

The banquet speaker on Wednesday, May 1 will be George E. Lindamood of
the Gartner Group. The title of his lecture will be "Why
Supercomputing Matters: A Perspective of the Proposed Federal High
Performance Computing and Communication Program".

There will also be a plenary session at 9:30am on Monday, April 29,
consisting of a panel discussion on the topic "Perspectives on
Advances in Supercomputing Performance" moderated by James. R.
Kasdorf from Westinghouse Corporate Computer Services and the
Pittsburgh Supercomputing Center. The panelists and the titles of
their prepared presentations are:

Gregory J. McRae (Carnegie Mellon University)
"Grand Challenges in Computational Science"

W. B. Barker (BBN Advanced Computers, Inc.)
"Parallel Computing: Past, Present, and Future"

M. L. Barton (Intel Supercomputing Systems Division)
"Technology Development for Supercomputing 2000"

Kenichi Miura (Fujitsu America, Inc.)
"Perspectives on Advances in Supercomputer Performance
- Fujitsu's View"

Steve Nelson (Cray Research, Inc.)
"Heterogeneous Supercomputing: Now and Then"

For a copy of the Technical Program or further information, please
contact the Program Chairman:

I. K. Abu-Shumays
RT-Mathematics, 34F
Bettis Atomic Power Laboratory
P. O. Box 79
West Mifflin, PA 15122-0079
(412) 476-6469, FAX (412) 476-5151


From: Biswa Datta <
Date: Mon, 15 Apr 91 22:40:59 CDT
Subject: Second NIU Conference

Second NIU Conference on Linear Algebra, Numerical Linear,
Algebra and Applications
May 3 - 5, 1991

Holmes Student Center
Northern Illinois University

Organizer and Chairman - Biswa Datta, Northern Illinois University
Advisor: Hans Schneider, University of Wisconsin - Madison

Sponsored By
The Institute for Mathematics and Its Applications (Minnesota)
and International Linear Algebra Society (ILAS)
Northern Illinois University

The purpose of the conference is to bring together researchers in linear
algebra, numerical linear algebra, and those working in various application
areas for an effective exchange of ideas and discussion of recent
developments and future directions of research.

All events will take place at the Holmes Student Center, Northern Illinois
University. Registration with begin at 7pm on Thursday, May 2. The
meeting itself with be Friday, May 3, through Sunday May 5.

Invited speakers include:

R. Thompson, University of California-Santa Barbara
R. Plemmons, Wake Forest University
Clyde Martin, Texas Tech University
Roger Horn, The Johns Hopkins University
Charles R. Johnson, College of William and Mary
Daniel Hershkowitz, Tecnion-Israel Institute of Technology
Robert Grossman, University of Illinois at Chicago
James R. Bunch, University of California-San Diego
Floyd Hanson, University of Illinois at Chicago
Chris Bischof, Argonne National Laboratory
William Gragg, Naval Post Graduate School at Monterey, California
Lothar Reichel, Naval Postgraduate School
Richard Brualdi, University of Wisconsin - Madison
Thomas Laffey, University College, Dublin, Ireland
Jose Dias de Silva, University of Lisbon, Lisbon, Portugal
S. Campbell, North Carolina State University
Kenneth Clark, Army Research Office
William Hager, University of Florida
George Cybenko, University of Illinois at Urbana-Champaign
Patricia Eberlein, University of Buffalo/SUNY
Kermit Sigmon, University of Florida
Daniel Boley, University of Minnesota
Homer Walker, Utah State University
Roland Freund, RIACS, NASA AMES Research Center
A. Yeremin, USSR Academy of Sciences
David Young, University of Texas, Austin
Michael Neumann, University of Connecticut
Avi Berman, Technion-Israel Institute of Technology
Chris Byrnes, Washington University
S. P. Bhattacharyya, Texas A & M University
Bijoy Ghosh, Washington University
Rama K. Yedavalli, Ohio State University
S. Friedland University of Illinois at Chicago
S. Wright, Argonne National Laboratory
Pradip Misra, Wright State University
Mohsen Pouramdi, Northern Illinois University

For more information about the program, contact faculty coordinator Biswa
Nath Datta, Department of Mathematical Sciences, at (815) 753-6759. For more
information about the logistics, contact Margaret Shaw, College of
Continuing Education, at (815) 753-1458.


From: J. C. Butcher <>
Date: Wed, 17 Apr 91 10:13:28 NZS
Subject: Position at the University of Auckland

Applied and Computational Mathematics Unit

The University of Auckland invites applications for a lectureship
in the Applied and Computational Mathematics Unit within the
Department of Mathematics and Statistics. Applicants should have
a proven record in teaching and research in some branch of
Applied or Computational Mathematics. Applications from
candidates with expertise in a field that will strengthen and
enhance the existing research interests of the Applied and
Computational Mathematics Unit in differential equations and
their numerical solutions and applications of scientific computing
to physical and other problems are particularly welcome.

The University of Auckland has a student population of about
17,000 and is the only university in Auckland, a city with about
one million inhabitants. Mathematics and Statistics is the largest

University positions are organised in a system similar to that of
the UK. Thus a lectureship is intended to be a permanent
appointment although, in the first instance, it is for a four year
term. Most appointments are continued after that period.

The salary scale for a lecturer begins at $NZ37,440 and increases
in annual increments. (Details of the salary scale and prospects
for promotion are available on request).

It is hoped that the successful applicant will take up his or her
duties by 1 August 1991 or soon after.

Applications, in accordance with the "Method of Application"
(available on request) should be submitted as soon as possible and
no later than 30 June 1991.

The University of Auckland is an Equal Opportunity Employer

For further information or enquiries about the position please
contact the Head of the Applied and Computational Mathematics Unit
Professor J. C. Butcher using email. or


From: Alan Craig <>
Date: Thu, 18 Apr 91 16:35:46 BST
Subject: Postgraduate Opportunity at University of Durham

The Department of Mathematical Sciences, University of Durham and the British
Gas Research Station in Northumberland have a CASE studentship in

Adaptive Finite Element Analysis
for Nonlinear Elastic Problems

The award is for three years, tenable from October 1991 and for a programme of
research leading to a PhD.

Outline of project: the finite element method is the computational tool for the
solution of structural analysis problems. The technique relies on replacing the
underlying differential equations by an approximating set of equations. However
the correct choice of this approximating set is by no means a trivial procedure.

In recent years adaptive methods have been developed which attempt to choose the
approximation in an automatic way linked to the problem. The analysis and
implementation of these methods for linear problems is now well understood. The
situation for nonlinear systems is however, less satisfactory and the interest
in adaptive techniques for general nonlinear response stems from the need to
maintain structural integrity in hazardous situations. In particular British
Gas are interested in techniques which can be used to analyse their offshore

The project will maintain a careful balance between theoretical analysis and
practical implementation. Experience has shown that in this area, as in many
others, the implementation can be guided and enriched by the analysis. The
ultimate aim is to produce high quality algorithms for nonlinear problems.

As part of the project the successful applicant will work in the British Gas
Research Station for periods totaling three months under the supervision of the
Industrial Supervisor Mrs. Jane Haswell.

The applicant will have, or will be about to obtain, a good honours degree in
Mathematics, Engineering or a related discipline. In addition to the standard
SERC award British Gas will make an extra payment of 1620 pounds per annum to
the student.

Further information can be obtained from Dr. Alan Craig at the address below and
informal enquiries are welcomed. Application should be made to the same address
and should include the names and addresses of at least two referees, one of whom
should be able to judge the applicants suitability for research.

Department of Mathematical Sciences
University of Durham
South Road
Durham DH1 3LE


From: Antonio Pizarro
Date: Fri, 19 Apr 91 08:35:56 IST
Subject: Position at Centenary College of Louisiana

Applications are invited for two tenure-track positions in MATH. A
Ph.D in Mathematics is required. Linear Algebra, Topology or Applied
Math are preferred fields, but all fields are considered. Beginning
Fall 1991. Send resume and three letters of recommendation to:

Antonio G. Pizarro, Chair
Math Department
Centenary College of Louisiana
Shreveport LA 71134, USA.


From: Anthony Skjellum <>
Date: Fri, 19 Apr 91 11:15:54 PDT
Subject: Position at Lawrence Livermore Laboratory

Postdoctoral Position in Parallel Computation - Lawrence Livermore

A postdoctoral position in scientific parallel algorithms research is
presently available in the Numerical Mathematics Group (NMG) at the
Lawrence Livermore National Laboratory, Computing and Mathematics
Researh Division. Candidates should be well versed in parallel
computation; particularly, distributed memory machines and issues.
Experience with large-scale parallel scientific algorithms and
applications is strongly desired. In-depth experience with the Unix
operating system and the C/C++ languages is also desired. Candidate
will participate in the growing NMG research effort in parallel
computation, but also will enjoy freedom to explore his or her own
allied research interests. A one year position is offered, renewable
for a second year if both NMG and the fellow concur. Candidate should
have completed his or her Ph.D. in an appropriate discipline prior to
October 1, 1991. Starting date is on or after October 1, 1991. US
citizenship is required. Send three letters of recommendation, resume,
and publications list to Dr. Anthony Skjellum, LLNL L-316, PO Box 808,
Livermore, CA 94550. (415)422-1161. e-mail:


End of NA Digest


From surfer.EPM.ORNL.GOV!nacomb Sat Apr 27 22:36:25 0400 1991