Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F68242080
dlaed4.f
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Jun 26, 12:43
Size
26 KB
Mime Type
text/html
Expires
Fri, Jun 28, 12:43 (2 d)
Engine
blob
Format
Raw Data
Handle
18362583
Attached To
rLAMMPS lammps
dlaed4.f
View Options
*>
\
brief
\
b
DLAED4
used
by
sstedc
.
Finds
a
single
root
of
the
secular
equation
.
*
*
===========
DOCUMENTATION
===========
*
*
Online
html
documentation
available
at
*
http
:
//
www
.
netlib
.
org
/
lapack
/
explore
-
html
/
*
*>
\
htmlonly
*>
Download
DLAED4
+
dependencies
*>
<
a
href
=
"http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed4.f"
>
*>
[
TGZ
]
</
a
>
*>
<
a
href
=
"http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed4.f"
>
*>
[
ZIP
]
</
a
>
*>
<
a
href
=
"http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed4.f"
>
*>
[
TXT
]
</
a
>
*>
\
endhtmlonly
*
*
Definition
:
*
===========
*
*
SUBROUTINE
DLAED4
(
N
,
I
,
D
,
Z
,
DELTA
,
RHO
,
DLAM
,
INFO
)
*
*
..
Scalar
Arguments
..
*
INTEGER
I
,
INFO
,
N
*
DOUBLE PRECISION
DLAM
,
RHO
*
..
*
..
Array
Arguments
..
*
DOUBLE PRECISION
D
(
*
),
DELTA
(
*
),
Z
(
*
)
*
..
*
*
*>
\
par
Purpose
:
*
=============
*>
*>
\
verbatim
*>
*>
This
subroutine
computes
the
I
-
th
updated
eigenvalue
of
a
symmetric
*>
rank
-
one
modification
to
a
diagonal
matrix
whose
elements
are
*>
given
in
the
array
d
,
and
that
*>
*>
D
(
i
)
<
D
(
j
)
for
i
<
j
*>
*>
and
that
RHO
>
0.
This
is
arranged
by
the
calling
routine
,
and
is
*>
no
loss
in
generality
.
The
rank
-
one
modified
system
is
thus
*>
*>
diag
(
D
)
+
RHO
*
Z
*
Z_transpose
.
*>
*>
where
we
assume
the
Euclidean
norm
of
Z
is
1.
*>
*>
The
method
consists
of
approximating
the
rational
functions
in
the
*>
secular
equation
by
simpler
interpolating
rational
functions
.
*>
\
endverbatim
*
*
Arguments
:
*
==========
*
*>
\
param
[
in
]
N
*>
\
verbatim
*>
N
is
INTEGER
*>
The
length
of
all
arrays
.
*>
\
endverbatim
*>
*>
\
param
[
in
]
I
*>
\
verbatim
*>
I
is
INTEGER
*>
The
index
of
the
eigenvalue
to
be
computed
.
1
<=
I
<=
N
.
*>
\
endverbatim
*>
*>
\
param
[
in
]
D
*>
\
verbatim
*>
D
is
DOUBLE PRECISION
array
,
dimension
(
N
)
*>
The
original
eigenvalues
.
It
is
assumed
that
they
are
in
*>
order
,
D
(
I
)
<
D
(
J
)
for
I
<
J
.
*>
\
endverbatim
*>
*>
\
param
[
in
]
Z
*>
\
verbatim
*>
Z
is
DOUBLE PRECISION
array
,
dimension
(
N
)
*>
The
components
of
the
updating
vector
.
*>
\
endverbatim
*>
*>
\
param
[
out
]
DELTA
*>
\
verbatim
*>
DELTA
is
DOUBLE PRECISION
array
,
dimension
(
N
)
*>
If
N
.GT.
2
,
DELTA
contains
(
D
(
j
)
-
lambda_I
)
in
its
j
-
th
*>
component
.
If
N
=
1
,
then
DELTA
(
1
)
=
1.
If
N
=
2
,
see
DLAED5
*>
for
detail
.
The
vector
DELTA
contains
the
information
necessary
*>
to
construct
the
eigenvectors
by
DLAED3
and
DLAED9
.
*>
\
endverbatim
*>
*>
\
param
[
in
]
RHO
*>
\
verbatim
*>
RHO
is
DOUBLE PRECISION
*>
The
scalar
in
the
symmetric
updating
formula
.
*>
\
endverbatim
*>
*>
\
param
[
out
]
DLAM
*>
\
verbatim
*>
DLAM
is
DOUBLE PRECISION
*>
The
computed
lambda_I
,
the
I
-
th
updated
eigenvalue
.
*>
\
endverbatim
*>
*>
\
param
[
out
]
INFO
*>
\
verbatim
*>
INFO
is
INTEGER
*>
=
0
:
successful
exit
*>
>
0
:
if
INFO
=
1
,
the
updating
process
failed
.
*>
\
endverbatim
*
*>
\
par
Internal
Parameters
:
*
=========================
*>
*>
\
verbatim
*>
Logical
variable
ORGATI
(
origin
-
at
-
i
?
)
is
used
for
distinguishing
*>
whether
D
(
i
)
or
D
(
i
+
1
)
is
treated
as
the
origin
.
*>
*>
ORGATI
=
.true.
origin
at
i
*>
ORGATI
=
.false.
origin
at
i
+
1
*>
*>
Logical
variable
SWTCH3
(
switch
-
for
-
3
-
poles
?
)
is
for
noting
*>
if
we
are
working
with
THREE
poles
!
*>
*>
MAXIT
is
the
maximum
number
of
iterations
allowed
for
each
*>
eigenvalue
.
*>
\
endverbatim
*
*
Authors
:
*
========
*
*>
\
author
Univ
.
of
Tennessee
*>
\
author
Univ
.
of
California
Berkeley
*>
\
author
Univ
.
of
Colorado
Denver
*>
\
author
NAG
Ltd
.
*
*>
\
date
September
2012
*
*>
\
ingroup
auxOTHERcomputational
*
*>
\
par
Contributors
:
*
==================
*>
*>
Ren
-
Cang
Li
,
Computer
Science
Division
,
University
of
California
*>
at
Berkeley
,
USA
*>
*
=====================================================================
SUBROUTINE
DLAED4
(
N
,
I
,
D
,
Z
,
DELTA
,
RHO
,
DLAM
,
INFO
)
*
*
--
LAPACK
computational
routine
(
version
3.4.2
)
--
*
--
LAPACK
is
a
software
package
provided
by
Univ
.
of
Tennessee
,
--
*
--
Univ
.
of
California
Berkeley
,
Univ
.
of
Colorado
Denver
and
NAG
Ltd
..
--
*
September
2012
*
*
..
Scalar
Arguments
..
INTEGER
I
,
INFO
,
N
DOUBLE PRECISION
DLAM
,
RHO
*
..
*
..
Array
Arguments
..
DOUBLE PRECISION
D
(
*
),
DELTA
(
*
),
Z
(
*
)
*
..
*
*
=====================================================================
*
*
..
Parameters
..
INTEGER
MAXIT
PARAMETER
(
MAXIT
=
30
)
DOUBLE PRECISION
ZERO
,
ONE
,
TWO
,
THREE
,
FOUR
,
EIGHT
,
TEN
PARAMETER
(
ZERO
=
0.0
D0
,
ONE
=
1.0
D0
,
TWO
=
2.0
D0
,
$
THREE
=
3.0
D0
,
FOUR
=
4.0
D0
,
EIGHT
=
8.0
D0
,
$
TEN
=
1
0.0
D0
)
*
..
*
..
Local
Scalars
..
LOGICAL
ORGATI
,
SWTCH
,
SWTCH3
INTEGER
II
,
IIM1
,
IIP1
,
IP1
,
ITER
,
J
,
NITER
DOUBLE PRECISION
A
,
B
,
C
,
DEL
,
DLTLB
,
DLTUB
,
DPHI
,
DPSI
,
DW
,
$
EPS
,
ERRETM
,
ETA
,
MIDPT
,
PHI
,
PREW
,
PSI
,
$
RHOINV
,
TAU
,
TEMP
,
TEMP1
,
W
*
..
*
..
Local
Arrays
..
DOUBLE PRECISION
ZZ
(
3
)
*
..
*
..
External
Functions
..
DOUBLE PRECISION
DLAMCH
EXTERNAL
DLAMCH
*
..
*
..
External
Subroutines
..
EXTERNAL
DLAED5
,
DLAED6
*
..
*
..
Intrinsic
Functions
..
INTRINSIC
ABS
,
MAX
,
MIN
,
SQRT
*
..
*
..
Executable
Statements
..
*
*
Since
this
routine
is
called
in
an
inner
loop
,
we
do
no
argument
*
checking
.
*
*
Quick
return
for
N
=
1
and
2.
*
INFO
=
0
IF
(
N
.EQ.
1
)
THEN
*
*
Presumably
,
I
=
1
upon
entry
*
DLAM
=
D
(
1
)
+
RHO
*
Z
(
1
)
*
Z
(
1
)
DELTA
(
1
)
=
ONE
RETURN
END IF
IF
(
N
.EQ.
2
)
THEN
CALL
DLAED5
(
I
,
D
,
Z
,
DELTA
,
RHO
,
DLAM
)
RETURN
END IF
*
*
Compute
machine
epsilon
*
EPS
=
DLAMCH
(
'Epsilon'
)
RHOINV
=
ONE
/
RHO
*
*
The
case
I
=
N
*
IF
(
I
.EQ.
N
)
THEN
*
*
Initialize
some
basic
variables
*
II
=
N
-
1
NITER
=
1
*
*
Calculate
initial
guess
*
MIDPT
=
RHO
/
TWO
*
*
If
||
Z
||_
2
is
not
one
,
then
TEMP
should
be
set
to
*
RHO
*
||
Z
||_
2
^
2
/
TWO
*
DO
10
J
=
1
,
N
DELTA
(
J
)
=
(
D
(
J
)
-
D
(
I
)
)
-
MIDPT
10
CONTINUE
*
PSI
=
ZERO
DO
20
J
=
1
,
N
-
2
PSI
=
PSI
+
Z
(
J
)
*
Z
(
J
)
/
DELTA
(
J
)
20
CONTINUE
*
C
=
RHOINV
+
PSI
W
=
C
+
Z
(
II
)
*
Z
(
II
)
/
DELTA
(
II
)
+
$
Z
(
N
)
*
Z
(
N
)
/
DELTA
(
N
)
*
IF
(
W
.LE.
ZERO
)
THEN
TEMP
=
Z
(
N
-
1
)
*
Z
(
N
-
1
)
/
(
D
(
N
)
-
D
(
N
-
1
)
+
RHO
)
+
$
Z
(
N
)
*
Z
(
N
)
/
RHO
IF
(
C
.LE.
TEMP
)
THEN
TAU
=
RHO
ELSE
DEL
=
D
(
N
)
-
D
(
N
-
1
)
A
=
-
C
*
DEL
+
Z
(
N
-
1
)
*
Z
(
N
-
1
)
+
Z
(
N
)
*
Z
(
N
)
B
=
Z
(
N
)
*
Z
(
N
)
*
DEL
IF
(
A
.LT.
ZERO
)
THEN
TAU
=
TWO
*
B
/
(
SQRT
(
A
*
A
+
FOUR
*
B
*
C
)
-
A
)
ELSE
TAU
=
(
A
+
SQRT
(
A
*
A
+
FOUR
*
B
*
C
)
)
/
(
TWO
*
C
)
END IF
END IF
*
*
It
can
be
proved
that
*
D
(
N
)
+
RHO
/
2
<=
LAMBDA
(
N
)
<
D
(
N
)
+
TAU
<=
D
(
N
)
+
RHO
*
DLTLB
=
MIDPT
DLTUB
=
RHO
ELSE
DEL
=
D
(
N
)
-
D
(
N
-
1
)
A
=
-
C
*
DEL
+
Z
(
N
-
1
)
*
Z
(
N
-
1
)
+
Z
(
N
)
*
Z
(
N
)
B
=
Z
(
N
)
*
Z
(
N
)
*
DEL
IF
(
A
.LT.
ZERO
)
THEN
TAU
=
TWO
*
B
/
(
SQRT
(
A
*
A
+
FOUR
*
B
*
C
)
-
A
)
ELSE
TAU
=
(
A
+
SQRT
(
A
*
A
+
FOUR
*
B
*
C
)
)
/
(
TWO
*
C
)
END IF
*
*
It
can
be
proved
that
*
D
(
N
)
<
D
(
N
)
+
TAU
<
LAMBDA
(
N
)
<
D
(
N
)
+
RHO
/
2
*
DLTLB
=
ZERO
DLTUB
=
MIDPT
END IF
*
DO
30
J
=
1
,
N
DELTA
(
J
)
=
(
D
(
J
)
-
D
(
I
)
)
-
TAU
30
CONTINUE
*
*
Evaluate
PSI
and
the
derivative
DPSI
*
DPSI
=
ZERO
PSI
=
ZERO
ERRETM
=
ZERO
DO
40
J
=
1
,
II
TEMP
=
Z
(
J
)
/
DELTA
(
J
)
PSI
=
PSI
+
Z
(
J
)
*
TEMP
DPSI
=
DPSI
+
TEMP
*
TEMP
ERRETM
=
ERRETM
+
PSI
40
CONTINUE
ERRETM
=
ABS
(
ERRETM
)
*
*
Evaluate
PHI
and
the
derivative
DPHI
*
TEMP
=
Z
(
N
)
/
DELTA
(
N
)
PHI
=
Z
(
N
)
*
TEMP
DPHI
=
TEMP
*
TEMP
ERRETM
=
EIGHT
*
(
-
PHI
-
PSI
)
+
ERRETM
-
PHI
+
RHOINV
+
$
ABS
(
TAU
)
*
(
DPSI
+
DPHI
)
*
W
=
RHOINV
+
PHI
+
PSI
*
*
Test
for
convergence
*
IF
(
ABS
(
W
)
.LE.
EPS
*
ERRETM
)
THEN
DLAM
=
D
(
I
)
+
TAU
GO
TO
250
END IF
*
IF
(
W
.LE.
ZERO
)
THEN
DLTLB
=
MAX
(
DLTLB
,
TAU
)
ELSE
DLTUB
=
MIN
(
DLTUB
,
TAU
)
END IF
*
*
Calculate
the
new
step
*
NITER
=
NITER
+
1
C
=
W
-
DELTA
(
N
-
1
)
*
DPSI
-
DELTA
(
N
)
*
DPHI
A
=
(
DELTA
(
N
-
1
)
+
DELTA
(
N
)
)
*
W
-
$
DELTA
(
N
-
1
)
*
DELTA
(
N
)
*
(
DPSI
+
DPHI
)
B
=
DELTA
(
N
-
1
)
*
DELTA
(
N
)
*
W
IF
(
C
.LT.
ZERO
)
$
C
=
ABS
(
C
)
IF
(
C
.EQ.
ZERO
)
THEN
*
ETA
=
B
/
A
*
ETA
=
RHO
-
TAU
ETA
=
DLTUB
-
TAU
ELSE IF
(
A
.GE.
ZERO
)
THEN
ETA
=
(
A
+
SQRT
(
ABS
(
A
*
A
-
FOUR
*
B
*
C
)
)
)
/
(
TWO
*
C
)
ELSE
ETA
=
TWO
*
B
/
(
A
-
SQRT
(
ABS
(
A
*
A
-
FOUR
*
B
*
C
)
)
)
END IF
*
*
Note
,
eta
should
be
positive
if
w
is
negative
,
and
*
eta
should
be
negative
otherwise
.
However
,
*
if
for
some
reason
caused
by
roundoff
,
eta
*
w
>
0
,
*
we
simply
use
one
Newton
step
instead
.
This
way
*
will
guarantee
eta
*
w
<
0.
*
IF
(
W
*
ETA
.GT.
ZERO
)
$
ETA
=
-
W
/
(
DPSI
+
DPHI
)
TEMP
=
TAU
+
ETA
IF
(
TEMP
.GT.
DLTUB
.OR.
TEMP
.LT.
DLTLB
)
THEN
IF
(
W
.LT.
ZERO
)
THEN
ETA
=
(
DLTUB
-
TAU
)
/
TWO
ELSE
ETA
=
(
DLTLB
-
TAU
)
/
TWO
END IF
END IF
DO
50
J
=
1
,
N
DELTA
(
J
)
=
DELTA
(
J
)
-
ETA
50
CONTINUE
*
TAU
=
TAU
+
ETA
*
*
Evaluate
PSI
and
the
derivative
DPSI
*
DPSI
=
ZERO
PSI
=
ZERO
ERRETM
=
ZERO
DO
60
J
=
1
,
II
TEMP
=
Z
(
J
)
/
DELTA
(
J
)
PSI
=
PSI
+
Z
(
J
)
*
TEMP
DPSI
=
DPSI
+
TEMP
*
TEMP
ERRETM
=
ERRETM
+
PSI
60
CONTINUE
ERRETM
=
ABS
(
ERRETM
)
*
*
Evaluate
PHI
and
the
derivative
DPHI
*
TEMP
=
Z
(
N
)
/
DELTA
(
N
)
PHI
=
Z
(
N
)
*
TEMP
DPHI
=
TEMP
*
TEMP
ERRETM
=
EIGHT
*
(
-
PHI
-
PSI
)
+
ERRETM
-
PHI
+
RHOINV
+
$
ABS
(
TAU
)
*
(
DPSI
+
DPHI
)
*
W
=
RHOINV
+
PHI
+
PSI
*
*
Main
loop
to
update
the
values
of
the
array
DELTA
*
ITER
=
NITER
+
1
*
DO
90
NITER
=
ITER
,
MAXIT
*
*
Test
for
convergence
*
IF
(
ABS
(
W
)
.LE.
EPS
*
ERRETM
)
THEN
DLAM
=
D
(
I
)
+
TAU
GO
TO
250
END IF
*
IF
(
W
.LE.
ZERO
)
THEN
DLTLB
=
MAX
(
DLTLB
,
TAU
)
ELSE
DLTUB
=
MIN
(
DLTUB
,
TAU
)
END IF
*
*
Calculate
the
new
step
*
C
=
W
-
DELTA
(
N
-
1
)
*
DPSI
-
DELTA
(
N
)
*
DPHI
A
=
(
DELTA
(
N
-
1
)
+
DELTA
(
N
)
)
*
W
-
$
DELTA
(
N
-
1
)
*
DELTA
(
N
)
*
(
DPSI
+
DPHI
)
B
=
DELTA
(
N
-
1
)
*
DELTA
(
N
)
*
W
IF
(
A
.GE.
ZERO
)
THEN
ETA
=
(
A
+
SQRT
(
ABS
(
A
*
A
-
FOUR
*
B
*
C
)
)
)
/
(
TWO
*
C
)
ELSE
ETA
=
TWO
*
B
/
(
A
-
SQRT
(
ABS
(
A
*
A
-
FOUR
*
B
*
C
)
)
)
END IF
*
*
Note
,
eta
should
be
positive
if
w
is
negative
,
and
*
eta
should
be
negative
otherwise
.
However
,
*
if
for
some
reason
caused
by
roundoff
,
eta
*
w
>
0
,
*
we
simply
use
one
Newton
step
instead
.
This
way
*
will
guarantee
eta
*
w
<
0.
*
IF
(
W
*
ETA
.GT.
ZERO
)
$
ETA
=
-
W
/
(
DPSI
+
DPHI
)
TEMP
=
TAU
+
ETA
IF
(
TEMP
.GT.
DLTUB
.OR.
TEMP
.LT.
DLTLB
)
THEN
IF
(
W
.LT.
ZERO
)
THEN
ETA
=
(
DLTUB
-
TAU
)
/
TWO
ELSE
ETA
=
(
DLTLB
-
TAU
)
/
TWO
END IF
END IF
DO
70
J
=
1
,
N
DELTA
(
J
)
=
DELTA
(
J
)
-
ETA
70
CONTINUE
*
TAU
=
TAU
+
ETA
*
*
Evaluate
PSI
and
the
derivative
DPSI
*
DPSI
=
ZERO
PSI
=
ZERO
ERRETM
=
ZERO
DO
80
J
=
1
,
II
TEMP
=
Z
(
J
)
/
DELTA
(
J
)
PSI
=
PSI
+
Z
(
J
)
*
TEMP
DPSI
=
DPSI
+
TEMP
*
TEMP
ERRETM
=
ERRETM
+
PSI
80
CONTINUE
ERRETM
=
ABS
(
ERRETM
)
*
*
Evaluate
PHI
and
the
derivative
DPHI
*
TEMP
=
Z
(
N
)
/
DELTA
(
N
)
PHI
=
Z
(
N
)
*
TEMP
DPHI
=
TEMP
*
TEMP
ERRETM
=
EIGHT
*
(
-
PHI
-
PSI
)
+
ERRETM
-
PHI
+
RHOINV
+
$
ABS
(
TAU
)
*
(
DPSI
+
DPHI
)
*
W
=
RHOINV
+
PHI
+
PSI
90
CONTINUE
*
*
Return
with
INFO
=
1
,
NITER
=
MAXIT
and not
converged
*
INFO
=
1
DLAM
=
D
(
I
)
+
TAU
GO
TO
250
*
*
End
for
the
case
I
=
N
*
ELSE
*
*
The
case
for
I
<
N
*
NITER
=
1
IP1
=
I
+
1
*
*
Calculate
initial
guess
*
DEL
=
D
(
IP1
)
-
D
(
I
)
MIDPT
=
DEL
/
TWO
DO
100
J
=
1
,
N
DELTA
(
J
)
=
(
D
(
J
)
-
D
(
I
)
)
-
MIDPT
100
CONTINUE
*
PSI
=
ZERO
DO
110
J
=
1
,
I
-
1
PSI
=
PSI
+
Z
(
J
)
*
Z
(
J
)
/
DELTA
(
J
)
110
CONTINUE
*
PHI
=
ZERO
DO
120
J
=
N
,
I
+
2
,
-
1
PHI
=
PHI
+
Z
(
J
)
*
Z
(
J
)
/
DELTA
(
J
)
120
CONTINUE
C
=
RHOINV
+
PSI
+
PHI
W
=
C
+
Z
(
I
)
*
Z
(
I
)
/
DELTA
(
I
)
+
$
Z
(
IP1
)
*
Z
(
IP1
)
/
DELTA
(
IP1
)
*
IF
(
W
.GT.
ZERO
)
THEN
*
*
d
(
i
)
<
the
ith
eigenvalue
<
(
d
(
i
)
+
d
(
i
+
1
))
/
2
*
*
We
choose
d
(
i
)
as
origin
.
*
ORGATI
=
.TRUE.
A
=
C
*
DEL
+
Z
(
I
)
*
Z
(
I
)
+
Z
(
IP1
)
*
Z
(
IP1
)
B
=
Z
(
I
)
*
Z
(
I
)
*
DEL
IF
(
A
.GT.
ZERO
)
THEN
TAU
=
TWO
*
B
/
(
A
+
SQRT
(
ABS
(
A
*
A
-
FOUR
*
B
*
C
)
)
)
ELSE
TAU
=
(
A
-
SQRT
(
ABS
(
A
*
A
-
FOUR
*
B
*
C
)
)
)
/
(
TWO
*
C
)
END IF
DLTLB
=
ZERO
DLTUB
=
MIDPT
ELSE
*
*
(
d
(
i
)
+
d
(
i
+
1
))
/
2
<=
the
ith
eigenvalue
<
d
(
i
+
1
)
*
*
We
choose
d
(
i
+
1
)
as
origin
.
*
ORGATI
=
.FALSE.
A
=
C
*
DEL
-
Z
(
I
)
*
Z
(
I
)
-
Z
(
IP1
)
*
Z
(
IP1
)
B
=
Z
(
IP1
)
*
Z
(
IP1
)
*
DEL
IF
(
A
.LT.
ZERO
)
THEN
TAU
=
TWO
*
B
/
(
A
-
SQRT
(
ABS
(
A
*
A
+
FOUR
*
B
*
C
)
)
)
ELSE
TAU
=
-
(
A
+
SQRT
(
ABS
(
A
*
A
+
FOUR
*
B
*
C
)
)
)
/
(
TWO
*
C
)
END IF
DLTLB
=
-
MIDPT
DLTUB
=
ZERO
END IF
*
IF
(
ORGATI
)
THEN
DO
130
J
=
1
,
N
DELTA
(
J
)
=
(
D
(
J
)
-
D
(
I
)
)
-
TAU
130
CONTINUE
ELSE
DO
140
J
=
1
,
N
DELTA
(
J
)
=
(
D
(
J
)
-
D
(
IP1
)
)
-
TAU
140
CONTINUE
END IF
IF
(
ORGATI
)
THEN
II
=
I
ELSE
II
=
I
+
1
END IF
IIM1
=
II
-
1
IIP1
=
II
+
1
*
*
Evaluate
PSI
and
the
derivative
DPSI
*
DPSI
=
ZERO
PSI
=
ZERO
ERRETM
=
ZERO
DO
150
J
=
1
,
IIM1
TEMP
=
Z
(
J
)
/
DELTA
(
J
)
PSI
=
PSI
+
Z
(
J
)
*
TEMP
DPSI
=
DPSI
+
TEMP
*
TEMP
ERRETM
=
ERRETM
+
PSI
150
CONTINUE
ERRETM
=
ABS
(
ERRETM
)
*
*
Evaluate
PHI
and
the
derivative
DPHI
*
DPHI
=
ZERO
PHI
=
ZERO
DO
160
J
=
N
,
IIP1
,
-
1
TEMP
=
Z
(
J
)
/
DELTA
(
J
)
PHI
=
PHI
+
Z
(
J
)
*
TEMP
DPHI
=
DPHI
+
TEMP
*
TEMP
ERRETM
=
ERRETM
+
PHI
160
CONTINUE
*
W
=
RHOINV
+
PHI
+
PSI
*
*
W
is
the
value
of
the
secular
function
with
*
its
ii
-
th
element
removed
.
*
SWTCH3
=
.FALSE.
IF
(
ORGATI
)
THEN
IF
(
W
.LT.
ZERO
)
$
SWTCH3
=
.TRUE.
ELSE
IF
(
W
.GT.
ZERO
)
$
SWTCH3
=
.TRUE.
END IF
IF
(
II
.EQ.
1
.OR.
II
.EQ.
N
)
$
SWTCH3
=
.FALSE.
*
TEMP
=
Z
(
II
)
/
DELTA
(
II
)
DW
=
DPSI
+
DPHI
+
TEMP
*
TEMP
TEMP
=
Z
(
II
)
*
TEMP
W
=
W
+
TEMP
ERRETM
=
EIGHT
*
(
PHI
-
PSI
)
+
ERRETM
+
TWO
*
RHOINV
+
$
THREE
*
ABS
(
TEMP
)
+
ABS
(
TAU
)
*
DW
*
*
Test
for
convergence
*
IF
(
ABS
(
W
)
.LE.
EPS
*
ERRETM
)
THEN
IF
(
ORGATI
)
THEN
DLAM
=
D
(
I
)
+
TAU
ELSE
DLAM
=
D
(
IP1
)
+
TAU
END IF
GO
TO
250
END IF
*
IF
(
W
.LE.
ZERO
)
THEN
DLTLB
=
MAX
(
DLTLB
,
TAU
)
ELSE
DLTUB
=
MIN
(
DLTUB
,
TAU
)
END IF
*
*
Calculate
the
new
step
*
NITER
=
NITER
+
1
IF
(
.NOT.
SWTCH3
)
THEN
IF
(
ORGATI
)
THEN
C
=
W
-
DELTA
(
IP1
)
*
DW
-
(
D
(
I
)
-
D
(
IP1
)
)
*
$
(
Z
(
I
)
/
DELTA
(
I
)
)
**
2
ELSE
C
=
W
-
DELTA
(
I
)
*
DW
-
(
D
(
IP1
)
-
D
(
I
)
)
*
$
(
Z
(
IP1
)
/
DELTA
(
IP1
)
)
**
2
END IF
A
=
(
DELTA
(
I
)
+
DELTA
(
IP1
)
)
*
W
-
$
DELTA
(
I
)
*
DELTA
(
IP1
)
*
DW
B
=
DELTA
(
I
)
*
DELTA
(
IP1
)
*
W
IF
(
C
.EQ.
ZERO
)
THEN
IF
(
A
.EQ.
ZERO
)
THEN
IF
(
ORGATI
)
THEN
A
=
Z
(
I
)
*
Z
(
I
)
+
DELTA
(
IP1
)
*
DELTA
(
IP1
)
*
$
(
DPSI
+
DPHI
)
ELSE
A
=
Z
(
IP1
)
*
Z
(
IP1
)
+
DELTA
(
I
)
*
DELTA
(
I
)
*
$
(
DPSI
+
DPHI
)
END IF
END IF
ETA
=
B
/
A
ELSE IF
(
A
.LE.
ZERO
)
THEN
ETA
=
(
A
-
SQRT
(
ABS
(
A
*
A
-
FOUR
*
B
*
C
)
)
)
/
(
TWO
*
C
)
ELSE
ETA
=
TWO
*
B
/
(
A
+
SQRT
(
ABS
(
A
*
A
-
FOUR
*
B
*
C
)
)
)
END IF
ELSE
*
*
Interpolation
using
THREE
most
relevant
poles
*
TEMP
=
RHOINV
+
PSI
+
PHI
IF
(
ORGATI
)
THEN
TEMP1
=
Z
(
IIM1
)
/
DELTA
(
IIM1
)
TEMP1
=
TEMP1
*
TEMP1
C
=
TEMP
-
DELTA
(
IIP1
)
*
(
DPSI
+
DPHI
)
-
$
(
D
(
IIM1
)
-
D
(
IIP1
)
)
*
TEMP1
ZZ
(
1
)
=
Z
(
IIM1
)
*
Z
(
IIM1
)
ZZ
(
3
)
=
DELTA
(
IIP1
)
*
DELTA
(
IIP1
)
*
$
(
(
DPSI
-
TEMP1
)
+
DPHI
)
ELSE
TEMP1
=
Z
(
IIP1
)
/
DELTA
(
IIP1
)
TEMP1
=
TEMP1
*
TEMP1
C
=
TEMP
-
DELTA
(
IIM1
)
*
(
DPSI
+
DPHI
)
-
$
(
D
(
IIP1
)
-
D
(
IIM1
)
)
*
TEMP1
ZZ
(
1
)
=
DELTA
(
IIM1
)
*
DELTA
(
IIM1
)
*
$
(
DPSI
+
(
DPHI
-
TEMP1
)
)
ZZ
(
3
)
=
Z
(
IIP1
)
*
Z
(
IIP1
)
END IF
ZZ
(
2
)
=
Z
(
II
)
*
Z
(
II
)
CALL
DLAED6
(
NITER
,
ORGATI
,
C
,
DELTA
(
IIM1
),
ZZ
,
W
,
ETA
,
$
INFO
)
IF
(
INFO
.NE.
0
)
$
GO
TO
250
END IF
*
*
Note
,
eta
should
be
positive
if
w
is
negative
,
and
*
eta
should
be
negative
otherwise
.
However
,
*
if
for
some
reason
caused
by
roundoff
,
eta
*
w
>
0
,
*
we
simply
use
one
Newton
step
instead
.
This
way
*
will
guarantee
eta
*
w
<
0.
*
IF
(
W
*
ETA
.GE.
ZERO
)
$
ETA
=
-
W
/
DW
TEMP
=
TAU
+
ETA
IF
(
TEMP
.GT.
DLTUB
.OR.
TEMP
.LT.
DLTLB
)
THEN
IF
(
W
.LT.
ZERO
)
THEN
ETA
=
(
DLTUB
-
TAU
)
/
TWO
ELSE
ETA
=
(
DLTLB
-
TAU
)
/
TWO
END IF
END IF
*
PREW
=
W
*
DO
180
J
=
1
,
N
DELTA
(
J
)
=
DELTA
(
J
)
-
ETA
180
CONTINUE
*
*
Evaluate
PSI
and
the
derivative
DPSI
*
DPSI
=
ZERO
PSI
=
ZERO
ERRETM
=
ZERO
DO
190
J
=
1
,
IIM1
TEMP
=
Z
(
J
)
/
DELTA
(
J
)
PSI
=
PSI
+
Z
(
J
)
*
TEMP
DPSI
=
DPSI
+
TEMP
*
TEMP
ERRETM
=
ERRETM
+
PSI
190
CONTINUE
ERRETM
=
ABS
(
ERRETM
)
*
*
Evaluate
PHI
and
the
derivative
DPHI
*
DPHI
=
ZERO
PHI
=
ZERO
DO
200
J
=
N
,
IIP1
,
-
1
TEMP
=
Z
(
J
)
/
DELTA
(
J
)
PHI
=
PHI
+
Z
(
J
)
*
TEMP
DPHI
=
DPHI
+
TEMP
*
TEMP
ERRETM
=
ERRETM
+
PHI
200
CONTINUE
*
TEMP
=
Z
(
II
)
/
DELTA
(
II
)
DW
=
DPSI
+
DPHI
+
TEMP
*
TEMP
TEMP
=
Z
(
II
)
*
TEMP
W
=
RHOINV
+
PHI
+
PSI
+
TEMP
ERRETM
=
EIGHT
*
(
PHI
-
PSI
)
+
ERRETM
+
TWO
*
RHOINV
+
$
THREE
*
ABS
(
TEMP
)
+
ABS
(
TAU
+
ETA
)
*
DW
*
SWTCH
=
.FALSE.
IF
(
ORGATI
)
THEN
IF
(
-
W
.GT.
ABS
(
PREW
)
/
TEN
)
$
SWTCH
=
.TRUE.
ELSE
IF
(
W
.GT.
ABS
(
PREW
)
/
TEN
)
$
SWTCH
=
.TRUE.
END IF
*
TAU
=
TAU
+
ETA
*
*
Main
loop
to
update
the
values
of
the
array
DELTA
*
ITER
=
NITER
+
1
*
DO
240
NITER
=
ITER
,
MAXIT
*
*
Test
for
convergence
*
IF
(
ABS
(
W
)
.LE.
EPS
*
ERRETM
)
THEN
IF
(
ORGATI
)
THEN
DLAM
=
D
(
I
)
+
TAU
ELSE
DLAM
=
D
(
IP1
)
+
TAU
END IF
GO
TO
250
END IF
*
IF
(
W
.LE.
ZERO
)
THEN
DLTLB
=
MAX
(
DLTLB
,
TAU
)
ELSE
DLTUB
=
MIN
(
DLTUB
,
TAU
)
END IF
*
*
Calculate
the
new
step
*
IF
(
.NOT.
SWTCH3
)
THEN
IF
(
.NOT.
SWTCH
)
THEN
IF
(
ORGATI
)
THEN
C
=
W
-
DELTA
(
IP1
)
*
DW
-
$
(
D
(
I
)
-
D
(
IP1
)
)
*
(
Z
(
I
)
/
DELTA
(
I
)
)
**
2
ELSE
C
=
W
-
DELTA
(
I
)
*
DW
-
(
D
(
IP1
)
-
D
(
I
)
)
*
$
(
Z
(
IP1
)
/
DELTA
(
IP1
)
)
**
2
END IF
ELSE
TEMP
=
Z
(
II
)
/
DELTA
(
II
)
IF
(
ORGATI
)
THEN
DPSI
=
DPSI
+
TEMP
*
TEMP
ELSE
DPHI
=
DPHI
+
TEMP
*
TEMP
END IF
C
=
W
-
DELTA
(
I
)
*
DPSI
-
DELTA
(
IP1
)
*
DPHI
END IF
A
=
(
DELTA
(
I
)
+
DELTA
(
IP1
)
)
*
W
-
$
DELTA
(
I
)
*
DELTA
(
IP1
)
*
DW
B
=
DELTA
(
I
)
*
DELTA
(
IP1
)
*
W
IF
(
C
.EQ.
ZERO
)
THEN
IF
(
A
.EQ.
ZERO
)
THEN
IF
(
.NOT.
SWTCH
)
THEN
IF
(
ORGATI
)
THEN
A
=
Z
(
I
)
*
Z
(
I
)
+
DELTA
(
IP1
)
*
$
DELTA
(
IP1
)
*
(
DPSI
+
DPHI
)
ELSE
A
=
Z
(
IP1
)
*
Z
(
IP1
)
+
$
DELTA
(
I
)
*
DELTA
(
I
)
*
(
DPSI
+
DPHI
)
END IF
ELSE
A
=
DELTA
(
I
)
*
DELTA
(
I
)
*
DPSI
+
$
DELTA
(
IP1
)
*
DELTA
(
IP1
)
*
DPHI
END IF
END IF
ETA
=
B
/
A
ELSE IF
(
A
.LE.
ZERO
)
THEN
ETA
=
(
A
-
SQRT
(
ABS
(
A
*
A
-
FOUR
*
B
*
C
)
)
)
/
(
TWO
*
C
)
ELSE
ETA
=
TWO
*
B
/
(
A
+
SQRT
(
ABS
(
A
*
A
-
FOUR
*
B
*
C
)
)
)
END IF
ELSE
*
*
Interpolation
using
THREE
most
relevant
poles
*
TEMP
=
RHOINV
+
PSI
+
PHI
IF
(
SWTCH
)
THEN
C
=
TEMP
-
DELTA
(
IIM1
)
*
DPSI
-
DELTA
(
IIP1
)
*
DPHI
ZZ
(
1
)
=
DELTA
(
IIM1
)
*
DELTA
(
IIM1
)
*
DPSI
ZZ
(
3
)
=
DELTA
(
IIP1
)
*
DELTA
(
IIP1
)
*
DPHI
ELSE
IF
(
ORGATI
)
THEN
TEMP1
=
Z
(
IIM1
)
/
DELTA
(
IIM1
)
TEMP1
=
TEMP1
*
TEMP1
C
=
TEMP
-
DELTA
(
IIP1
)
*
(
DPSI
+
DPHI
)
-
$
(
D
(
IIM1
)
-
D
(
IIP1
)
)
*
TEMP1
ZZ
(
1
)
=
Z
(
IIM1
)
*
Z
(
IIM1
)
ZZ
(
3
)
=
DELTA
(
IIP1
)
*
DELTA
(
IIP1
)
*
$
(
(
DPSI
-
TEMP1
)
+
DPHI
)
ELSE
TEMP1
=
Z
(
IIP1
)
/
DELTA
(
IIP1
)
TEMP1
=
TEMP1
*
TEMP1
C
=
TEMP
-
DELTA
(
IIM1
)
*
(
DPSI
+
DPHI
)
-
$
(
D
(
IIP1
)
-
D
(
IIM1
)
)
*
TEMP1
ZZ
(
1
)
=
DELTA
(
IIM1
)
*
DELTA
(
IIM1
)
*
$
(
DPSI
+
(
DPHI
-
TEMP1
)
)
ZZ
(
3
)
=
Z
(
IIP1
)
*
Z
(
IIP1
)
END IF
END IF
CALL
DLAED6
(
NITER
,
ORGATI
,
C
,
DELTA
(
IIM1
),
ZZ
,
W
,
ETA
,
$
INFO
)
IF
(
INFO
.NE.
0
)
$
GO
TO
250
END IF
*
*
Note
,
eta
should
be
positive
if
w
is
negative
,
and
*
eta
should
be
negative
otherwise
.
However
,
*
if
for
some
reason
caused
by
roundoff
,
eta
*
w
>
0
,
*
we
simply
use
one
Newton
step
instead
.
This
way
*
will
guarantee
eta
*
w
<
0.
*
IF
(
W
*
ETA
.GE.
ZERO
)
$
ETA
=
-
W
/
DW
TEMP
=
TAU
+
ETA
IF
(
TEMP
.GT.
DLTUB
.OR.
TEMP
.LT.
DLTLB
)
THEN
IF
(
W
.LT.
ZERO
)
THEN
ETA
=
(
DLTUB
-
TAU
)
/
TWO
ELSE
ETA
=
(
DLTLB
-
TAU
)
/
TWO
END IF
END IF
*
DO
210
J
=
1
,
N
DELTA
(
J
)
=
DELTA
(
J
)
-
ETA
210
CONTINUE
*
TAU
=
TAU
+
ETA
PREW
=
W
*
*
Evaluate
PSI
and
the
derivative
DPSI
*
DPSI
=
ZERO
PSI
=
ZERO
ERRETM
=
ZERO
DO
220
J
=
1
,
IIM1
TEMP
=
Z
(
J
)
/
DELTA
(
J
)
PSI
=
PSI
+
Z
(
J
)
*
TEMP
DPSI
=
DPSI
+
TEMP
*
TEMP
ERRETM
=
ERRETM
+
PSI
220
CONTINUE
ERRETM
=
ABS
(
ERRETM
)
*
*
Evaluate
PHI
and
the
derivative
DPHI
*
DPHI
=
ZERO
PHI
=
ZERO
DO
230
J
=
N
,
IIP1
,
-
1
TEMP
=
Z
(
J
)
/
DELTA
(
J
)
PHI
=
PHI
+
Z
(
J
)
*
TEMP
DPHI
=
DPHI
+
TEMP
*
TEMP
ERRETM
=
ERRETM
+
PHI
230
CONTINUE
*
TEMP
=
Z
(
II
)
/
DELTA
(
II
)
DW
=
DPSI
+
DPHI
+
TEMP
*
TEMP
TEMP
=
Z
(
II
)
*
TEMP
W
=
RHOINV
+
PHI
+
PSI
+
TEMP
ERRETM
=
EIGHT
*
(
PHI
-
PSI
)
+
ERRETM
+
TWO
*
RHOINV
+
$
THREE
*
ABS
(
TEMP
)
+
ABS
(
TAU
)
*
DW
IF
(
W
*
PREW
.GT.
ZERO
.AND.
ABS
(
W
)
.GT.
ABS
(
PREW
)
/
TEN
)
$
SWTCH
=
.NOT.
SWTCH
*
240
CONTINUE
*
*
Return
with
INFO
=
1
,
NITER
=
MAXIT
and not
converged
*
INFO
=
1
IF
(
ORGATI
)
THEN
DLAM
=
D
(
I
)
+
TAU
ELSE
DLAM
=
D
(
IP1
)
+
TAU
END IF
*
END IF
*
250
CONTINUE
*
RETURN
*
*
End
of
DLAED4
*
END
Event Timeline
Log In to Comment