Parallel computation of pseudospectra by fast descent.
The pseudospectrum descent method (PsDM) is proposed, a new parallel method for the computation of pseudospectra. The idea behind the method is to use points from an already existing pseudospectrum level curve  to generate in parallel the points of a new level curve  such that  > . This process can be continued for several steps to approximate several pseudospectrum level curves lying inside the original curve. It is showed via theoretical analysis and experimental evidence that PsDM is embarrassingly parallel, like GRID, and that it adjusts to the geometric characteristics of the pseudospectrum; in particular it captures disconnected components. Results obtained on a parallel system using MPI validate the theoretical analysis and demonstrate interesting load-balancing issues.
a region of the complex plane and ii) compute  min (zI A) for every node
zh . The boundary curves @  (A) are obtained as the contour plots of the
smallest singular values on the
mesh
h . The obvious advantage of GRID is its
straightforward simplicity and robustness. If we let for the moment C min be
This work has been partially supported by the Greek General Secretariat for Research
and Development, Project ENE 99-07
a measure of the average cost for the computation of  min (zI A), it is easily
seen that the total cost of GRID can be approximated by
where
denotes the number of nodes
of
h . The total cost quickly becomes
prohibitive with the increase of either the number of nodes or the size of A.
It has been observed that the cost formula (2) readily indicates two major
methods for accelerating the computation: a) reducing the number of nodes z
and hence the number of evaluations of  min , and b) reducing the cost of each
evaluation of  min (zI A). Both approaches are the subject of active research;
see [15] for a comprehensive survey of recent eorts.
The use of path following, a powerful tool in many areas of applied mathe-
matics, in order to compute a single boundary curve @  (A) was suggested by
Kostin in [8]. It was M. Bruhl in [5] who presented an algorithm to that end
and showed that signicant savings can be achieved compared to GRID when
seeking a small number of boundary curves. The key is that tracing a single
boundary curve, drastically reduces the number of  min evaluations. Bekas and
Gallopoulos carried this work further in Cobra [2], a method that ameliorated
two weaknesses of the original path following method, in particular i) its lack of
large-grain parallelism and ii) its frequent failure near sharp turns or neighboring
curves. In particular, parallelism was introduced by incorporating multiple,
corrections. More recently, D. Mezher and B. Philippe suggested
PAT, a new path following method that reliably traces contours [13]. Despite
its advantages over the original method of [5], Cobra maintains two
of the path following approach vis-a-vis GRID. These are a) that in each run, a
single boundary curve @  (A) is computed, b) that is closed. Therefore, only
one curve is computed at a time and its disconnected components are not cap-
tured, at least in a single run. While we can consider the concurrent application
of multiple path following procedures to remedy these problems, the solution is
less straightforward than it sounds, especially for (b). In the remainder of this
paper, when our results do not depend on the exact version of path following
that we choose to use, we denote these methods collectively by PF.
In this paper we propose the Pseudospectrum Descent Method (PsDM) that
takes an approach akin to PF but results in a set of points dening pseu-
dospectrum boundaries for several values of . PsDM starts from an initial
contour @  (A), approximated by N points z k computed by a PF method.
These points are corrected towards directions of steepest descent, to N points
. The process can be repeated recursively, without the need
to reuse PF, and computes pseudospectra contours. PsDM can be viewed as a
\dynamic" version of GRID where the stride from mesh point to mesh point has
been replaced by a PF Prediction-Correction step. To illustrate this fact and
thus provide the reader with an immediate feeling of the type of information
obtained using the method, we show in Figure 1 the results from the application
of PsDM to matrix kahan of order 100 which has been obtained from the
Test Matrix Toolbox ([7]) and is a typical example of matrices with interesting
Y

Figure

1: Pseudospectrum contours and trajectories of points computed by PsDM
for @  (A); kahan of order 100. Arrows show the
directions used in preparing the outermost curve with path following and the
directions used in marching from the outer to the inner curves with PsDM. See
section 4 for further results with this matrix.
pseudospectra. The plot shows i) the trajectories of the points undergoing the
steepest descent and ii) the corresponding level curves. The intersections are
the actual points computed by PsDM.
We note that the idea for plotting the pseudospectrum using descent owes
to an original idea of I. Koutis for the parallel computation of eigenvalues using
descent described in [3].
The rest of this paper is organized as follows. Section 2 brie
y reviews path
following. Section 3 describes PsDM. Section 4 illustrates the characteristics
of PsDM as well as its parallel performance using an MPI implementation. We
also introduce adaptivity and show that the method can reveal disconnected
components. Section 5 presents our conclusions.
2 Review of Path Following
Consider the function G
According to denition of the pseudospectrum, the zeros of G(z) are points of
the boundary @  (A). Allgower and Georg in [1] describe a generic procedure

Figure

2: A generic Prediction-Correction PF scheme.
to numerically trace solution curves of equations such as (3). In Table 1 we
outline the algorithm and illustrate one step in Figure 2. See [2, 5] for more
details.
z 0 on @  (A).
for k=1,.
(* Prediction phase *)
1.1 Determine a prediction direction p k .
1.2 Choose a steplength h and predict the point ~
(* Correction phase *)
2.1 Determine a correction direction c k .
2.2 Correct ~ z k to z k using Newton iteration.

Table

1: PF generic procedure for the computation of a single contour.
In the sequel, as in [5], we would be identifying the complex plane C with
R 2 and frequently use the notation G(z) for G(x; y). Critical to the eective use
of PF is the availability of the gradient rG(x; y); that this becomes available at
little cost follows from the following important result (cf. [5, 6]):
Theorem 2.1 Let z 2 C n(A). Then G(x,y) is real analytic in a neighborhood
of is a simple singular value. Then the gradient of
G at z is
min u min ); =(v
min u min )): (4)
where denote the left and right singular vectors corresponding to
min (zI A).
Therefore, assuming that the minimum singular triplet
A (in the sequel we would be referring to it simply as \the triplet at z") has
been computed, the gradient is available with one inner product.
From the above it follows that the dominant eort in one step of the PF
scheme proposed in [5] amounts to i) the estimation of the prediction direction
and ii) the Newton iteration of the correction procedure. Regarding (ii) we
note that it is su-cient (cf. [5, 2]) to use a single Newton step, requiring one
triplet evaluation at ~
z k . Then, as proposed in [5], ~
z k is corrected towards the
direction of steepest descent c k , i.e.
z k
min
where the triplet ( min ; u min ; v min ) is associated with ~
z k . Regarding (i), note
that selecting p k to be tangential to the curve at z k 1 requires the triplet at
z k 1 . Since this would double the cost, it has been shown acceptable to take p k
orthogonal to the previous correction direction c k 1 ; cf. [5, 2]. Therefore the
overall cost of a single step of the original method of [5] is approximately equal
to the cost for computing the triplet.
3 The Pseudospectrum Descent Method
Let us now assume that an initial contour @  (A) is available in the form of
a some approximation (e.g. piecewise linear) based on N points z k previously
computed using some version of PF. In order to obtain a new set of points that
dene an inner level curve we proceed in two steps:
1: Start from z k and compute an intermediate point ~
w k by a single modied
Newton step towards a steepest descent direction d k obtained earlier.
Step 2: Correct ~
w k to w k using a Newton step along the direction l k of steepest
descent at ~

Figure

3 illustrates the basic idea for a single initial point. Applying one Newton
step at z k would require r min ((x k therefore, a triplet evaluation
at z k . To avoid this extra cost, we apply the same idea as PF and use instead
min q min ), that is already available from
the original path following procedure. In essence we approximated the gradient
based at z k with the gradient based at ~ z k . Vectors g min ; q min are the right and
left singular vectors associated with  min (~z k I A). Applying correction as in
(5) it follows that
~
min g min

Figure

3: Computing @ - (A); - < .
is the new pseudospectrum boundary we are seeking.
Once we have computed ~
we perform a second Newton step that yields w
min ( ~
where the triplet used is associated with ~
w k . These steps can be applied to all N
points in what we call one sweep of PsDM; we denote it by PsDM and outline it
in the following Table. Starting from an initial contour @  (A) we have shown
of the initial contour @  (A).
on the target contour @ - (A); - < .
for k=1,. ,N
1. Compute the intermediate point ~
according to (6).
2. Compute the target point w k using (7).

Table

2: PsDM: One sweep of PsDM.
how to compute points that approximate a nearby contour @ - (A); - < .
Assume now that the new points computed with one sweep of PsDM dene
satisfactory approximations of @ - (A) (cf. end of section). We ask whether it
would it be practical to use these points to march one further step to approximate
another
As noted in the previous discussion, the
application of the sweep PsDM uses r min ((~x k i.e. the triplet at
~
z k . This is readily available when the curve @  (A) is obtained via PF. Observe
now that as the sweep proceeds to compute @ ~ - (A) from @ - (A), it also computes
the triplet at ~
w k . Therefore enough derivative information is available
for the sweep PsDM to proceed one more step with starting points computed
Method PsDM
approximating a contour @  (A).
approximating M contours
for i=1,. ,M
Compute points of @ - i by PsDM on the N points of @ - i 1

Table

3: The PsDM method.
via the previous application of PsDM. Therefore, it is not necessary to run PF
again. Continuing with this repeated application of sweeps PsDM, we obtain the
promised pseudospectrum descent method, outlined in Table 3 and illustrated
in

Figure

4.
It is worth noting that each sweep can be viewed as a map that takes as
input N in points approximating  - i (A) and produces N out points approximating
In the description so far N in = N out , but as we show in
Section 4.4, this is not necessarily an optimal strategy.
We nally note that the two step process described above was found to be
necessary. By contrast, the more straightforward approach in which compute the
exact steepest descent directions from each point at the starting curve and use
those values directly to approximate the next set of points in one step required
a much smaller stepsize to be successful. We did not experiment further with
this procedure.
Cost considerations It is evident that the cost for the computations of the
intermediate points ~
relatively small since the derivatives
at ~ z have already been computed by PF or a previous sweep of
PsDM. Furthermore, we have assumed that  min (z k I
the points z k approximate @  (A). On the other hand, computing the nal
points evaluations. To avoid proliferation
of symbols, we use C min to now denote the average cost for computing the
triplet. Therefore, we approximate the cost of a single sweep of PsDM by C
NC min . The computation of each target point
from the computation of all other target points w k. On a system with P
processors, we can assign the computation of at most dN=P e target points to
each processor; one sweep will then proceed with no need for synchronization
and communication and its total cost is approximated by C
We discuss this issue in more detail in Section 4.3. We emphasize that typically,
the number of points N on each curve is expected to be large, therefore the
algorithm is scalable. This is better than the PF methods described in the
introduction: Cobra allows only a moderate number of independent calculations
of  min per step while in the original version of PF, coarse-grain parallelism

Figure

4: Descent process for a single point
becomes available only if we attempt to compute dierent curves in parallel.
analysis for a single sweep of PsDM In order to gauge the quality of
the approximation of the curve obtained in a single sweep of PsDM, it is natural to
use as a measure the value which we call \stepsize of the sweep". Let
be the point approximated starting from a point z 2 @ - i (A).
Dene the function G 0
be the point obtained using a single Newton step at ~
as depicted in Figure
real analytic in a domain D if the minimum singular value
min (wI A) is simple for all w 2 D. Let all ve points ~ z; z; ~
Fig. 5) lie in the interior of such a domain. From the relations
z
and ~
z
it follows that w 0 can also be considered as the outcome of two exact Newton
steps originating from ~ z. From the identity (w 0 I
and standard singular value inequalities, it follows that
min (w 0 I
and therefore From results of
Levin and Ben-Israel (see [12]) regarding Newton's method for underdetermined
systems, under standard assumptions there will be a region in which the Newton
iteration used to produce w 0 from ~
w will converge quadratically; we next assume
that we are within this region. Then jw 0 wj   1 jw 0 ~
w y )k); note that 1=krG 0 ( ~
w y )k is at least 1, since G 0 ( ~
This gure can be considered as an enlargement of Figure 3 that reveals that the point
computed with a single Newton step in the sweep produces only an approximation of w k .

Figure

5: The transformation of z to w 0 via PsDM.
is the inner product of singular vectors (cf. Theorem 2.1). Therefore, using
relations (7) and (9) it follows that
wI A). It follows that
if this minimum is not much smaller than 1 and the assumptions made above
hold, then the error induced by one sweep of PsDM will be bounded
by a moderate multiple of the square of the stepsize of the sweep. A full scale
analysis of the global error, for multiple sweeps of PsDM, is the subject of current
work.
4 Renements and numerical experiments
We conducted experiments with matrices that have been used in the literature
to benchmark pseudospectra algorithms. Our results, presented below, indicate
that PsDM returns the same levels of accuracy as GRID at a fraction of the cost.
We then describe a parallel implementation of PsDM. The speedups obtained underscore
the parallel nature of the algorithm. We also show that the application
of PsDM to large matrices for which the triplets have to be computed via some
iterative method is likely to suer from load imbalance and speedups that are
than what one would expect from an embarassingly parallel algorithm;
we describe a simple heuristic to address this problem.
Finally we discuss some other properties that underline the
exibility of
PsDM. We show in particular that: i) adaptation of the number of points computed
in each sweep of PsDM can lead to signicant cost savings, and ii) PsDM
can capture disconnected components of the pseudospectrum lying inside the
initial boundary computed via PF.
4.1 System conguration
We performed our experiments on a SGI Origin 2000 system with 8 MIPS
R10000 processors. The system had a total of 768 MB RAM and 1MB cache
per processor, running IRIX 6.5. The codes were written in Fortran-90 using
F90/77 MIPSpro version 7.2.1 compilers. For the parallelization we used the
MPI programming paradigm, implemented by SGI's MPT 1.2.1.0. We used
ARPACK [11] to approximate the triplets and SPARSKIT [14], suitably modied to
handle double precision complex arithmetic, for sparse matrix{vector multiplies.
All our experiments were conducted in single-user mode.
4.2 Numerical experiments with PsDM
We remind the reader that for matrices with real elements, the pseudospectrum
curves are symmetric with respect to the real axis, and it has become standard
in the pseudospectrum literature to measure the success of methods for pseudospectra
by direct comparison of the one half of the gure computed with a
new method and the other half with GRID. Therefore, in subsequent experiments
with any method we only compute points and curves lying on the upper or lower
half of the complex plane. We start with the same matrix used to obtain Figure
1, that is kahan of order 100. We used the parallel version of PF, namely Cobra,
to obtain points to approximate (one half of) the pseudospectrum curve
corresponding to asked PsDM to compute 60 contours corresponding
to values of  ranging from - using a
stepsize that remained equal
6. The upper half of gure 6 illustrates the contours
corresponding to while the the lower half illustrates the
same contours computed using GRID on a 100  100 mesh of equidistant points.
We chose this resolution in order to allow GRID to oer a level of detail that is
between the minimum and maximum resolution oered by PsDM. In particular,
the distance between neighboring points of GRID was selected near the median of
the smallest and largest distance of neighboring points computed in the course
of PsDM. The pictures are virtually indistinguishable and indicate that PsDM
achieves an accuracy comparable to GRID. On the other hand, the cost is much
smaller; in particular PsDM approximates the contours using 2700 points while
GRID used 10000 mesh points. It is also worth noting that despite the fact that
GRID does not require the computation of singular vectors, in the context of
iterative methods such as ARPACK the extra cost is not signicant and PsDM is
expected to be far less expensive.
In order to further analyze the accuracy of PsDM, we computed the relative
error jmin
at each computed point z of the (approximate) curves
Y
PsDM
GRID

Figure

Selected pseudospectrum contours @  (A);
kahan of order 100 computed by PsDM (top) and GRID (bottom).
produced by the algorithm; from these values, we obtained the maximum
and mean relative errors for each contour point and show the results in

Figure

7. Note that the best way to read this gure is from left to right, as
this shows how the maximum and mean errors per curve develop as we consider
the curves from the outermost to the innermost. The maximum relative
error appears to be satisfactory except in a few cases, where the maximum error
even then, however, the mean error remains two orders
of magnitude smaller. This highlights one observation we made in our experi-
ments, namely that only a very limited number of points suer from increased
relative error. We also note that there is an apparent increase in the error as we
approach very small values of . This does not re
ect a weakness of PsDM, but
the di-culty of the underlying SVD method, ARPACK in this case, to compute
approximations of the minimum singular value with very small relative error.
4.3 Parallel performance
One important advantage of PsDM is that, by construction, it is embarrassingly
parallel: Each sweep can be split into a number of tasks equal to the number of
points it handles and each task can proceed independently with its work, most
of it being triplet computations. Furthermore, assuming that no adaptation is
used, no communication is needed between sweeps. Let us assume that this is
the case and that the number of points computed in each sweep and therefore the
number of tasks is constant, say N . It is thus natural to use static partitioning
e
Relative
error
Maximum relative error
Mean relative error

Figure

7: Maximum and mean relative errors for each curve (60 total) computed
via PsDM for
(*Input*) N points z k on the starting number of processors
1. For each processor, select and assign the sweep computations for
dN=P e points.
2. Each processor proceeds with PsDM on its own set of points.
3. A marked processor gathers the results of all processors.

Table

4: Parallel PsDM with static partitioning.
of tasks, allocating to each of the P processor those tasks handling the sweep for
approximately dN=P e points. Table 4 outlines the method. The next question
is how to allocate tasks to processors. One natural idea is to use \static block
partitioning" and allocate to each processor the computations corresponding to
dN=P e consecutive points, say z
We applied PsDM on matrix kahan(100) from Section 4.2, starting from
points on the initial curve @ 0:1 (A) and computed
curves. Therefore, the total number of points that are computed by the end of
the run was In these and subsequent performance results we
did not take into account the time taken by Cobra to approximate the initial
curve.

Table

5 depicts the corresponding execution times and speedups. The
improvements are sunstantial and the speedups re
ect the parallel nature of the
# of Processors
Time (secs) 264 135 70 36

Table

5: Performance of parallel implementation of PsDM for kahan(100)).

Table

Number of points computed by each processor by the time the rst
processor nishes its share of the workload under static block partitioning.
processors are used to compute curves for kahan(100).
Numbers in boldface denote the number of points computed by the processor
that nished rst.
algorithm. In order to better understand the eect of the static task allocation
policy, we marked which processor would nish rst and then examined how
much work had been accomplished by then in each of the remaining processors.
The results are shown in Table 6, where in boldface are shown the number of
points that have been completed by the processor that nished rst. Line 1, for
instance, shows that processor 0 nished rst (the number 600 is in boldface)
while, by that time, processor 1 had accomplished the computation of 591 points.
Since each processor had to deal with this meant that
processor 1 still had to accomplish work for 9 points. As each row shows, the
work is reasonably well balanced, thus justifying the good speedups reported in

Table

5.
We next applied PsDM on matrix gre 1107 (1107  1107 sparse, real and
unsymmetric) from the Harwell-Boeing collection. The initial curve @ 0:1 (A)
was approximated by 64 points computed by Cobra. We computed 20 pseu-
dospectrum curves corresponding to log therefore each
processor was allocated d64  20=P e consecutive points. Rows 3 and 4 of Table
7 depict the times and speedups. Even though computing time is reduced, the
speedups are far inferior than those reported for kahan in Table 5. To explain
this phenomenon, we prepared for gre 1107 a table similar to Table 6. Results
are tabulated in rows 3-5 of Table 8 and reveal severe load imbalance. For
example, notice that when using 2 processors in PsDM, (2nd row of Table 8),
by the time processor 1 has nished with all its allocated points, processor 0
has nished with only 75% them. Similar patterns hold when using 4 and 8
processors. This load imbalance is due to the varying level of di-culty that an
iterative method, such as ARPACK, has when computing the triplet corresponding
to  min (zI A) as z moves from point to point.
Two ways to resolve this problem are a) a system-level approach to dispatch
Static block partitioning
Time (secs) 6000 3750 2500 1390
4.3
Static cyclic partitioning
Time (secs) 6000 3060 1560 840

Table

7: Performance of parallel implementation of PsDM using static block and
cyclic partitionings for gre 1107.
proc. / proc. id.
Static block partitioning
Static cyclic partitioning

Table

8: As Table 6 for matrix gre 1107 using static block and cyclic partitionings
to compute 20 pseudospectrum cuves.
tasks to the processors from a queue as processors are freed, and/or b) a problem-
level approach in which we estimate the work involved in each task and partition
the tasks so as to achieve acceptable load balance. The former approach has the
potential for better load balance, at additional system-level overhead; the latter
approach has little overhead and is simpler to implement, but appears to require
a priori estimates for the workload of each task. Even though this is di-cult
to know beforehand when using iterative methods, it helps to note that when
interested in load balancing, what is important is not information regarding the
amount of time taken by each task, but information regarding work dierential
between dierent tasks. Based on this idea, we tried the heuristic that the
number of iterations required for a triplet based at z is likely to be similar
for neighboring values of z. It is therefore reasonable to allocate points in an
interleaved fashion. This leads to \static cyclic partitioning", in which processor
are initialized with the points z Therefore, by
interleaving the initial points we aim to shue the \di-cult" points, along the
descent of SpDM and distribute them to all the processors. Rows 6-7 of Table
7 and rows 7-9 of Table 8 depict the speedups and load distribution obtained
from static cyclic partitioning. The improvement in load balance is evident and
leads to much better speedups.
We also experimented with \static block cyclic partitioning" in which points
were partitioned in dN=be blocks of b consecutive points each and then the blocks
were allocated in a cyclic fashion. We experimented with blocks of size 2 and 4
and found that performance was becoming worse with increasing blocksize, and
were all inferior than the static cyclic case.
We nally show, in Figure 8, the contours obtained for gre 1107.
Y

Figure

8: gre 1107. @  (A) contours for log
4.4 Adapting to decreasing contour lengths
It is known that for a given matrix and varying values of , the -pseudospectrum
forms a family of nested sets on the complex plane. Consequently, for any
given  and smaller -, the area of  - (A) is likely to be smaller than the area
of   (A); similarly, the length of the boundary is likely to change; it would
typically become smaller, unless there is separation and creation of disconnected
components whose total perimeter exceeds that of the original curve, in which
case it might increase. Let us assume now that any contour is approximated
by the polygonal path dened by the points computed by PsDM. Then we can
readily compute the lengths of the approximating paths. See for example the
lengths corresponding to dierent values of  for matrix kahan in Figure 6 and
the lengths of the corresponding polygonal paths in Table 9 (row 2). It is clear,
log
approx. lengths 4.48 2.63 1.85 1.48 1.29 1.18 1.1
# points 46 29 24 22 22 21 20

Table

9: Approximate contour lengths of the pseudospectra of kahan(100) and
number of points used to approximate each curve following the point reduction
policy described in this section.
in this case, that the polygonal path lengths become smaller as we move inwards.
Nevertheless, in the PsDM algorithm we described so far, the same number of
points (46) was used to approximate both curves.
It would seem appropriate, then, to monitor any signicant length reduction
or increase, and adapt accordingly the number of points computed by PsDM.
In the remainder of this section we examine the case of reduced path lengths,
noting that similar policies could also be adapted if we needed to increase, rather
than reduce, the number of points. The following general scheme describes the
transition from @ -m (A) to @ -m+1 (A). We denote by N -m the points dening
the contour @ -m (A).
1. Exclude K points of @ -m (A) according to a curve length criterion.
2. Proceed to compute the next contour @ -m+1 (A); - m+1 < - m , starting
from the N -m K remaining points of @ -m (A).
The scheme allows one to use a variety of curve length criteria, possibly chosen
dynamically as PsDM proceeds. In Table 10, we show a single step of PsDM
together with the implementation of such a strategy, based on the variation
of the minimum distance between consecutive points on subsequent contours.
Steps 1 to 6 implement the point reduction and also compute the next minimum
distance, l m+1 , between consecutive points of @ -m (A), while step 7 is the standard
sweep described in Table 2. First, the algorithm computes the distances
consecutive points on @ -m (A), and then drops any
point z k from the curve if it nds that d k + d k+1 < 2l m , unless the previous
point z k 1 had already been dropped. In the algorithm presented in Table 10,
this test is implemented by means of the boolean variables  k ;  k and a boolean
recurrence specied in line 4.1. Our implementation keeps xed the rst and
last points on @ -m (A), though this is easily modied. We also take advantage
of the fact that A is real so that we need to compute only one half of the curve.

Figure

9 illustrates the underlying idea. In the top curve, all points, z k ,
inside a box or circle are candidates for dropping, because they satisfy jz k 1
z (inside circle) are
dropped because they are the only ones for which the preceding point, z k 1 , was
not dropped. The remaining points are shown in the bottom curve, renumbered
to indicate their relative position within the curve.
distance between consecutive points on @ -m 1 (A);
Points
distance between consecutive points on @ -m (A);
Points
1. for
2. for
3.
4. for
4.1
4.2 if (
4.3
5. N
6. l
z
7. Apply single sweep of PsDM on f~z i 2 @ -m (A);
to compute fw g.

Table

10: Single sweep of PsDM with adaptive point reduction.
We implemented PsDM with the point reduction policy described in Table
on our parallel platform. We specically required the parallel version to
produce the same points that would have been produced, had the above policy
ran serially. To achieve this, processors synchronize so that the sweep outlined
in

Table

starts only when all input data is available to all. Due to the low
complexity of this procedure relative to the remaining work, we decided to use
a simple approach in which steps 1 to 6 are executed by a single, master, processor
instead of a more complicated parallel reduction policy. In particular,
processors send back points to the master processor who applies point reduc-
tion. The new points are then ready for allocation to the processors and the
application of PsDM (step 7 of Table 10). The question arises then, how to allocate
these points. As before, we could use a queue or a static approach. Given
the success of the latter in the previous experiments, and its lack of overhead, we
experimented with static approaches described earlier, assigning approximately
e points to each processor.
We rst applied this adaptive version of PsDM on kahan(100). This almost
halved the number of triplet evaluations, reducing them from 2760 to 1466.

Table

9 (3d row) depicts how the number of points varied for selected values
of . The contours computed using this strategy are depicted in Figure 10 vs.

Figure

9: Adaptive point reduction: Inscribed points are candidates for dropping
but only the encircled ones drop successfully.
kahan gre 1107
static block static cyclic static block static cyclic

Table

Execution times in seconds and speedups (in parentheses) for parallel
PsDM with point reduction for matrices kahan(100) (left) and gre 1107 (right).
the same contours computed using PsDM but holding the number of points per
contour constant. It is clear that adaptation is very eective without visibly
aecting the quality of the output. Columns 2 and 3 of Table 11 depict the
execution times and speedups (in parentheses). Once again, static cyclic is
superior to block assignment. Furthermore, the corresponding speedups are
similar, though, naturally, not as high as PsDM without point reduction (cf.

Table

5).
We also performed the same experiment on matrix gre 1107 and show the
times and corresponding speedups in columns 4 and 5 of Table 11. We started
with 74 points on @ 0:1 (A) and computed 20 curves as deep as @ 0:001 (A),
where we concluded with points. The total number of triplet evaluations was
reduced to 1300 compared to that would have been required if no
point reduction had been performed. In this case too, cyclic allocation performs
better than block. Furthermore, the achieved speedups are satisfactory, as they
are very close to those reported for the no drop policy in rows 7-9 of Table 8.
The above discussion shows that PsDM can be enhanced to adapt to the geo-
-0.4
Y
PsDM
REDUCED PsDM

Figure

10: Pseudospectra contours, log
Up: normal PsDM. Down: PsDM with point reduction.
metric features of the contours, in which case there is a corresponding reduction
in cost as we move inwards. It is clear that this idea can serve as springboard
for the design of alternative adaptation strategies.
4.5 Capturing disconnected components
One documented weakness of path following methods is that they cannot readily
capture disconnected components of the pseudospectrum. One question is how
does PsDM handle this di-culty? In this section we show that PsDM manages to
trace disconnected curves, as long as they lie inside the initial curve @  (A).
In that respect, the performance is the same as that of GRID, where the level
curves traced are only those that lie within the area discretized
by
h . Consider
matrix grcar(50) from the Test Matrix Toolbox. We begin from 192 points
proceed with a xed step 0.1 to compute curves. Thus
the inner curve is @ 10 6(A). The pseudospectra for disconnected
components. Figure 11 demonstrates that PsDM manages to retrieve
Y

Figure

Capturing disconnected components for grcar(50). Curves from
outer to inner:
this intricate structure. Of course this procedure does not return components
that were already outside the initial @  (A) curve. If we want to assure that all
components are captured, we can start with an initial curve that is large enough
using techniques such as those presented in [4].
Conclusions
We presented PsDM, a new method for the computation of the pseudospectrum
of a matrix A, that combines the appealing characteristics of the traditional
method GRID and the versatility of path following methods. We saw that PsDM
automatically generates new curves starting, for instance, from one application
of PF. The method is such that the curves adapt to the geometric properties of
the pseudospectrum and is able to capture disconnected components. Experimental
results showed that the method achieves a signicant reduction of the
number of necessary triplet evaluations vs. GRID even though, like GRID, it also
computes the pseudospectrum for several values of . Given the parallel nature
of PsDM, we implemented it using MPI. These experiments also revealed the advantages
of a simple heuristic designed to achieve load balance. An OpenMP
implementation of PsDM is currently underway and is expected to serve as a platform
for the investigation of a variety of alternative policies. It is worth noting
that PsDM computes successively, using path following but in the direction of
steepest descent, points dening the nested curves @  (A). This idea holds great
promise for pseudospectra and eigenvalues as described in [3, 9, 10]. Overall,
we believe that the approach used in PsDM will be useful in the construction of
an adaptive algorithm for the computation of pseudospectra that will be based
on path following. We note that the codes used in this paper are available from
URL http://www.hpclab.ceid.upatras.gr/scgroup/pseudospectra.html.

Acknowledgments

We thank Bernard Philippe, for his astute comments regarding the paper. We also
thank the referees for their valuable suggestions that helped improve the paper as well
as our colleague, E Kokiopoulou, for her comments and support.



--R

Numerical Continuation Methods: An Introduction
Cobra: Parallel path following for computing the matrix pseudospectrum.
Parallel algorithms for the computation of pseudospectra.
Using the

An algorithm for computing the distance to uncontrollability.
The Test Matrix Toolbox for MATLAB (version 3.0).

Iterations on domains for computing the matrix (pseudo) spectrum.
Hermitian methods for computing eigenvalues.
Arpack User's Guide: Solution of Large-Scale Eigenvalue Problems With Implicitly Restarted Arnoldi Methods
Directional Newton methods in n variables.

SPARSKIT: A basic tool-kit for sparse matrix computations (version 2)
Computation of pseudospectra.
--TR
Numerical continuation methods: an introduction
An algorithm for computing the distance to uncontrollability
Directional Newton Methods in n Variables

--CTR
C. Bekas , E. Kokiopoulou , E. Gallopoulos, The design of a distributed MATLAB-based environment for computing pseudospectra, Future Generation Computer Systems, v.21 n.6, p.930-941, June 2005
pseudospectra;parallel computation;newton's method;ARPACK