Subject: NA Digest, V. 91, # 17

From: mathworks!apple!surfer.EPM.ORNL.GOV!nacomb (NA-NET)
Date: Sun, 28 Apr 91 01:08:27 -0400

NA Digest   Saturday, April 27, 1991   Volume 91 : Issue 17

Today's Editor: Cleve Moler

Today's Topics:

     Netlib/NA-Net Machine Will Be Down a Few Days
     Stone's Strongly Implicit Procedure
     Solving a Nonlinear Initial BVP
     References Sought for Parallel SOR
     2nd Annual Midwest NA Day
     New Measures of Precision in Floating-Point Arithmetic

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

From: Jack Dongarra <[email-masked]>
Date: Sat, 27 Apr 91 12:20:40 -0400
Subject: Netlib/NA-Net Machine Will Be Down a Few Days

The Mathematical Sciences Section at Oak Ridge National Laboratory
is moving to a new building next week. As a result, the computer
running netlib and na-net at ORNL will be down starting Wednesday, 
May 1, 1991, for a few days.
We are automatically rerouting mail for [email-masked] to 
[email-masked], so netlib service will not be interrupted.
However, the na-net will be off the air until the move has been 
completed.  Mail sent to na.<lastname>@na-net.ornl.gov will be queued 
until the machine comes back up.  We hope to have the machine back 
by Friday, May 3.

   Jack

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

From: Glenn McKee <[email-masked]>
Date: Tue, 23 Apr 1991 10:30 EST
Subject: Stone's Strongly Implicit Procedure

Dear colleagues,

I am interested in speeding up a PDE algorithm that involves
solving a large, sparse matrix system.  I understand that 
a useful algorithm is a modified version of Stone's Strongly
Implicit Proceedure.

Our library lacks references to the this algorithm, and I was
wondering if someone had the algorithm already coded in either
C or FORTRAN and would consider sending me a copy.

Sincerely,

Glenn McKee
E-mail:  [email-masked]

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

From: R. Kaasschieter <[email-masked]>
Date: Tue, 23 Apr 91 16:44:39 +0200
Subject: Solving a Nonlinear Initial BVP

Dear colleagues,

It happens that we have to solve for our department of chemistry
the following initial boundary value problem on a three-dimensional
domain:

du/dt = div ( D(u) grad u )     in the domain,
- n . ( D(u) grad u ) = f(u)    on the boundary,
u = 1                           for t = 0.

D is uniformly bounded from below by a positive constant and
f(u) >= 0 for all u. 
A typical example for D is D(u) = exp(-u). Since the solution is known
to be positive, D(u) > 1.

We have planned to use finite elements for the space discretisation. 
Our main problem is a sensible choice for the time discretisation.

Do you have any experience in solving this problem or related problems,
or do you know relevant literature?

Rik Kaasschieter and Bob Mattheij

Please send your suggestions and comments to
  Dr. E.F. Kaasschieter
  Department of Mathematics and Computing Science
  Eindhoven University of Technology
  P.O. Box 513
  5600 MB  Eindhoven
  The Netherlands
  E-mail: [email-masked]

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

From: Renato De Leone <[email-masked]>
Date: Thu, 25 Apr 1991 07:33:03 CDT
Subject: References Sought for Parallel SOR

Hello na-netters,
I am looking for references on parallel Successive Overrelaxation
algorithms for SIMD machine. In particular for solving system of
linear equations with  specially structured (5-diagonals) matrices.
Any ideas or pointers to references would be greatly appreciated.

Thanks a lot
   Renato     
 
    Renato De Leone
    Computer Sciences Dept.
    University of Wisconsin-Madison
    1210 W. Dayton   Madison, WI 53705
    [email-masked]
    (608) 263-2677

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

From: Steven Lee <[email-masked]>
Date: Sat, 27 Apr 91 12:33:56 CDT
Subject: 2nd Annual Midwest NA Day

                        THIRD ANNOUNCEMENT
                             for the 
                      2nd Annual Midwest NA Day    

                      FINAL CALL FOR PRESENTERS
                   
Location: Champaign-Urbana, IL 
          Digital Computer Laboratory
          1304 W. Springfield Ave, Urbana		

Date:     Saturday, May 11 1991
          8:30AM - 4:30PM

The 2nd Annual Midwest NA Day is scheduled for Saturday, May 11 1991
on the University of Illinois, Urbana-Champaign campus. The first NA
Day was held in April of last year to take special note of the
retirement of Bill Gear from the University of Illinois.  This year,
the conference coincides with the May graduation ceremonies in which
Gene Golub is to receive an honorary doctorate from the University.

This is an informal conference with an emphasis on scheduling a
variety of talks from speakers with interests in numerical analysis
and scientific computing.

Those who wish to be included on the electronic mailing list for news
on this event (schedule information, travel directions, etc.) should
send email to:

[email-masked]

Again, for those who plan to attend and would also like to contribute
a 20-minute or 40-minute presentation, please send the title and a
brief abstract of your proposed talk to the e-mail address listed
above.  We have already had numerous responses to the previous
announcement in NA Digest; we hope to finalize the program and list of
speakers within the next few days.

Organizers: Paul Saylor, Robert Skeel, Faisal Saied, Ahmed Sameh,
Michael J. Holst and Steven Lee

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

From: Bernard Danloy <[email-masked]>
Date: Tue, 23 Apr 91 and Fri, 26 Apr 91
Subject: New Measures of Precision in Floating-Point Arithmetic

In the last NA-Digest (21 April ; Vol.91/16), two communications were devoted to
the precision of the floating point arithmetic of a computer.
I would like to thank Nick Higham and Cleve Moler for their attempt to define
unambiguous concepts ; the lack of well defined concepts is too often the source
of confusion and errors. In addition, I am very happy to have learned that
the number mu defined by Higham and Moler is not always a negative power of the
base b : I never payed enough attention to that quantity by itself and was
ready to claim the opposite ; I just hope I was not the only one.

However it seems to me that Higham and Moler did not go as far as needed : they
refer to the basic concept of "floating point number" : this is meaningful in a
classical environment where the arithmetical unit and the memory use exactly
the same representation for real numbers. If this is the case, it is irrelevant
to store the intermediate results of a process in a variable or to keep them
within the arithmetic unit ; the quantities eps,mu and u are unambiguous if we
define the floating point numbers as the possibles values a variable of real
type can take.

But it should be reminded that not all computers work that way : it may happen
that the arithmetical unit can handle quantities it is impossible to save in a
variable without loss of information : such a situation was mentioned by Higham
and Moler and arises when the memory uses the basic IEEE standard ( 53 bits )
while the arithmetical unit conforms to the extended IEEE standard ( 64 bits ).
How should we define eps,mu and u in such a context ?

We are then facing two possible values for eps and u ; the values derived from
the characteristics of the memory appear as a reasonable choice but the values
related to the arithmetical unit are not a useless information.

The case of mu is much more difficult because the number of possible values is
almost unpredictable : mu is the result of some computation and therefore its
value depends on the whole sequence of operations : it is then crucial to know
exactly which computation has to be performed and how it is implemented. On one
single machine, the number of computations resulting in different values for mu
can be far greater than 2 ! I show now what you obtain on a Sun 3/60 computer
equipped with an optional arithmetical coprocessor Motorola 68881 when you are
looking for mu through Pascal programs. Here follows a summary of my results
( x and y denote real variables )

Without using the arithmetical coprocessor :

1. The boolean expression 1+x>1 is true iff  x >= 2^(-53) + 2^(-105)
  
2. When y is precomputed by y := 1+x
   the boolean expression  y>1  is true iff  x >= 2^(-53) + 2^(-105)

Using the optionnal arithmetical coprocessor :

3. The boolean expression 1+x>1 is true iff  x >= 2^(-64) + 2^(-116)

4. When y is precomputed by y := 1+x
   the boolean expression  y>1  is true iff  x >= 2^(-53) + 2^(-64) + 2^(-105)

5. The boolean expression 1+(x+y)>1 is true in many situations, for example,
     a.  if x > 2^(-64) and y >= 0
     b.  if x = 2^(-64) and y >= 2^(-127)
     c.  if x = 2^(-64) - 2^(-117) and y >= 2^(-117) + 2^(-128) + 2^(-169)

The above results clearly show that defining mu as the smallest quantity such
that 1+mu>1 is ambiguous if we don't say precisely :
- how mu may be obtained and used
- where the result of the addition has to be stored before it is compared to 1

I am far from being an expert but I suggest to adopt the following definition :
  if x is a variable of real type, mu is the smallest value that x can take if
  we want the computer to give the value true to the boolean expression  1+x>1

That corresponds to the situations 1 and 3 above ; the fact that mu does not
reflect the representation of the floating point numbers within the memory is
not disturbing at all : comparing 1 and 1+x is part of a tricky way to get
information about that representation, assuming a unique standard is valid
everywhere. Maybe, it is time to realize that, when this is not the case, the
errors generated by a single computation do have two independent sources :
	- the arithmetical operation itself
	- the transfer of numbers to ( and from ? ) the memory
It is then perfectly normal to use two independent parameters mu and eps to
characterize the precision of a floating point arithmetic : the value of mu
says how precisely a result can be obtained within the arithmetical processor
and eps says how precisely that result can be stored in the memory.

  ************

Date: Fri, 26 Apr 91 20:15:11 +0200

I would like to add new comments and proposals to my first note on the measure
of precision of a floating-point arithmetic.

I just realized that if we want to characterize a t-digit base b arithmetic, we
need four and not only three parameters. Two of them are concerned with "static"
properties ( Cleve Moler did use very properly the term " geometric " ) and the
two others are concerned with " dynamic " properties. The problem we are facing
now is just that, as a result of historical evolution, a confusion has been
created between " static " and " dynamic " ; the best example is perhaps the
quantity u ( unit roundoff ) : it is clearly a " dynamic " concept but it has
been defined until now by a  " static " value : I suggest we try to get rid of
our habits and focus on what we really need. Here are my proposals :

" STATIC " ( or " GEOMETRIC " ) ASPECT

b, t and eps = b^(1-t) are clearly the basic parameters ; two of them define
the third. I suggest to characterize the " static " properties of the floating-
point arithmetic by
   (1) eps = b^(1-t) : this defines ( within the memory ) the spacing of the
	floating-point numbers between 1 and b
	Important remark : eps has nothing to do with the equivalent spacing
	   within the arithmetical unit.
   (2) b : the base ( the knowledge of b is only needed if we want to get error
	bounds as precise as possible : for a same value of eps, the precision
	attached to some computations is a bit higher if the base b is smaller )

" DYNAMIC " ASPECT

Any arithmetic computation usually requires three things :
	- a transfer of the data from the memory to the processor
	- the execution of the operation itself ( within the processor )
	- a transfer of the result from the processor to the memory
With a " dynamic " point of view, the precision may be characterized by two
independant parameters :
    (3) u is the unit roundoff of the memory and gives information about the
	rounding which occurs when transfering a number from the processor to
	the memory ( let us remember that reading data on a peripheral usually
	implies a conversion to the base b and involves thus the processor ).
	I suggest to define u as the smallest floating-point number x such that
	true is the computed value of the boolean expression fl(y>1) if y has
	been previously obtained by fl(1+x), stored in the memory and picked
	back up ( this double transfer is crucial )

	N.B. The old definition of u as 0.5*b^(1-t) is a very special case of
	the new one and is equivalent in the case of rounding away from zero.
	( What was considered as perfect twenty years ago, before rounding to
	even appeared ).

	Actually the new value is equal to u = rho * eps where rho is
	    1  in the case of pure chopping
	   1/2 in the case of rounding away from zero
	   something slightly bigger than 1/2 in the case of rounding to even
	The new definition of u makes it equal to the quantity mu introduced by
	Nick Higham provided that the computations are performed in two steps
	as required. ( By the way, Nick did the two steps since he asked the
	computer to perform the addition and made himself the comparison ).

    (4) nu is the unit roundoff of the arithmetical processor when it performs
	a computation on data extracted from the memory. It is defined as the
	smallest floating-point number x such that true is the computed value of
	the boolean expression fl(1+x>1) where the unusual notation ( I wrote
	fl(1+x>1) and not fl(1+x)>1 ) means that the user does not force the
	intermediate result fl(1+x) to be sent back to the memory ( the computer
	is given complete freedom and gets the due credit for its choice ).

	Actually nu is equal to u if the computer ( I should probably have to
	say the software ) always transfers intermediate results back to the
	memory and never tries to avoid it. The value of nu is equivalent to
	that of mu if the computation is performed in one single step : we just
	ask the computer to evaluate the complete expression as a whole.

As I did already mention it in my first note, I ran some computations on a Sun3
with a Motorola chip 68881. I found

Using Pascal with the otion -f68881 activated :
	       eps = 2^(-52)
	       b = 2
	       u = 2^(-53) + 2^(-64) + 2^(-105)
	       nu = 2^(-64) + 2^(-116)

Using Pascal without the otion -f68881 :
	       eps = 2^(-52)
	       b = 2
	       u = 2^(-53) + 2^(-105)
	       nu = 2^(-53) + 2^(-105)

Using MATLAB ( and the Motorola chip ) : 
	       eps = 2^(-52)
	       b = 2
	       u = 2^(-53) + 2^(-64) + 2^(-105)
	       nu = 2^(-53) + 2^(-64) + 2^(-105)

Partial tests have been run on a PC and confirm the last values for MATLAB.
Further trials will performed later.

It turns out that the " static " parameters are independant of the software
and characterize a standard like IEEE, as they should ; on the opposite, the
" dynamic " parameters do vary very strongly : even on a single machine, they
may depend on the software but also on the options activated by the user.

At this point, one thing is clear for me : we really need to agree on precise
concepts and the topics are worth a much deeper discussion. I hope NA-NET will
contribute to it.

Bernard Danloy
Institut de Mathematique
University of Louvain-la-Neuve
Belgium
Email : [email-masked]

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

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