dlaebz.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 574 行 · 第 1/3 页
HTML
574 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>dlaebz.f</title>
<meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!--
body { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default a { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.comment { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.comment a { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
--></style>
</head>
<body>
<pre>
SUBROUTINE <a name="DLAEBZ.1"></a><a href="dlaebz.f.html#DLAEBZ.1">DLAEBZ</a>( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
$ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
$ NAB, WORK, IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK auxiliary routine (version 3.1) --
</span><span class="comment">*</span><span class="comment"> Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment"> November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Scalar Arguments ..
</span> INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * )
DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
$ WORK( * )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Purpose
</span><span class="comment">*</span><span class="comment"> =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> <a name="DLAEBZ.22"></a><a href="dlaebz.f.html#DLAEBZ.1">DLAEBZ</a> contains the iteration loops which compute and use the
</span><span class="comment">*</span><span class="comment"> function N(w), which is the count of eigenvalues of a symmetric
</span><span class="comment">*</span><span class="comment"> tridiagonal matrix T less than or equal to its argument w. It
</span><span class="comment">*</span><span class="comment"> performs a choice of two types of loops:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IJOB=1, followed by
</span><span class="comment">*</span><span class="comment"> IJOB=2: It takes as input a list of intervals and returns a list of
</span><span class="comment">*</span><span class="comment"> sufficiently small intervals whose union contains the same
</span><span class="comment">*</span><span class="comment"> eigenvalues as the union of the original intervals.
</span><span class="comment">*</span><span class="comment"> The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
</span><span class="comment">*</span><span class="comment"> The output interval (AB(j,1),AB(j,2)] will contain
</span><span class="comment">*</span><span class="comment"> eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IJOB=3: It performs a binary search in each input interval
</span><span class="comment">*</span><span class="comment"> (AB(j,1),AB(j,2)] for a point w(j) such that
</span><span class="comment">*</span><span class="comment"> N(w(j))=NVAL(j), and uses C(j) as the starting point of
</span><span class="comment">*</span><span class="comment"> the search. If such a w(j) is found, then on output
</span><span class="comment">*</span><span class="comment"> AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output
</span><span class="comment">*</span><span class="comment"> (AB(j,1),AB(j,2)] will be a small interval containing the
</span><span class="comment">*</span><span class="comment"> point where N(w) jumps through NVAL(j), unless that point
</span><span class="comment">*</span><span class="comment"> lies outside the initial interval.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Note that the intervals are in all cases half-open intervals,
</span><span class="comment">*</span><span class="comment"> i.e., of the form (a,b] , which includes b but not a .
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> To avoid underflow, the matrix should be scaled so that its largest
</span><span class="comment">*</span><span class="comment"> element is no greater than overflow**(1/2) * underflow**(1/4)
</span><span class="comment">*</span><span class="comment"> in absolute value. To assure the most accurate computation
</span><span class="comment">*</span><span class="comment"> of small eigenvalues, the matrix should be scaled to be
</span><span class="comment">*</span><span class="comment"> not much smaller than that, either.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
</span><span class="comment">*</span><span class="comment"> Matrix", Report CS41, Computer Science Dept., Stanford
</span><span class="comment">*</span><span class="comment"> University, July 21, 1966
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Note: the arguments are, in general, *not* checked for unreasonable
</span><span class="comment">*</span><span class="comment"> values.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Arguments
</span><span class="comment">*</span><span class="comment"> =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IJOB (input) INTEGER
</span><span class="comment">*</span><span class="comment"> Specifies what is to be done:
</span><span class="comment">*</span><span class="comment"> = 1: Compute NAB for the initial intervals.
</span><span class="comment">*</span><span class="comment"> = 2: Perform bisection iteration to find eigenvalues of T.
</span><span class="comment">*</span><span class="comment"> = 3: Perform bisection iteration to invert N(w), i.e.,
</span><span class="comment">*</span><span class="comment"> to find a point which has a specified number of
</span><span class="comment">*</span><span class="comment"> eigenvalues of T to its left.
</span><span class="comment">*</span><span class="comment"> Other values will cause <a name="DLAEBZ.70"></a><a href="dlaebz.f.html#DLAEBZ.1">DLAEBZ</a> to return with INFO=-1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> NITMAX (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The maximum number of "levels" of bisection to be
</span><span class="comment">*</span><span class="comment"> performed, i.e., an interval of width W will not be made
</span><span class="comment">*</span><span class="comment"> smaller than 2^(-NITMAX) * W. If not all intervals
</span><span class="comment">*</span><span class="comment"> have converged after NITMAX iterations, then INFO is set
</span><span class="comment">*</span><span class="comment"> to the number of non-converged intervals.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> N (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The dimension n of the tridiagonal matrix T. It must be at
</span><span class="comment">*</span><span class="comment"> least 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> MMAX (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The maximum number of intervals. If more than MMAX intervals
</span><span class="comment">*</span><span class="comment"> are generated, then <a name="DLAEBZ.85"></a><a href="dlaebz.f.html#DLAEBZ.1">DLAEBZ</a> will quit with INFO=MMAX+1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> MINP (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The initial number of intervals. It may not be greater than
</span><span class="comment">*</span><span class="comment"> MMAX.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> NBMIN (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The smallest number of intervals that should be processed
</span><span class="comment">*</span><span class="comment"> using a vector loop. If zero, then only the scalar loop
</span><span class="comment">*</span><span class="comment"> will be used.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ABSTOL (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> The minimum (absolute) width of an interval. When an
</span><span class="comment">*</span><span class="comment"> interval is narrower than ABSTOL, or than RELTOL times the
</span><span class="comment">*</span><span class="comment"> larger (in magnitude) endpoint, then it is considered to be
</span><span class="comment">*</span><span class="comment"> sufficiently small, i.e., converged. This must be at least
</span><span class="comment">*</span><span class="comment"> zero.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> RELTOL (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> The minimum relative width of an interval. When an interval
</span><span class="comment">*</span><span class="comment"> is narrower than ABSTOL, or than RELTOL times the larger (in
</span><span class="comment">*</span><span class="comment"> magnitude) endpoint, then it is considered to be
</span><span class="comment">*</span><span class="comment"> sufficiently small, i.e., converged. Note: this should
</span><span class="comment">*</span><span class="comment"> always be at least radix*machine epsilon.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> PIVMIN (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> The minimum absolute value of a "pivot" in the Sturm
</span><span class="comment">*</span><span class="comment"> sequence loop. This *must* be at least max |e(j)**2| *
</span><span class="comment">*</span><span class="comment"> safe_min and at least safe_min, where safe_min is at least
</span><span class="comment">*</span><span class="comment"> the smallest number that can divide one without overflow.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> D (input) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment"> The diagonal elements of the tridiagonal matrix T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> E (input) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment"> The offdiagonal elements of the tridiagonal matrix T in
</span><span class="comment">*</span><span class="comment"> positions 1 through N-1. E(N) is arbitrary.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> E2 (input) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment"> The squares of the offdiagonal elements of the tridiagonal
</span><span class="comment">*</span><span class="comment"> matrix T. E2(N) is ignored.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> NVAL (input/output) INTEGER array, dimension (MINP)
</span><span class="comment">*</span><span class="comment"> If IJOB=1 or 2, not referenced.
</span><span class="comment">*</span><span class="comment"> If IJOB=3, the desired values of N(w). The elements of NVAL
</span><span class="comment">*</span><span class="comment"> will be reordered to correspond with the intervals in AB.
</span><span class="comment">*</span><span class="comment"> Thus, NVAL(j) on output will not, in general be the same as
</span><span class="comment">*</span><span class="comment"> NVAL(j) on input, but it will correspond with the interval
</span><span class="comment">*</span><span class="comment"> (AB(j,1),AB(j,2)] on output.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> AB (input/output) DOUBLE PRECISION array, dimension (MMAX,2)
</span><span class="comment">*</span><span class="comment"> The endpoints of the intervals. AB(j,1) is a(j), the left
</span><span class="comment">*</span><span class="comment"> endpoint of the j-th interval, and AB(j,2) is b(j), the
</span><span class="comment">*</span><span class="comment"> right endpoint of the j-th interval. The input intervals
</span><span class="comment">*</span><span class="comment"> will, in general, be modified, split, and reordered by the
</span><span class="comment">*</span><span class="comment"> calculation.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> C (input/output) DOUBLE PRECISION array, dimension (MMAX)
</span><span class="comment">*</span><span class="comment"> If IJOB=1, ignored.
</span><span class="comment">*</span><span class="comment"> If IJOB=2, workspace.
</span><span class="comment">*</span><span class="comment"> If IJOB=3, then on input C(j) should be initialized to the
</span><span class="comment">*</span><span class="comment"> first search point in the binary search.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> MOUT (output) INTEGER
</span><span class="comment">*</span><span class="comment"> If IJOB=1, the number of eigenvalues in the intervals.
</span><span class="comment">*</span><span class="comment"> If IJOB=2 or 3, the number of intervals output.
</span><span class="comment">*</span><span class="comment"> If IJOB=3, MOUT will equal MINP.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> NAB (input/output) INTEGER array, dimension (MMAX,2)
</span><span class="comment">*</span><span class="comment"> If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
</span><span class="comment">*</span><span class="comment"> If IJOB=2, then on input, NAB(i,j) should be set. It must
</span><span class="comment">*</span><span class="comment"> satisfy the condition:
</span><span class="comment">*</span><span class="comment"> N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
</span><span class="comment">*</span><span class="comment"> which means that in interval i only eigenvalues
</span><span class="comment">*</span><span class="comment"> NAB(i,1)+1,...,NAB(i,2) will be considered. Usually,
</span><span class="comment">*</span><span class="comment"> NAB(i,j)=N(AB(i,j)), from a previous call to <a name="DLAEBZ.160"></a><a href="dlaebz.f.html#DLAEBZ.1">DLAEBZ</a> with
</span><span class="comment">*</span><span class="comment"> IJOB=1.
</span><span class="comment">*</span><span class="comment"> On output, NAB(i,j) will contain
</span><span class="comment">*</span><span class="comment"> max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
</span><span class="comment">*</span><span class="comment"> the input interval that the output interval
</span><span class="comment">*</span><span class="comment"> (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
</span><span class="comment">*</span><span class="comment"> the input values of NAB(k,1) and NAB(k,2).
</span><span class="comment">*</span><span class="comment"> If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
</span><span class="comment">*</span><span class="comment"> unless N(w) > NVAL(i) for all search points w , in which
</span><span class="comment">*</span><span class="comment"> case NAB(i,1) will not be modified, i.e., the output
</span><span class="comment">*</span><span class="comment"> value will be the same as the input value (modulo
</span><span class="comment">*</span><span class="comment"> reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
</span><span class="comment">*</span><span class="comment"> for all search points w , in which case NAB(i,2) will
</span><span class="comment">*</span><span class="comment"> not be modified. Normally, NAB should be set to some
</span><span class="comment">*</span><span class="comment"> distinctive value(s) before <a name="DLAEBZ.174"></a><a href="dlaebz.f.html#DLAEBZ.1">DLAEBZ</a> is called.
</span><span class="comment">*</span><span class="comment">
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?