Update utils.jl

Enhance robustness of Arpack usage by catching NCV-related errors

Catch the XYAUPD_Exception error in Arpack calls and dynamically increase NCV as needed to handle convergence issues. This modification addresses instability in LSR1 computations due to insufficient NCV.

Author: MohamedLaghdafHABIBOULLAH

src/utils.jl

@@ -11,36 +11,68 @@ function opnorm_eig(B; max_attempts::Int = 3)
   n = size(B, 1)
   nev = 1
   ncv = max(20, 2 * nev + 1)
-  while !(have_eig || attempt > max_attempts)
-    attempt += 1
-    (d, nconv, niter, nmult, resid) = eigs(B; nev = nev, ncv = ncv, which = :LM, ritzvec = false, check = 1)
-    have_eig = nconv == 1
-    if (have_eig)
-      λ = abs(d[1])
-    else
-      ncv = min(2 * ncv, n)
-    end
+
+  while !(have_eig || attempt >= max_attempts)
+      attempt += 1
+      try
+          # Perform eigendecomposition
+          d, nconv, niter, nmult, resid = eigs(B; nev = nev, ncv = ncv, which = :LM, ritzvec = false, check = 1)
+          
+          # Check if eigenvalue has converged
+          have_eig = nconv == 1
+          if have_eig
+              λ = abs(d[1])  # Take absolute value of the largest eigenvalue
+              break  # Exit loop if successful
+          else
+              # Increase NCV for the next attempt if convergence wasn't achieved
+              ncv = min(2 * ncv, n)
+          end
+      catch e
+          if occursin("XYAUPD_Exception", string(e))
+              @warn "Arpack error: $e. Increasing NCV to $ncv and retrying."
+              ncv = min(2 * ncv, n)  # Increase NCV but don't exceed matrix size
+          else
+              rethrow(e)  # Re-raise if it's a different error
+          end
+      end
   end
+
   return λ, have_eig
 end
 
 function opnorm_svd(J; max_attempts::Int = 3)
   have_svd = false
   attempt = 0
   σ = zero(eltype(J))
-  n = min(size(J)...)
+  n = min(size(J)...)  # Minimum dimension of the matrix
   nsv = 1
   ncv = 10
-  while !(have_svd || attempt > max_attempts)
-    attempt += 1
-    (s, nconv, niter, nmult, resid) = svds(J, nsv = nsv, ncv = ncv, ritzvec = false, check = 1)
-    have_svd = nconv == 1
-    if (have_svd)
-      σ = maximum(s.S)
-    else
-      ncv = min(2 * ncv, n)
-    end
+
+  while !(have_svd || attempt >= max_attempts)
+      attempt += 1
+      try
+          # Perform singular value decomposition
+          s, nconv, niter, nmult, resid = svds(J; nsv = nsv, ncv = ncv, ritzvec = false, check = 1)
+          
+          # Check if singular value has converged
+          have_svd = nconv >= 1
+          if have_svd
+              σ = maximum(s.S)  # Take the largest singular value
+              break  # Exit loop if successful
+          else
+              # Increase NCV for the next attempt if convergence wasn't achieved
+              ncv = min(2 * ncv, n)
+          end
+      catch e
+          if occursin("XYAUPD_Exception", string(e))
+              @warn "Arpack error: $e. Increasing NCV to $ncv and retrying."
+              ncv = min(2 * ncv, n)  # Increase NCV but don't exceed matrix size
+          else
+              rethrow(e)  # Re-raise if it's a different error
+          end
+      end
   end
+
   return σ, have_svd
 end