NA Digest Sunday, December 17, 1989 Volume 89 : Issue 49
Today's Editor: Cleve Moler
Today's Topics:
-------------------------------------------------------
From: Warren Ferguson <smu!ferguson@uunet.UU.NET>
Date: Mon, 11 Dec 89 12:19:12 CST
Subject: Ad Hoc Shift in EISPACK HQR
While playing around with the Eispack routine HQR I noted that some simple
eigenvalue problems cause nearly fatal convergence problems. I used HQR
to find the eigenvalues of the companion matrices:
Case 1:
0 0 0 -4 characteristic polynomial (x^2-2)^2
1 0 0 0 29 iterations to find first eigenvalue
0 1 0 4
0 0 1 0
Case 2:
0 0 0 -1 characteristic polynomial (x+1)^4
1 0 0 -4 20 iterations to find first eigenvalue
0 1 0 -6
0 0 1 -4
Case 3:
0 0 0 -16 characteristic polynomial (x+2)^4
1 0 0 -32 22 iterations to find first eigenvalue
0 1 0 -24
0 0 1 -8
HQR correctly found the eigenvalues in all cases. However, if it were not
for the detection of two consecutive small subdiagonal entries HQR would
not have been able to solve case 1. I considered these problems because I
thought they would give HQR problems and in this sense I was not
disappointed. My question is twofold. First, Wilkinson's HACLA article
on HQR does not seem to match the ALGOL code with regards to the handling
of exceptional shifts? HQR seems to explicitly shift the matrix and uses
some unexplained constants like 0.75 and 0.4375 when determining the
exceptional shift. Secondly, does anyone know what strategy will be used
in the LAPACK project? My problem is to understand what HQR is doing
when computing the exceptional shift and to find out what further progress
has been made in determing better exceptional shifts. Can anyone point me
in the right direction?
Warren Ferguson
ferguson%smu.edu
wferguson.nanet
------------------------------
From: Cleve Moler <moler@mathworks.com>
Date: Sun Dec 17 15:05:29 PST 1989
Subject: Reply to Ferguson on EISPACK HQR
In the above note, Warren Ferguson is correct that such matrices
stress EISPACK's routines HQR and HQR2. To generalize his
examples, let H(p,q,r) denote the Hessenberg companion matrix
of the polynomial (x^p - r)^q . The eigenvalues of H(p,q,r)
are the p-th roots of r, each repeated q times. Ferguson's
first example is H(2,2,2). Here is a Fortran code segment
that generates H(p,q,r):
integer p,q,i,j,n
real*.. r,s,H(..,..)
c
n = p*q
do 20 i = 1, n
do 10 j = 1, n
H(i,j) = 0.
10 continue
if (i .gt. 1) H(i,i-1) = 1.
20 continue
c
s = 1.
do 30 j = 1, q
s = -r*s*(q-j+1)/j
H(n+1-j*p,n) = -s
30 continue
When q > 1, the multiple roots cause linear convergence of the
basic Francis double shift QR algorithm used by HQR and HQR2. The
exceptional, or "ad hoc" shift, is intended to break up iterations
which don't converge at all. H(p,1,1) happens to be left unchanged
by the basic Francis iteration, and so needs the ad hoc shift
to get off dead center. But if r is not equal to 1, H(2,2,r)
will converge, though slowly, without the ad hoc shift.
For these matrices, it is even possible on some machines that
subtleties of roundoff error and the convergence test may cause
HQR to reach the iteration limit and report no convergence.
In this nonsymmetric situation, we do not have a complete analysis
of convergence that includes the effects of the ad hoc shift and
roundoff error.
I don't know what the LAPACK developers have planned, but I
suggest they stick with what is used in EISPACK unless they
develop a scheme which they can analyze completely and which
is *guaranteed* to be an improvement.
By the way, the ad hoc shift used in EISPACK's HQR is the same
as that used by Wilkinson in his Algol code in the Handbook.
I'm pretty sure Wilkinson's article describes the same thing,
but some algebra is required to go from one to the other.
The 0.75 factor in the code is one-half the 1.5 factor in the
article. Deriving the 0.4375 factor is left as an exercise.
(If you want to see a really obscure factor, look at the ad hoc
shift in the QZ algorithm. But that's another story.)
-- Cleve
------------------------------
From: Jens Tingleff <jensting@diku.dk)>
Date: 12 Dec 89 13:23:03 GMT
Subject: Zeros of Two Variable Polynomials
I am working on some image analysis, and have found it necessary
to use polynomials in two variables. The functions that I work
on are of the form
N - 1 N - 1
+---+ +---+
\ \ m n
f (x, y) = > > a x y
/ / mn
+---+ +---+
m = 0 n = 0
x, y and a Complex.
mn
I need to find the stationary points, and thus the zeros of a
two-polynomial system.
Any information on the properties of the zeros of polynomials
of this type would be appreciated. What I need are some simple
rules, and algorithms, such as those available for one-variable
polynomials.
I have found several articles on `finding all zeros of systems of
polynomial equations` (such as those by Li, Yorke, Carcia and others),
but these articles are vague on the kind of info that I want.
Info, pointers to literature, etc are needed.
Thank You in advance
Jens Tingleff
DIKU, U of Copenhagen, DK
------------------------------
From: T. R. Hopkins <trh%ukc.ac.uk@nsfnet-relay.ac.uk>
Date: Wed, 13 Dec 89 13:48:32 GMT
Subject: Netlib Service in the UK
Following recommendations made at the EASE Working Party on Numerical
and Mathematical Software at Warwick University in February 1988, the
University of Kent has been funded for two years to provide a pilot
numerical software distribution service. The service will formally start
on January 1st, 1990.
The software will be distributed over the UK's Janet via electronic mail
using the netlib system currently providing a similar services in the United
States. The software available is almost identical with that distributed
by Eric Grosse (thanks for the tape :-)) and Jack Dongarra.
A full list of the available libraries may be obtainable by sending a
mail message to netlib@uk.ac.ukc whose body contains the single line
send index
Any comments or suggestions may be sent to
netlib-suggest@uk.ac.ukc
Note: netlib-suggest is already available; netlib itself will not
officially come on-line until January 1st, 1990 although it is running as
a beta release at the moment.
Tim Hopkins and Maggie Bowman
Computing Laboratory
University of Kent
Canterbury, Kent, CT2 7NF
UK
------------------------------
From: Hans Schneider <hans@math.wisc.edu>
Date: Wed, 13 Dec 89 19:24:30 cst
Subject: International Linear Albebra Society
ILAS MEETING
The International Linear Algebra Society is pleased to announce its second
conference "LINEAR ALGEBRA, NUMERICAL LINEAR ALGEBRA AND APPLICATIONS" to
be held at Northern Illinois university, De Kalb, Illinois, U.S.A. The
tentative dates are April 25-28,1991.
The conference has been sponsored by Northern Illinois University.
Professors Biswa Nath Datta of Northern Illinois University and Hans
Schneider of University of Wisconsin are the co-chairmen of the organizing
committee.
The conference will cover all active research areas of core and numerical
linear algebra and their applications to systems and control theory,
signal processing, statistics, combinatorics, mathematical economics and
modeling, etc.
The purpose of the conference is to bring together researchers in linear
algebra, numerical linear algebra, and those working in various
applications areas for an effective exchange of ideas and discussions on
recent developments and future directions.
All (and particularly members of ILAS) are encouraged to send suggestions
for the meeting to either of the co-organizers.
------------------------------
End of NA Digest
**************************
-------