Партнёры:
Хостинг:
NAME
slaebz - contain the iteration loops which compute and use
the function N(w), which is the count of eigenvalues of a
symmetric tridiagonal matrix T less than or equal to its
argument w
SYNOPSIS
SUBROUTINE SLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN,
ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL, AB, C,
MOUT, NAB, WORK, IWORK, INFO )
INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
REAL ABSTOL, PIVMIN, RELTOL
INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * )
REAL AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), WORK( *
)
#include <sunperf.h>
void slaebz(int ijob, int nitmax, int n, int mmax, int minp,
int nbmin, float abstol, float reltol, float piv-
min, float *d, float *e, float *e2, int *nval,
float *ab, float *sc, int *mout, int *nab, int
*info) ;
PURPOSE
SLAEBZ contains the iteration loops which compute and use
the function N(w), which is the count of eigenvalues of a
symmetric tridiagonal matrix T less than or equal to its
argument w. It performs a choice of two types of loops:
IJOB=1, followed by
IJOB=2: It takes as input a list of intervals and returns a
list of
sufficiently small intervals whose union contains
the same
eigenvalues as the union of the original intervals.
The input intervals are (AB(j,1),AB(j,2)],
j=1,...,MINP.
The output interval (AB(j,1),AB(j,2)] will contain
eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <=
MOUT.
IJOB=3: It performs a binary search in each input interval
(AB(j,1),AB(j,2)] for a point w(j) such that
N(w(j))=NVAL(j), and uses C(j) as the starting
point of
the search. If such a w(j) is found, then on output
AB(j,1)=AB(j,2)=w. If no such w(j) is found, then
on output
(AB(j,1),AB(j,2)] will be a small interval contain-
ing the
point where N(w) jumps through NVAL(j), unless that
point
lies outside the initial interval.
Note that the intervals are in all cases half-open inter-
vals, i.e., of the form (a,b] , which includes b but not
a .
To avoid underflow, the matrix should be scaled so that its
largest element is no greater than overflow**(1/2) * under-
flow**(1/4) in absolute value. To assure the most accurate
computation of small eigenvalues, the matrix should be
scaled to be
not much smaller than that, either.
See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiago-
nal Matrix", Report CS41, Computer Science Dept., Stanford
University, July 21, 1966
Note: the arguments are, in general, *not* checked for
unreasonable values.
ARGUMENTS
IJOB (input) INTEGER
Specifies what is to be done:
= 1: Compute NAB for the initial intervals.
= 2: Perform bisection iteration to find eigen-
values of T.
= 3: Perform bisection iteration to invert N(w),
i.e., to find a point which has a specified number
of eigenvalues of T to its left. Other values
will cause SLAEBZ to return with INFO=-1.
NITMAX (input) INTEGER
The maximum number of "levels" of bisection to be
performed, i.e., an interval of width W will not
be made smaller than 2^(-NITMAX) * W. If not all
intervals have converged after NITMAX iterations,
then INFO is set to the number of non-converged
intervals.
N (input) INTEGER
The dimension n of the tridiagonal matrix T. It
must be at least 1.
MMAX (input) INTEGER
The maximum number of intervals. If more than
MMAX intervals are generated, then SLAEBZ will
quit with INFO=MMAX+1.
MINP (input) INTEGER
The initial number of intervals. It may not be
greater than MMAX.
NBMIN (input) INTEGER
The smallest number of intervals that should be
processed using a vector loop. If zero, then only
the scalar loop will be used.
ABSTOL (input) REAL
The minimum (absolute) width of an interval. When
an interval is narrower than ABSTOL, or than REL-
TOL times the larger (in magnitude) endpoint, then
it is considered to be sufficiently small, i.e.,
converged. This must be at least zero.
RELTOL (input) REAL
The minimum relative width of an interval. When
an interval is narrower than ABSTOL, or than REL-
TOL times the larger (in magnitude) endpoint, then
it is considered to be sufficiently small, i.e.,
converged. Note: this should always be at least
radix*machine epsilon.
PIVMIN (input) REAL
The minimum absolute value of a "pivot" in the
Sturm sequence loop. This *must* be at least max
|e(j)**2| * safe_min and at least safe_min, where
safe_min is at least the smallest number that can
divide one without overflow.
D (input) REAL array, dimension (N)
The diagonal elements of the tridiagonal matrix T.
E (input) REAL array, dimension (N)
The offdiagonal elements of the tridiagonal matrix
T in positions 1 through N-1. E(N) is arbitrary.
E2 (input) REAL array, dimension (N)
The squares of the offdiagonal elements of the
tridiagonal matrix T. E2(N) is ignored.
NVAL (input/output) INTEGER array, dimension (MINP)
If IJOB=1 or 2, not referenced. If IJOB=3, the
desired values of N(w). The elements of NVAL will
be reordered to correspond with the intervals in
AB. Thus, NVAL(j) on output will not, in general
be the same as NVAL(j) on input, but it will
correspond with the interval (AB(j,1),AB(j,2)] on
output.
AB (input/output) REAL array, dimension (MMAX,2)
The endpoints of the intervals. AB(j,1) is a(j),
the left endpoint of the j-th interval, and
AB(j,2) is b(j), the right endpoint of the j-th
interval. The input intervals will, in general,
be modified, split, and reordered by the calcula-
tion.
C (input/output) REAL array, dimension (MMAX)
If IJOB=1, ignored. If IJOB=2, workspace. If
IJOB=3, then on input C(j) should be initialized
to the first search point in the binary search.
MOUT (output) INTEGER
If IJOB=1, the number of eigenvalues in the inter-
vals. If IJOB=2 or 3, the number of intervals
output. If IJOB=3, MOUT will equal MINP.
NAB (input/output) INTEGER array, dimension (MMAX,2)
If IJOB=1, then on output NAB(i,j) will be set to
N(AB(i,j)). If IJOB=2, then on input, NAB(i,j)
should be set. It must satisfy the condition:
N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
which means that in interval i only eigenvalues
NAB(i,1)+1,...,NAB(i,2) will be considered. Usu-
ally, NAB(i,j)=N(AB(i,j)), from a previous call to
SLAEBZ with IJOB=1. On output, NAB(i,j) will con-
tain max(na(k),min(nb(k),N(AB(i,j)))), where k is
the index of the input interval that the output
interval (AB(j,1),AB(j,2)] came from, and na(k)
and nb(k) are the the input values of NAB(k,1) and
NAB(k,2). If IJOB=3, then on output, NAB(i,j)
contains N(AB(i,j)), unless N(w) > NVAL(i) for all
search points w , in which case NAB(i,1) will not
be modified, i.e., the output value will be the
same as the input value (modulo reorderings -- see
NVAL and AB), or unless N(w) < NVAL(i) for all
search points w , in which case NAB(i,2) will not
be modified. Normally, NAB should be set to some
distinctive value(s) before SLAEBZ is called.
WORK (workspace) REAL array, dimension (MMAX)
Workspace.
IWORK (workspace) INTEGER array, dimension (MMAX)
Workspace.
INFO (output) INTEGER
= 0: All intervals converged.
= 1--MMAX: The last INFO intervals did not con-
verge.
= MMAX+1: More than MMAX intervals were gen-
erated.
FURTHER DETAILS
This routine is intended to be called only by other
LAPACK routines, thus the interface is less user-friendly.
It is intended for two purposes:
(a) finding eigenvalues. In this case, SLAEBZ should have
one or
more initial intervals set up in AB, and SLAEBZ should
be called
with IJOB=1. This sets up NAB, and also counts the
eigenvalues.
Intervals with no eigenvalues would usually be thrown
out at
this point. Also, if not all the eigenvalues in an
interval i
are desired, NAB(i,1) can be increased or NAB(i,2)
decreased.
For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
eigenvalue. SLAEBZ is then called with IJOB=2 and MMAX
no smaller than the value of MOUT returned by the call
with
IJOB=1. After this (IJOB=2) call, eigenvalues
NAB(i,1)+1
through NAB(i,2) are approximately AB(i,1) (or AB(i,2))
to the
tolerance specified by ABSTOL and RELTOL.
(b) finding an interval (a',b'] containing eigenvalues
w(f),...,w(l).
In this case, start with a Gershgorin interval (a,b).
Set up
AB to contain 2 search intervals, both initially (a,b).
One
NVAL element should contain f-1 and the other should
contain l
, while C should contain a and b, resp. NAB(i,1) should
be -1
and NAB(i,2) should be N+1, to flag an error if the
desired
interval does not lie in (a,b). SLAEBZ is then called
with
IJOB=3. On exit, if w(f-1) < w(f), then one of the
intervals --
j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-
1, while
if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0
and r
>= 0, then the interval will have
N(AB(j,1))=NAB(j,1)=f-k and
N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and
w(l-r)=...=w(l+k) are handled similarly.
|
Закладки на сайте Проследить за страницей |
Created 1996-2024 by Maxim Chirkov Добавить, Поддержать, Вебмастеру |