NA Digest Sunday, December 11, 1994 Volume 94 : Issue 50

Today's Editor:

Cleve Moler

The MathWorks, Inc.

moler@mathworks.com

Today's Topics:
• Greengard Method
• Random Arithmetic and x*(1/x)
• The Bad Scenario for FDIV--Simplified Account
• Sparse Matrix Times a Packed Vector
• Despertly Seeking Matrices
• Reminder of Deadline for Benelux Meeting
• Conference on Computer Methods in Water Resources
• Scientific Computation and Differential Equations
• Symposium on Systems Analysis and Simulation
• Short Course in Numerical Computing with MATLAB
• Interface 95 Announcement
• Bristol-Bath Numerical Analysis Day
• Positions at Bath University
• Positions at Florida State University
• Position at Cornell Theory Center
• Position at Colorado State University
• Position at University of Delaware
• Position at University of Tennessee

Submissions for NA Digest:

Information about NA-NET:

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

From: Leslie Greengard <greengar@greengard.cims.nyu.edu>
Date: Thu, 8 Dec 94 11:43:31 -0500
Subject: Greengard Method

A message on the NA-NET last week referred to "the Greengard method"
and asked for a reference. The algorithm being referred to is
most likely the fast multipole method (FMM). It was initiated
in work by V. Rokhlin

V. Rokhlin (1985), Rapid solution of integral equations of classical
potential theory, J. Comput. Phys. 60, pp. 187-207,

and generalized in a series of subsequent papers

L. Greengard and V. Rokhlin (1987),
A fast algorithm for particle simulations,
J. Comput. Phys. 73, pp. 325-348.

J. Carrier, L. Greengard, and V. Rokhlin (1988),
A fast adaptive multipole algorithm for particle simulations,
SIAM J. Sci. Statist. Comput. 9, pp. 669-686.

L. Greengard and V. Rokhlin (1988),
Rapid Evaluation of Potential Fields in Three Dimensions,
in Vortex Methods, C. Anderson and C. Greengard (eds.),
Lecture Notes in Mathematics, vol. 1360,
Springer-Verlag, pp. 121-141.

The work described in these papers is also available in book form:

L. Greengard (1988),
The Rapid Evaluation of Potential Fields in Particle Systems,
MIT Press, Cambridge.

Because many people seem to have learned of the FMM from the last source,
it is occasionally referred to as "the Greengard method." It should
not be. On the book jacket, the name "Rokhlin-Greengard" algorithm
is used, but "fast multipole method" is both more descriptive and
more common usage.

Leslie Greengard


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

From: Alan Edelman <edelman@math.mit.edu>
Date: Tue, 6 Dec 94 13:42:11 EST
Subject: Random Arithmetic and x*(1/x)

In discussing the Pentium FDIV bug Cleve Moler wrote:

>> There are serious and subtle numerical difficulties associated
>> with using even a correct multiply and subtract to decide if a
>> proposed quotient is correct. For example, simply replacing
>> N/D by N*(1/D) won't work if you after the correctly rounded,
>> IEEE floating point result.

For a fun brief education on such subtle issues, I whipped up a little
note on when x*(1/x)=1 which may be clicked on from

http://theory.lcs.mit.edu/~edelman/comprehensive.html

On my own personal wishlist is that every theory of computation
student would hear one lecture on the IEEE standard for floating
point computation. Better yet he or she should be obliged to
take a numerical analysis course.

And, on the same subject, Manny Blum wrote

>> We must demand radically stronger standards of reliability.
>> But how is such reliability to be achieved? We believe that the
>> answer is via run-time result-checking and self-correcting. In
>> particular, we believe that the chip of the future will have efficient
>> checkers and correctors embedded in its hardware. For example, every
>> time a chip divides, it will do a multiplication to check the answer,
>> and, if the answer is incorrect, will call a randomized self-corrector
>> such as the one described above. Note that testing the chip's
>> division before release will then be a simple matter: all the company
>> need do is run the chip on a number of random inputs and verify the
>> correctness of each answer. It is thus assured that the chip can only
>> be incorrect on at most a small portion of all possible inputs; and a
>> randomized self-corrector will be capable of entirely suppressing such
>> occasional bad cases.

To follow Einstein, I don't like playing dice with arithmetic. (Yes,
I know Einstein was wrong on the "dice" issue.) My fear would be that
we would be lulled into producing sloppy arithmetic because we expect
the randomizer to fix it for us. Given that correct arithmetic is
achievable by careful human beings, I am not in favor of bogging down
chips with randomized error fixers. Some things are worth just getting
right.


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

From: Vaughan Pratt <pratt@Sunburn.Stanford.EDU>
Date: 4 Dec 1994 07:56:12 GMT
Subject: The Bad Scenario for FDIV--Simplified Account

In article <D09Bv0.23G@utnetw.utoledo.edu>, Tim Arheit
<tarheit@cse.utoledo.edu> wrote:

>The problem is that my software doesn't use random number as data. One
>of the postes in this group shows that simple division of ~integer numbers
>shows a significant error (something on the order of 100 times less accurate)
>4.99999999/14.999999 (5/15) (I think these are some of the example
>numbers given, it's in this group somewhere) generated significant error.
>The point is that integer constants in a floating point equation are far
>from random.


Let me correct a minor detail and at the same time repackage what I
said in two earlier messages in a hopefully more digestible and
memorizable form.

The exact example here subtracted one millionth (easy to remember) from
each of 5 and 15, i.e. 4.999999/14.999999. The correct value is
0.33333329, Pentium math makes it 0.333329, wrong in the fifth decimal
place, specifically a relative error of 1.2 in 10^5.

Other such rationals are 7/48 and 9/54. For each of these, if you
subtract one millionth from both numerator and denominator before
dividing, you again obtain a relative error of about one in a hundred
thousand.

Among rationals with numerator and denominator both bounded by 1000,
each decremented by a millionth, there are 14 that behave as badly as
the three examples above, and an additional 412 that the Pentium
evaluates with a relative error between 10^-5 and 10^-7. (For
comparison, using single precision instead of double yields a maximum
relative error of 6*10^-8, so all these errors are worse than would
result merely from using single precision instead of double.)

Hence if you compute with randomly chosen approximate (in the above
sense) rationals from that population, the probability of an error in
the fifth place is one in 70,000, while the probability of an error in
the sixth or seventh place is one in 2,000.

The P90 can do better than two divisions per microsecond, at which
rate, assuming the above population of approximate rationals, the small
error (sixth or seventh place) occurs once a millisecond on average
while the large error (fifth place) occurs once every 35 milliseconds.

When the division operands are bounded by 100 instead of 1000, the
1/2,000 figure for the small error drops to 1/200 (a factor of 10
worse) while the 1/70,000 rate for the large error drops to 1/2,000 (a
factor of 35 worse), i.e. the large errors now happen once a
millisecond.

Detailed tables of such errors, along with a program for verifying my
conclusions for yourself, can be found in my two earlier messages on
this subject on comp.arch,comp.sys.intel. Search for the phrase
"ibruis" (a variable in my program) or the subject "A natural scenario
with high FDIV bug probability (was: In Intel's Defense...)".

Vaughan Pratt http://boole.stanford.edu/boole.html


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

From: Hassane Sadok <Hassane.Sadok@univ-lille1.fr>
Date: Tue, 6 Dec 1994 14:08:38 --100
Subject: Sparse Matrix Times a Packed Vector

Dear Colleagues,
I would like to know what is the best way for implementing in Fortran
a product of a sparse matrix stored in compressed sparse row (CSR) format by
a packed vector.
Is the best way, the one achieved by first expanding the packed vector into a full-lenght vector?
Many thanks for your help.

Hassane SADOK
sadok@ano.univ-lille1.fr
Laboratoire d Analyse Numerique et d Optimisation
Universite des Sciences et Technologies de Lille.
Bat M3. 59655 Villeneuve d Ascq cedex.
France.


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

From: George Miminis <george@triton.cs.mun.ca>
Date: Tue, 6 Dec 1994 19:42:51 -0330
Subject: Despertly Seeking Matrices

Has anybody encountered the following class of (2^k)x(2^k) matrices, (for
some natural number k), defined recursively as follows:

(i): If k=1 then

(a -b)
(b a)

is such a matrix.

(ii): If A, B are such matrices of size 2^k each, then

(A -B)
(B A)

is such a matrix of size 2^(k+1).

It is well known that complex numbers are defined as 2x2 such matrices. As
an example of (2^3)x(2^3) such matrix consider

                         ( 1|-2|-3  4|-5  6  7 -8)

                         (--|--|     |           )

                         ( 2| 1|-4 -3|-6 -5  8  7)

                         (--|--|-----|           )

                         ( 3|-4| 1 -2|-7  8 -5  6)

                         (--|--|     |           )

                         ( 4| 3| 2  1|-8 -7 -6 -5)

                         (--|--|-----|-----------)

                         ( 5|-6|-7  8| 1 -2 -3  4)

                         (--|--|     |           )

                         ( 6| 5|-8 -7| 2  1 -4 -3)

                         (--|--|-----|           )

                         ( 7|-8| 5 -6| 3 -4  1 -2)

                         (--|--|     |           )

                         ( 8| 7| 6  5| 4  3  2  1)