NA Digest Saturday, August 24, 1991 Volume 91 : Issue 34

Today's Editor: Cleve Moler

Today's Topics:

Submissions for NA Digest:

Mail to

Information about NA-NET:

Mail to


From: Herman J. Woltring <UGDIST@NICI.KUN.NL>
Date: Mon, 19 Aug 91 07:59 MET
Subject: Cubic & Quartic Equation Solvers

Dear NA-NET readers,

I should be grateful for FORTRAN code to solve (cubic and) quartic equations.
It seems that the algorithm in section 3.8.2 in my edition of Abramovitz &
Stegun is incorrect, while I am getting good results with the algorithm in
section 5.5 of "Numerical Recipes" by Press et al. However, rounding errors
in the quartic equation algorithm of section 3.8.3 in Abramovitz and Stegun
make it difficult to take reliable branching decisions.

Thanks in advance!

Herman J. Woltring <na.woltring>, Eindhoven, The Netherlands


From: John Richardson <>
Date: 23 Aug 91 15:36:48 GMT
Subject: GMRES Fortran Code

Does anyone have, or know of an ftp site that has, f77 or other source code
the implements the Generalized Minimum RESidual (any flavor) algorithm?


John Richardson
Thinking Machines Corporation, Cambridge Mass., USA


From: Glen Clark <...!gatech!dscatl!opto!glen>
Date: 23 Aug 91 13:29:00 GMT
Subject: UGLY Optimization Problems Wanted

Over the course of several years, we have written a
suite of optimizers for constrained, non-linear
problems. The original application was electromag-
netics, but the core algorithms are not application-
specific and we believe they may be useful for a
broad spectrum of non-linear problems.

Immunity to cheap, early gains is high and execution
time compared to SA is very favorable.

We'd like to find out just how good these algorithms are.
Every mother thinks her child is beautiful, but there's
nothing like a little bit of competition to give you
a quick reality check. To that end, we're soliciting UGLY,
non-linear, optimization problems (parametric, not
combinatorial) in an effort to stump the system.

If you have a particularly challenging optimization
problem that you'd like to see humbled, we'd like to see
if we can break our pick on it, at our expense. If
interested, please send us a paragraph or two giving a
general outline of the problem. Minimally, it should
include subject area (aerodynamics, network loading,
etc.) and a brief description of the particular application.
Any information regarding the present state of the art in
that application would be helpful, but is not essential.
All replies will be acknowledged.

C or FORTRAN problems are preferred, but interesting
problems already coded in ada or pascal will also be

All exchanges will be held strictly confidential.
Binding non-disclosure agreements available where

Glen Clark, P.E. gatech!dscatl!opto!glen

Glen Clark & Associates
Consulting Engineers
Alpharetta, GA


Date: Tue, 20 Aug 91 11:57 EDT
Subject: SIAM Conference in Minneapolis

Dear Colleague:

Just a friendly reminder.....if you have not yet registered for
the conference on applied linear algebra, September 11-14, 1991,
Minneapolis, Minnesota, please take a few minutes of your valuable
time ....and register now -- before the deadline -- September 4.
You can register by telephone -- 215-382-9800, by fax -
- 215-386-7999, or by e-mail to:

Thanks and see you in Minneapolis.

SIAM Conferences


From: John Reid <>
Date: Tue, 20 Aug 91 10:19:28 BST
Subject: Does IEEE Arithmetic Stink?

This is a response to Joe Grcar's article in NA Digest 91/33 (Aug 18). I
think that he is attacking the wrong target. It is not IEEE arithmetic that
stinks. It is the implementation of Fortran 77 that stinks. The Standard
requires (page 15-8, lines 8-11 and page 15-10, lines 36-39) that the actual
arguments in a procedure reference must agree in order, number, and type with
the corresponding dummy arguments of the referenced procedure. Most
implementations have no mechanism for checking this condition. Everyone
makes silly mistakes from time to time and it is highly desirable that
as many as possible are properly diagnosed, preferably before execution
commences. As far as I can tell from Joe's description, his engineer's faulty
program could have been diagnosed at link time.
The situation will be better in Fortran 90. There are simple mechanisms for
making the properties of the dummy arguments visible at compile time and I am
expecting that libraries will make these available to their users by providing
modules instead of procedures. An example is a module that contains nothing
more than a single Fortran 77 procedure; the only change that a Fortran 77
programmer would need to make is to add a USE statement that names the module
and thereby makes the interface visible to the compiler.
John Reid,
Atlas Centre,
Rutherford Appleton Laboratory,
Didcot, Oxon OX11 0QX


Date: Mon, 19 Aug 91 11:13:49 EDT
Subject: IEEE Arithmetic Stinks Less Than ...

I trust that others who were involved in the IEEE development
and approval process will respond to Joe Grcar's complaints at
a technical level. The purpose of this note is to give the
viewpoint of someone who -- though a minor participant in IEEE
committees -- is mainly a user of their product in developing
modest scientific computer programs. Grcar's points are taken
as he presented them.

(1) The standard makes simple things complicated.

"In contrast, the SUN OS documentation for IEEE arithmetic is
one hundred and sixty pages long." Surely this is a complaint
about SUN documentation. The IEEE 754 standard is one of the
shortest standards documents in the computer arena. Moreover,
it DOES make things simple(r) for the user, because the
properties of the arithmetic are reproducible. I recall Data
General changing the radix from 2 to 16 between two releases of

(2) The standard makes precision indefinite.

Implementers (of compilers) should give the user a fairly
sensible set of defaults e.g. underflow gracefully to zero,
generate error exception on overflow, etc. Grcar's complaint
here is one that such defaults are not provided.

(3) The standard makes error analysis impossible.

Someone else can field this one. (see (2) above).

(4) The standard doesn't prevent failure, it only hides

The programmer is not absolved of the responsibility of
checking the data against the real world. The IEEE arithmetic
is a set of tools to keep the alligators down to waist level.

(5) The standard makes debugging difficult.

"The standard allows programs to run with garbage for data.
These are the hardest programs to debug." If you (or the
compiler writer) turn off the traps (see (2) above), like
the operators at Chernobyl or Bhopal turned off the safety
devices, you can expect to have accidents. Try driving
along a mountain road without steering wheel or brakes.

Grcar's example: This example is akin to a family disaster
which occurred when the cook made zabaglione with salt instead
of fine sugar. If you attach the wrong libraries to a linker,
you will have trouble, but it is not the fault of the
arithmetic. In this regard, I like programs which come with a
modest test output one can try to reproduce before using the
program seriously.

In my experience on PCs, compilers have become a whole lot
nicer to use since the writers adopted IEEE as an option. More
recently it is often the default, and for those who don't have
the 80x87 or a comparable co-processor, there are software
emulations which are usually invoked automatically. I am
tempted to infer that Joe Grcar has been the victim of an
implementation which is user-UNFRIENDLY.

A useful exercise, if one is worried about the arithmetic or
its default settings, is to run PARANOIA (available on NETLIB
courtesy of many workers, in particular Kahan and Karpinski).
Before responding here, I tried recently acquired Microsoft
Quick C (1 flaw, 1 defect reported), and Turbo C++ (gave
compile errors, which I will have to look into).

John Nash, Faculty of Administration, University of Ottawa


From: F. W. Olver <>
Date: Thu, 22 Aug 91 16:24:57 EDT
Subject: IEEE Arithmetic

Joe Grcar's severe criticisms of IEEE arithmetic in his message of Aug.14
may well touch resonant chords in the computing community. We, too,
would be very interested in hearing details about similar complaints.

The problem lies not so much with the IEEE standard as with the floating-
point system itself. This system was widely adopted in the 1950's when the
problem of rescaling to keep within the representable ranges of the
fixed-point system became irksome AND when suitable hardware was constructed
by the engineers. Although the floating-point system has proved to be a
vast improvement on the fixed-point system, nevertheless it did not solve
the rescaling problem completely. All too frequently programmers and
software designers encounter the phenomena of overflow and underflow---
the "betes noires of portability" as W.S. Brown once dubbed them.

The situation can be even worse than Joe Grcar depicts. In a problem
encountered this summer at NIST (formerly NBS) concerning turbulent
combustion, IEEE arithmetic yielded not the recognizable, but useless,
NaN's or infinities, but quite reasonable-looking---and seriously erroneous---
graphs without an error message or other warning. The cause had nothing
to do with the standard graphics packages being used. It was eventually
traced to floating-point underflow. Although IEEE arithmetic was used it
would have occurred with any currently available floating-point hardware.
Neither the touted IEEE gradual underflow nor the use of double precision
made any difference. (And even when the erroneous nature of the output
was recognised it was clear that no simple rescaling would remedy the

In recent years a new form of computer arithmetic has been under development
as an alternative to the floating-point system. The new system is called
the symmetric level-index (sli) system, and it has many advantages. First,
the sli system is closed (in a natural way, not with artificial machine symbols)
under arithmetic operations in finite-precision arithmetic, other than
division by zero, of course. This abolishes the problems of overflow and
underflow completely. Next, because the representation is mathematically
determined it is natural and its description is brief---no gigantic committee
was needed to create a specification. Thirdly, there is a natural error
measure that is uniform in nature for all numbers.This error measure
reduces to absolute precision or relative precision in appropriate
circumstances, and is therefore called generalised precision. In short, all
five forms of woe listed in Joe Grcar's message would be taken care of
satisfactorily in the sli system.

Predictably (because so much is now at stake), the big guns of the defenders
of floating-point and the IEEE standard have already been trained on the
new system. The odd thing, however, is that the arguments they are using---
slowness of execution speed, loss of precision (in their own carefully
defined sense), and lack of requisite error analysis---are exactly the same
as those used forty years ago by fixed-point diehards in attacking the
then proposed floating-point system.

Several papers and articles have already been published on the sli system
and its companion, the level-index (li) system. A brief introduction
appeared in the April issue of the Notices of the American Mathematical
Society. For an indepth account see pp.95-168 of Lecture Notes in
Mathematics No.1397 "Numerical Analysis and Parallel Processing", edited by
P.R.Turner, Springer-Verlag, 1989. Software implementations are also available.

Daniel Lozier (, Frank Olver (,
and Peter Turner(


From: David K. Kahaner <>
Date: Mon, 19 Aug 91 15:44:44 EDT
Subject: Japanese Receive US Researchers

ABSTRACT. Summary of NSF "Directory of Japanese Company Laboratories
Willing to Receive American Researchers".

The National Science Foundation sponsors an active Japan program for US
scientists. As part of that program they have prepared a directory of
approximately 150 Japanese companies that have said they are willing to
receive American researchers at their laboratories. The list was compiled
by mailing to 553 companies with more than 30 researchers listed in the
"All Japan Research Institutes Directory". Of these, 284 companies
responded and 154 provided positive responses. NSF's Directory lists each
company as well as specific information about their activities,
personnel, facilities, and the type of research that they are interested
in supporting. Names and addresses of contacts are also given for each.
This Directory, as well as other information can be obtained from the
following address,
Japan Programs
Division of International Programs
National Science Foundation
Washington DC 20550
or by fax at (202) 357-5839.


From: David Walker <walker@rios2.EPM.ORNL.GOV>
Date: Tue, 20 Aug 91 08:10:04 -0400
Subject: Postdoc Position at ORNL




A postdoctoral research position is available at Oak Ridge National
Laboratory in parallel algorithms for linear algebra. The research will
focus initially on the development of a library of portable, efficient
routines for performing dense linear algebra computations on MIMD distributed
memory concurrent computers. A subsequent goal will be the development of
software for sparse matrix computations on these types of computer.
The practical use of the library software will be demonstrated by using it
to build large-scale applications that are transparently portable between
different concurrent computers.

The research to be conducted is part of a program involving Jack Dongarra
and David Walker of Oak Ridge National Laboratory, Jack Demmel of University
of California at Berkeley, Michael Heath of University of Illinois, and
Danny Sorensen of Rice University.

The position requires an in-depth knowledge of linear algebra algorithms
for scientific computation, and practical experience with MIMD distributed
memory concurrent computers. The appointment will be for two years, preferably
beginning October 1, 1991, or shortly thereafter. Benefits of the position
include a competitive salary, travel opportunities, access to state-of-the-art
computational facilities, and collaborative research opportunities in a very
active research program.

Inquiries should be directed to:

David W. Walker
Bldg. 6012 / MS-6367
P. O. Box 2008
Oak Ridge Nat


From: Frank Luk <luk@jacobi.EE.CORNELL.EDU>
Date: Tue, 20 Aug 91 16:56:31 -0400
Subject: Workshop on Total Least Squares

A Workshop on Total Least Squares was held on August 12-15
in Leuven, Belgium, on the campus of the Catholic University.
The two major organizing figures were Sabine Van Huffel and Gene Golub.
There were 60 participants from many countries, from as far away as China.
My three hour layovers at the JFK and Syracuse airports became insignificant
after I learned from the Romanian speaker about her 36 hour sit-up train
ride each way.

Two distinguishing features of this workshop were the hospitality
of the locals (all visitors were picked up at the Brussels airport
on arrival, and returned there on departure), and a mix of participants
from wide ranging areas such as chemistry, geophysics, medical imaging,
scientific computing, signal processing, and statistics.

The list of participants was studded with stars. To NA-NETers,
three familiar names among many would be Bjorck, Golub, and Stewart.
Participants from other disciplines included leaders of similar lofty statures.

This workshop has achieved its goal of promoting scientific exchanges.
I, for one, was quite ignorant of the many important practical applications
of total least squares. It was gratifying to observe how several speakers
tried to make their talks comprehensible to all by using notations
from another discipline. Laughter roared when one speaker described
how confused he quickly became, when he tried to use unfamiliar notations
to prepare his slides; he finally decided that it would be safer for him
to confuse the audience instead.

A nontypical presentation was made by Richard Roy, who heads a small company
in California. He revealed that, unlike other participants, he was there
to try to make money out of the ideas presented, and he referred
to the $VD, whereupon Golub asked if he had used a spell checker.

Good food was abundant: lunches in a student cafeteria, and dinners
in a restaurant. The renowned Belgian beer was always a topic of discussion.
Some were amazed that a 9% beer could be selected as a lunch beverage
in the student cafeteria.

Everyone appreciated the huge amount of work put in by Van Huffel.
particularly since she was due to become a mother five days
after the workshop. Nonetheless, she gave two excellent talks,
and attended other presentations and all social activities.
Lastly, I would express my envy of the achievements of Joos Vandewalle
at the Catholic University: not only has he built up
a strong research group of almost ten students in matrix computing,
he has also created an atmosphere permeated with friendliness
and collegiality. His work was leadership by example at its best:
he was probably the busiest driver shuttling between Leuven
and the Brussels airport on the arrival and departure dates.

Epilogue: Sabine gave birth to a boy on August 19.


From: Mike Bieterman <>
Date: Thu, 22 Aug 91 13:42:23 PDT
Subject: Pacific N.W. Numerical Analysis Seminar

This is to announce the fifth annual

SEPTEMBER 21, 1991




Ian A. Cavers, University of British Columbia,
- chosen as "best student abstract" -
``Complete Spectrum Determination for Symmetric
Sparse Matrices"

Adel Faridani, Oregon State University,
``Local Tomography"

Tom Grandine, Boeing Computer Services,
``Some Numerical Computations Arising in CAD"

Robert O'Malley, University of Washington,
``Asymptotic/Numerical Methods for
Differential-Algebraic Equations"

Brian Tomas, University of Washington \& Boeing Aerospace,
``An Introduction to Wavelets for Computation"

Brian Wetton, University of British Columbia,
``Finite Difference Methods for Navier Stokes
with Boundaries: Vorticity and Primitive Variable Methods"


Cleve Ashcraft, Boeing Computer Services,
``On formulating general sparse and multifrontal as extremal
members of one algorithm family".

John T. Betts and Paul D. Frank, Boeing Computer Services,
``Sparse Nonlinear Programming and its Application to
Trajectory Optimization".

Forrester Johnson, Robin Melvin, Dave Young, Dave Foutch,
John Bussoletti, and Mike Bieterman, The Boeing Company,
``Propulsion System Exhaust Simulation Using Full Potential
and Euler Formulations".

Ping Lin, University of British Columbia,
``Remarks on a Mildly Nonlinear Turning Point Problem".

Dan Pierce and John G. Lewis, Boeing Computer Services,
``Sparse Rank Revealing QR Factorization".

Andre Weideman, Oregon State University,
``Singualrities in an Integrable Semi-discretization
of the Nonlinear Schrodinger Equation".

Laurence B. Wigton, Boeing Commercial Airplane Group,
``High Quality Grid Generation Using Laplacian Sweeps".

Dan Curtis and Elizabeth Kirby, Boeing Computer Services,
``Some Examples of Scientific Visualization at Boeing".

Additional poster presentations are welcome. If you would like to
contribute one, please contact the organizers to assure space.

In addition to the talks and poster session, the day will include
a catered dinner at Jim Phillips' home. The dinner will feature Northwest
delicacies (poached salmon among others) and microbrewery beer. The cost of
the dinner will be $15 for adults and $10 for children. Please indicate on
your registration if you will be attending the dinner and make checks for the
dinner payable to Roger Grimes.

During the seminar, coffee, drinks and lunch will be provided by
The Boeing Company.


Register by September 16. We need adequate time for Boeing security to
issue visitor's passes. Please be aware that late registrations will cause
us some inconvenience. You may register for the seminar by e-mail or
or regular mail. Please, in your mail, indicate how many people will be
attending the dinner. Make checks for the dinner payable to Roger Grimes.
(You may pay either when you send in your registration or at the seminar.)
If you do not plan to attend the dinner, there will be no fees. Note:
Persons under eighteen will not be permitted to attend the seminar.


Boeing Computer Services
33-07 Building
2760 160th Avenue South East
Bellevue, Washington

A map will be sent by mail. If you receive this announcement via electronic
mail and do not receive a copy (with included map) via regular mail within
one week's time, you can request a regular mail copy by contacting one of the

If you have any questions or concerns, please do not hesitate to contact a
member of the organizing committee. They are:

Mike Bieterman - (206) 865-6272 -
Dave Ferguson - (206) 865-3521 -
Roger Grimes - (206) 865-3585 -
John Lewis - (206) 865-3510 -

Their common mailing address is:

Boeing Computer Services
Attn: --------- MS 7L-21
P. O. Box 24346
Seattle, Washington 98124-0346

See you on September 21.


From: Eric Grosse <>
Date: Sat, 24 Aug 91 11:42 EDT
Subject: Netlib News: Big Files

\title{Netlib news: big files}
\author{Eric Grosse}{}

In the early days of netlib, most files were of modest size. A typical
request, say for the ``real symmetric'' driver from eispack, produced a
reply of 30 kilobytes or less. The occasional Fortran program of more
than 100 kilobytes was split arbitrarily into pieces for mailing; these
had to be glued back together by the recipient.

Now the collection includes large differential equation and
optimization codes, packages written in C where specific filenames
matter, and large bibliographies. A variety of mail size limits are
appearing in the network. The use of Unix ``shell archives'' for
transferring collections of files has become widely understood and
accepted. In light of all this, it is clearly time for netlib to
adapt. This column announces three changes for dealing with the
problem of big files: ftp, revised mail splitting, and user-specified
mail limits.

We allow ftp access to a copy of the netlib files. Ftp is a file
transfer protocol widely available on machines connected to the
Internet. Using whatever syntax is appropriate for your local system,
connect to (also called and login using
the name ``netlib''. You'll be asked for your email address, so that we
can keep an accurate log of activity and send you bug fixes if, heavens
forfend, serious problems should be uncovered. Set ``binary'' mode, if
your version of ftp doesn't do that automatically. Then use the
``dir'', ``cd'', and ``get'' or ``mget'' commands to retrieve desired
files. Many have a ``.Z'' suffix indicating that they need to be
uncompressed after the transfer. If this all sounds too difficult or
confusing, just stick to email.

Netlib now uses a new mechanism for splitting up its reply into
multiple mail messages. The messages you will now receive are a shell
archives, scripts to be executed using the Unix ``sh'' command.
(People on other operating systems will have to edit manually.)
After executing each message, in arbitrary order, the netlib files will
be recreated on your system with their original names. In the case of
files too big to send in a single message, this is done by sending
temporary partial files; each script checks to see if all the partial
files have arrived and, if so, assembles them into the original big
file. But you don't need to worry about that; it is all automatic.

In some cases, an author's contribution to netlib is already in the
form of a shell archive. So whether you get it by ftp or by executing
mail messages, you still need to execute the file another time to do
the final unpacking. Yes, recursion can be confusing.

By default, netlib splits its reply into messages of at most 100
kilobytes. This satisfies most machines on the Internet. To change
the limit to, say, 20 kilobytes, you can give netlib the command
``mailsize 20k''.

Initially these changes are installed only at, so that
I can respond to suggestions for improvements. Assuming the new mail
splitting meets with general approval, it will be installed on the
other netlib nodes. As an alternative to ftp, Jack Dongarra will soon
be releasing an X-windows interface, which will allow downloading files

A final word regarding big files: if your mail connection to netlib still
depends somewhere along the path on dialup links, be sure you don't abuse
someone else's hospitality, lest your service be cut off.

\section*{Recent additions}

PCON61 (the latest version of Rheinboldt and Burkardt's well-known
program for continuation) has been installed. ``Program options
include the ability to search for solutions for which a given component
has a specified value. Another option is a search for a limit or
turning point with respect to a given component.''

In the {\bf opt} library, you will find Toint's {\bf ve08ad},
a subroutine which minimizes a partially separable objective function
with optional bound constraints.

Though the code is unchanged, the recent appearance of
Total Least Squares in SIAM's Frontiers book series merits notice.
For details, ``send index for {\bf vanhuffel}''.

{\bf Confpack}, by Gutknecht and Hough, implements conformal mapping
by numerical solution of Symm's integral equation. ``send index
for conformal.''

Gerard Holzman's {\bf spin} package for protocol validation is out of
the customary (scientific computing) scope of netlib but should
interest people in distributed computing.

The {\bf pvm} library has been installed. From the abstract:
``PVM stands for Parallel Virtual Machine. It is a software package
that allows the utilization of a heterogeneous network of parallel
and serial computers as a single computational resource. PVM consists
of two parts: a daemon process that any user can install on a machine,
and a user library that contains routines for initiating processes on
other machines, for communicating between processes, and synchronizing

TOMS algorithms 685 (separable elliptic PDEs) and 686 (updating
QR factorizations) arrived in February.

F2c, cheney-kincaid, kincaid-cheney, dassl, gcv, a/loess, fn/alngam,
fn/ci, misc/netlib, picl, pltmg, port/linpar, and numerous other files
have been updated.

In general, if you're wondering whether a file has been added in a
particular directory or whether it has been modified recently, you can
try a request like ``send directory from eispack''. That file is
automatically recreated nightly, and hence tells the truth even when
the netlib administrator is too busy to keep the ``index'' and
``changes'' files up to date.

{\em If you're new to netlib, send e-mail containing the line ``help''
to one of the Internet addresses
or uucp address {\tt uunet!research!netlib}. A few minutes later, assuming
you have speedy mail connections, you will receive information on how to use
netlib and an overview of the many mathematical software libraries and
databases in the collection.}

{\em Eric Grosse can be reached at the Computing Science Research
Center, AT\&T Bell Laboratories, Murray Hill NJ 07974, USA or by email
at {\tt}. Unix is a registered trademark of AT\&T
Bell Laboratories. This column was written August 24, 1991.}

}% end sloppy


End of NA Digest