ClearAll[DPR1FindEigenvalue];
Options[DPR1FindEigenvalue] = {
"Tolerance" -> 100 $MachineEpsilon,
"MaxIterations" -> 20
};
DPR1FindEigenvalue[diag_?VectorQ, z_?VectorQ, OptionsPattern[]] :=
Module[{p, zz, a0, b0},
p = Ordering[diag];
zz = (z z)[[p]];
a0 = diag[[p]];
b0 = Append[Rest[a0], a0[[-1]] + Total[zz]];
Reverse[
cDPR1FindEigenvalue[a0, zz, a0, b0, OptionValue["MaxIterations"],
OptionValue["Tolerance"]]]
];
cDPR1FindEigenvalue = With[{
eps = N[$MachineEpsilon]
},
Compile[{{diag, _Real, 1}, {zz, _Real, 1}, {a0, _Real}, {b0, _Real},
{maxiter, _Integer}, {RelTOL, _Real}},
Module[{n, a, fa, b, fb, c, fc, Dfc, TOL, cold, fcold,
x, \[Delta]c, iter, \[Epsilon], sum, \[Delta]},
n = Length[diag];
a = a0;
b = b0;
\[Epsilon] = 10. eps;
(*TOL=RelTOL(b-a);*)
TOL = RelTOL Max[Abs[b], Abs[a]];
c = 0.5 (a + b);
(*Evaluate the characteristic polynomial f at c.*)
sum = 0.;
Do[
\[Delta] = Compile`GetElement[diag, i] - c;
sum += Compile`GetElement[zz, i]/\[Delta];
, {i, 1, n}
];
fc = 1. + sum;
cold = c;
fcold = fc;
(*First we have to bracket the solution.*)
(*Trying to find {a,b} with f[a] < 0 < f[b].*)
If[fc >= 0.,
(
While[fc >= 0.,
cold = c; fcold = fc;
c = 0.5 (a + c);
(*Evaluate the characteristic polynomial f at c.*)
sum = 0.;
Do[
\[Delta] = Compile`GetElement[diag, i] - c;
sum += Compile`GetElement[zz, i]/\[Delta];
, {i, 1, n}
];
fc = 1. + sum;
];
a = c; fa = fc;
b = cold; fb = fcold;
)
,
(
While[fc < 0.,
cold = c; fcold = fc;
c = 0.5 (c + b);
(*Evaluate the characteristic polynomial f at c.*)
sum = 0.;
Do[
\[Delta] = Compile`GetElement[diag, i] - c;
sum += Compile`GetElement[zz, i]/\[Delta];
, {i, 1, n}
];
fc = 1. + sum;
];
a = cold; fa = fcold;
b = c; fb = fc;
)
];
(*Refinement.*)
c = 0.5 (a + b);
(*Evaluate the characteristic polynomial f at c.*)
sum = 0.;
Do[
\[Delta] = Compile`GetElement[diag, i] - c;
sum += Compile`GetElement[zz, i]/\[Delta];
, {i, 1, n}
];
fc = 1. + sum;
iter = 0;
While[iter < maxiter,
++iter;
(*First try a Newton step.*)
(*\[Delta]c=-fc/cDf[diag,zz,c];*)
(*Evaluate the derivative Dfc of the characteristic polynomial \
f at c.*)
Dfc = 0.;
Do[
\[Delta] = Compile`GetElement[diag, i] - c;
Dfc += Compile`GetElement[zz, i]/(\[Delta] \[Delta]);
, {i, 1, n}
];
\[Delta]c = -fc/Dfc;
x = c + \[Delta]c;
If[a < x < b,
(
(*Newton's method gives a valid iterate.*)
If[Abs[\[Delta]c] < Max[\[Epsilon], TOL],
(*Print["A"->a<x<b];*)
Return[c];
];
c = x;
(*Evaluate the characteristic polynomial f at c.*)
sum = 0.;
Do[
\[Delta] = Compile`GetElement[diag, i] - c;
sum += Compile`GetElement[zz, i]/\[Delta];
, {i, 1, n}
];
fc = 1. + sum;
)
,
(
(*Newton's method leaves the bracketing interval;
we continue with the interval's midpoint.*)
If[fc < 0.,
a = c; fa = fc;
,
b = c; fb = fc;
];
If[Abs[b - a] < Max[\[Epsilon], TOL],
Return[c];
];
c = 0.5 (a + b);
(*Evaluate the characteristic polynomial f at c.*)
sum = 0.;
Do[
\[Delta] = Compile`GetElement[diag, i] - c;
sum += Compile`GetElement[zz, i]/\[Delta];
, {i, 1, n}
];
fc = 1. + sum;
)
];
];
(*We ran out of iterations. Return the result anyways.*)
Return[c]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]
];
ClearAll[DPR1Eigensystem];
Options[DPR1Eigensystem] = {
"Tolerance" -> 100 $MachineEpsilon,
"MaxIterations" -> 20
};
DPR1Eigensystem[diag_?VectorQ, z_?VectorQ, OptionsPattern[]] :=
Module[{\[Lambda], U},
\[Lambda] = DPR1Eigenvalues[diag, z,
"MaxIterations" -> OptionValue["MaxIterations"],
"Tolerance" -> OptionValue["Tolerance"]
];
U = DPR1EigenvectorFromEigenvalue[diag, z, \[Lambda]];
{\[Lambda], U}
];
This has grown in a more complex project than I thought. Thus for maintenance reasons, I decided to create the following github repository:
https://github.com/HenrikSchumacher/DPR1Eigensystem
Needs["DPR1Eigensystem`"];
n = 10000;
diag = RandomReal[{1, n}, n];
z = RandomReal[{-1, 1}, n];
{\[Lambda]True, UTrue} = Eigensystem[DiagonalMatrix[diag] + KroneckerProduct[z, z]]; // AbsoluteTiming // First
{\[Lambda], U} = DPR1Eigensystem[diag, z, MaxIterations -> 20, Tolerance -> 10^-14]; // AbsoluteTiming // First
Max[Abs[\[Lambda]True - \[Lambda]]]/Max[Abs[\[Lambda]True]]
Max[Table[ Max[Abs[DPR1Apply[diag, z, U[[k]]] - \[Lambda][[k]] U[[k]]]], {k, 1, n}]]
