slaebz.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 574 行 · 第 1/3 页
HTML
574 行
</span><span class="comment">*</span><span class="comment"> WORK (workspace) REAL array, dimension (MMAX)
</span><span class="comment">*</span><span class="comment"> Workspace.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IWORK (workspace) INTEGER array, dimension (MMAX)
</span><span class="comment">*</span><span class="comment"> Workspace.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> INFO (output) INTEGER
</span><span class="comment">*</span><span class="comment"> = 0: All intervals converged.
</span><span class="comment">*</span><span class="comment"> = 1--MMAX: The last INFO intervals did not converge.
</span><span class="comment">*</span><span class="comment"> = MMAX+1: More than MMAX intervals were generated.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Further Details
</span><span class="comment">*</span><span class="comment"> ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> This routine is intended to be called only by other LAPACK
</span><span class="comment">*</span><span class="comment"> routines, thus the interface is less user-friendly. It is intended
</span><span class="comment">*</span><span class="comment"> for two purposes:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> (a) finding eigenvalues. In this case, <a name="SLAEBZ.194"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> should have one or
</span><span class="comment">*</span><span class="comment"> more initial intervals set up in AB, and <a name="SLAEBZ.195"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> should be called
</span><span class="comment">*</span><span class="comment"> with IJOB=1. This sets up NAB, and also counts the eigenvalues.
</span><span class="comment">*</span><span class="comment"> Intervals with no eigenvalues would usually be thrown out at
</span><span class="comment">*</span><span class="comment"> this point. Also, if not all the eigenvalues in an interval i
</span><span class="comment">*</span><span class="comment"> are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
</span><span class="comment">*</span><span class="comment"> For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
</span><span class="comment">*</span><span class="comment"> eigenvalue. <a name="SLAEBZ.201"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> is then called with IJOB=2 and MMAX
</span><span class="comment">*</span><span class="comment"> no smaller than the value of MOUT returned by the call with
</span><span class="comment">*</span><span class="comment"> IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1
</span><span class="comment">*</span><span class="comment"> through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
</span><span class="comment">*</span><span class="comment"> tolerance specified by ABSTOL and RELTOL.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
</span><span class="comment">*</span><span class="comment"> In this case, start with a Gershgorin interval (a,b). Set up
</span><span class="comment">*</span><span class="comment"> AB to contain 2 search intervals, both initially (a,b). One
</span><span class="comment">*</span><span class="comment"> NVAL element should contain f-1 and the other should contain l
</span><span class="comment">*</span><span class="comment"> , while C should contain a and b, resp. NAB(i,1) should be -1
</span><span class="comment">*</span><span class="comment"> and NAB(i,2) should be N+1, to flag an error if the desired
</span><span class="comment">*</span><span class="comment"> interval does not lie in (a,b). <a name="SLAEBZ.213"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> is then called with
</span><span class="comment">*</span><span class="comment"> IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals --
</span><span class="comment">*</span><span class="comment"> j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
</span><span class="comment">*</span><span class="comment"> if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
</span><span class="comment">*</span><span class="comment"> >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and
</span><span class="comment">*</span><span class="comment"> N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and
</span><span class="comment">*</span><span class="comment"> w(l-r)=...=w(l+k) are handled similarly.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Parameters ..
</span> REAL ZERO, TWO, HALF
PARAMETER ( ZERO = 0.0E0, TWO = 2.0E0,
$ HALF = 1.0E0 / TWO )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,
$ KLNEW
REAL TMP1, TMP2
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC ABS, MAX, MIN
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Check for Errors
</span><span class="comment">*</span><span class="comment">
</span> INFO = 0
IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN
INFO = -1
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Initialize NAB
</span><span class="comment">*</span><span class="comment">
</span> IF( IJOB.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute the number of eigenvalues in the initial intervals.
</span><span class="comment">*</span><span class="comment">
</span> MOUT = 0
CDIR$ NOVECTOR
DO 30 JI = 1, MINP
DO 20 JP = 1, 2
TMP1 = D( 1 ) - AB( JI, JP )
IF( ABS( TMP1 ).LT.PIVMIN )
$ TMP1 = -PIVMIN
NAB( JI, JP ) = 0
IF( TMP1.LE.ZERO )
$ NAB( JI, JP ) = 1
<span class="comment">*</span><span class="comment">
</span> DO 10 J = 2, N
TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP )
IF( ABS( TMP1 ).LT.PIVMIN )
$ TMP1 = -PIVMIN
IF( TMP1.LE.ZERO )
$ NAB( JI, JP ) = NAB( JI, JP ) + 1
10 CONTINUE
20 CONTINUE
MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 )
30 CONTINUE
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Initialize for loop
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> KF and KL have the following meaning:
</span><span class="comment">*</span><span class="comment"> Intervals 1,...,KF-1 have converged.
</span><span class="comment">*</span><span class="comment"> Intervals KF,...,KL still need to be refined.
</span><span class="comment">*</span><span class="comment">
</span> KF = 1
KL = MINP
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> If IJOB=2, initialize C.
</span><span class="comment">*</span><span class="comment"> If IJOB=3, use the user-supplied starting point.
</span><span class="comment">*</span><span class="comment">
</span> IF( IJOB.EQ.2 ) THEN
DO 40 JI = 1, MINP
C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
40 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Iteration loop
</span><span class="comment">*</span><span class="comment">
</span> DO 130 JIT = 1, NITMAX
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Loop over intervals
</span><span class="comment">*</span><span class="comment">
</span> IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Begin of Parallel Version of the loop
</span><span class="comment">*</span><span class="comment">
</span> DO 60 JI = KF, KL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute N(c), the number of eigenvalues less than c
</span><span class="comment">*</span><span class="comment">
</span> WORK( JI ) = D( 1 ) - C( JI )
IWORK( JI ) = 0
IF( WORK( JI ).LE.PIVMIN ) THEN
IWORK( JI ) = 1
WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
END IF
<span class="comment">*</span><span class="comment">
</span> DO 50 J = 2, N
WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI )
IF( WORK( JI ).LE.PIVMIN ) THEN
IWORK( JI ) = IWORK( JI ) + 1
WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
END IF
50 CONTINUE
60 CONTINUE
<span class="comment">*</span><span class="comment">
</span> IF( IJOB.LE.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IJOB=2: Choose all intervals containing eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span> KLNEW = KL
DO 70 JI = KF, KL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Insure that N(w) is monotone
</span><span class="comment">*</span><span class="comment">
</span> IWORK( JI ) = MIN( NAB( JI, 2 ),
$ MAX( NAB( JI, 1 ), IWORK( JI ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Update the Queue -- add intervals if both halves
</span><span class="comment">*</span><span class="comment"> contain eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span> IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> No eigenvalue in the upper interval:
</span><span class="comment">*</span><span class="comment"> just use the lower interval.
</span><span class="comment">*</span><span class="comment">
</span> AB( JI, 2 ) = C( JI )
<span class="comment">*</span><span class="comment">
</span> ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> No eigenvalue in the lower interval:
</span><span class="comment">*</span><span class="comment"> just use the upper interval.
</span><span class="comment">*</span><span class="comment">
</span> AB( JI, 1 ) = C( JI )
ELSE
KLNEW = KLNEW + 1
IF( KLNEW.LE.MMAX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Eigenvalue in both intervals -- add upper to
</span><span class="comment">*</span><span class="comment"> queue.
</span><span class="comment">*</span><span class="comment">
</span> AB( KLNEW, 2 ) = AB( JI, 2 )
NAB( KLNEW, 2 ) = NAB( JI, 2 )
AB( KLNEW, 1 ) = C( JI )
NAB( KLNEW, 1 ) = IWORK( JI )
AB( JI, 2 ) = C( JI )
NAB( JI, 2 ) = IWORK( JI )
ELSE
INFO = MMAX + 1
END IF
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?