Use shifted Arnoldi method for finding just 1 eval/vec closest to guess

Author: duetosymmetry

Kernel/SpinWeightedSpheroidalHarmonics.m

@@ -99,6 +99,10 @@
   ,{nUp+nDown+1,nUp+nDown+1}];
 
 
+(* Mma ought to provide a sparse identity matrix but... *)
+SparseIdentityMatrix[n_Integer?NonNegative]:=SparseArray[{{x_,x_}->1},{n,n}];
+
+
 (* ::Subsection::Closed:: *)
 (*Continued fraction*)
 
@@ -132,12 +136,22 @@
 (*Spherical expansion method*)
 
 
-Options[SWSHEigenvalueSpectral] = {"NumTerms" -> Automatic};
+Options[SWSHEigenvalueSpectral] = {"InitialGuess" -> Automatic, "NumTerms" -> Automatic};
 
 
-SWSHEigenvalueSpectral[s_, l_, m_, \[Gamma]_, OptionsPattern[]]:=
- Module[{nDown, nUp, Matrix, Eigens, lmin},
-  (* FIXME: Improve the estimate of nmax. It should depend on the accuarcy sought. *)
+SWSHEigenvalueSpectral[s_, l_, m_, \[Gamma]_, opts:OptionsPattern[]]:=
+ Module[{nDown, nUp, AMatrix, lmin, AGuess},
+  Switch[OptionValue["InitialGuess"],
+   Automatic,
+     (* NOTE this is currently an initial value of A, not \[Lambda] *)
+     AGuess = l(l+1) - s(s+1);,
+   _?NumericQ,
+     AGuess = OptionValue["InitialGuess"];,
+   _,
+    Message[SpinWeightedSpheroidalEigenvalue::optx, Method -> {"SphericalExpansion", opts}];
+    Return[$Failed];
+  ];
+  (* FIXME: Improve the estimate of nmax. It should depend on the accuracy sought. *)
   Switch[OptionValue["NumTerms"],
    Automatic,
     nUp = Ceiling[Abs[3/2\[Gamma]]]+5;
@@ -151,13 +165,13 @@
   lmin=Max[Abs[s],Abs[m]];
   nDown = Min[nUp, l-lmin];
 
-  Matrix=SparseSWSHMatrix[s,l,m,\[Gamma],nDown,nUp];
-
-  (* To choose the eigenvalue corrsponding to the desired l, we assume that the real
-     part of eigenvalue is a monotonic function of l. *)
-  Eigens=-Sort[Quiet[Eigenvalues[Matrix],Eigenvalues::arhm]];
-
-  Eigens[[-(nDown+1)]]-s(s+1)
+  AMatrix = SparseSWSHMatrix[s,l,m,\[Gamma],nDown,nUp];
+  (* Because the matrix built up above has an annoying convention... *)
+  (* After this step, the eigenvalues are A's *)
+  AMatrix = -AMatrix - s(s+1)*SparseIdentityMatrix[nUp+nDown+1];
+  (* The above manipulation can be absorbed into the shift below, and flipping the sign of the eigenvalue *)
+  (* Do we want to do that? *)
+  Eigenvalues[AMatrix, -1, Method->{"Arnoldi", "Shift"->AGuess}][[1]]
 ];
 
 
@@ -326,11 +340,21 @@
 (*Spherical expansion method*)
 
 
-Options[SWSHSSpectral] = {"NumTerms" -> Automatic};
+Options[SWSHSSpectral] = {"InitialGuess"->Automatic, "NumTerms" -> Automatic};
 
 
 SWSHSSpectral[s_Integer, l_Integer, m_Integer, \[Gamma]_, opts:OptionsPattern[]] :=
- Module[{lmin, nUp, nDown, A, esys,evec,eval,sign,pos},
+ Module[{lmin, nUp, nDown, AMatrix, evec,sign,AGuess},
+   Switch[OptionValue["InitialGuess"],
+   Automatic,
+     (* NOTE this is currently an initial value of A, not \[Lambda] *)
+     AGuess = l(l+1) - s(s+1);,
+   _?NumericQ,
+     AGuess = OptionValue["InitialGuess"];,
+   _,
+    Message[SpinWeightedSpheroidalEigenvalue::optx, Method -> {"SphericalExpansion", opts}];
+    Return[$Failed];
+  ];
   (* FIXME: Improve the estimate of nmax. It should depend on the accuracy sought. *)
   Switch[OptionValue["NumTerms"],
    Automatic,
@@ -345,11 +369,13 @@
   lmin = Max[Abs[s],Abs[m]];
   nDown = Min[l-lmin,nUp];
 
-  A = SparseSWSHMatrix[s,l,m,\[Gamma],nDown,nUp];
-  esys = Quiet[Eigensystem[A], Eigenvalues::arhm];
-  eval = -Sort[esys[[1]]][[-(nDown+1)]];
-  pos  = Position[esys[[1]], -eval][[1]];
-  evec = First[esys[[2,pos]]];
+  AMatrix = SparseSWSHMatrix[s,l,m,\[Gamma],nDown,nUp];
+  (* Because the matrix built up above has an annoying convention... *)
+  (* After this step, the eigenvalues are A's *)
+  AMatrix = -AMatrix - s(s+1)*SparseIdentityMatrix[nUp+nDown+1];
+  (* The above manipulation can be absorbed into the shift below, and flipping the sign of the eigenvalue *)
+  (* Do we want to do that? *)
+  evec = Eigenvectors[AMatrix, -1, Method->{"Arnoldi", "Shift"->AGuess}][[1]];
 
   sign=Sign[evec[[Min[l-lmin+1,(nUp+nDown)/2+1]]]];
   SpinWeightedSpheroidalHarmonicSFunction[s, l, m, \[Gamma], {"SphericalExpansion", sign*evec, nDown, nUp}]