Subject: NA Digest, V. 91, # 16

From: [email-masked] (NA-NET)
Date: Sun, 21 Apr 91 14:36:07 -0400

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

Today's Editor: Cleve Moler

Today's Topics:

     Last Call for SIAG/LA Prize
     Sparse Solvers for Overdetemrined Systems
     Block Tridiagonal Matrices
     Shampine/Gordon Integrator in C
     C Source for a Stiff ODE Integrator
     Three Measures of Precision in Floating Point Arithmetic
     Meeting on Mathematics, Computations, and Reactor Physics
     Second NIU Conference
     Position at the University of Auckland
     Postgraduate Opportunity at University of Durham
     Position at Centenary College of Louisiana
     Position at Lawrence Livermore Laboratory

-------------------------------------------------------

From: Mike Heath <[email-masked]>
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
Email:	[email-masked]

------------------------------

From: Rob MacLeod <[email-masked]>
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?

Thanks,

Rob.

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 <[email-masked]>
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
configurations.

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

e-mail: [email-masked]
        [email-masked]
        [email-masked]

------------------------------

From: Roland England <[email-masked]>
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
[email-masked]
[email-masked]

Tel. +44-908-65-2329

The Open University
Faculty of Mathematics
Milton Keynes
MK7 6AA
Great Britain

------------------------------

From: Steve Nichols <[email-masked]>
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
[email-masked] 

------------------------------

From: Nick Higham <[email-masked]>
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
predict.

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
found:

     >> 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-
                         .111
                        -----
                         .001

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

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 <[email-masked]>
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:

ADVANCES IN MATHEMATICS, COMPUTATIONS, AND REACTOR PHYSICS

		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 <[email-masked].
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 <[email-masked]>
Date: Wed, 17 Apr 91 10:13:28 NZS
Subject: Position at the University of Auckland

            THE UNIVERSITY OF AUCKLAND, NEW ZEALAND
                        LECTURESHIP
           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
department.

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.

[email-masked]        or      [email-masked]

------------------------------

From: Alan Craig <[email-masked]>
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
rigs.

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.

Dr.A.W.Craig
Department of Mathematical Sciences
University of Durham
South Road
Durham DH1 3LE
England

[email-masked]

------------------------------

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 <[email-masked]>
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: [email-masked].

------------------------------

End of NA Digest
**************************
-------