From eeed18aeffc7175b16c80ee1f5f402d91ab10393 Mon Sep 17 00:00:00 2001
From: farhadrclass <31899325+farhadrclass@users.noreply.github.com>
Date: Thu, 21 Aug 2025 17:06:13 -0400
Subject: [PATCH 1/2] Add robust opnorm implementation and tests

Introduces a new opnorm function that efficiently computes the operator 2-norm for matrices and linear operators, dispatching to LAPACK, ARPACK, or TSVD as appropriate. Adds comprehensive tests for opnorm covering both Float64 and BigFloat types, and integrates the new test file into the test suite.

Update Project.toml
---
 Project.toml        |   4 ++
 src/utilities.jl    | 122 ++++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl    |   1 +
 test/test_opnorm.jl |  31 +++++++++++
 4 files changed, 158 insertions(+)
 create mode 100644 test/test_opnorm.jl

diff --git a/Project.toml b/Project.toml
index 6f31907..0bdd1f2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,10 +4,12 @@ version = "2.11.0"
 
 [deps]
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [weakdeps]
@@ -31,6 +33,7 @@ AMDGPU = "2"
 CUDA = "5"
 ChainRulesCore = "1"
 FastClosures = "0.2, 0.3"
+GenericLinearAlgebra = "0.3.18"
 JLArrays = "0.1, 0.2"
 LDLFactorizations = "0.9, 0.10"
 LinearAlgebra = "1"
@@ -38,6 +41,7 @@ Metal = "1.1"
 Printf = "1"
 Requires = "1"
 SparseArrays = "1"
+TSVD = "0.4.4"
 TimerOutputs = "^0.5"
 julia = "^1.6.0"
 
diff --git a/src/utilities.jl b/src/utilities.jl
index 1252d95..111a520 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -1,6 +1,11 @@
 export check_ctranspose,
   check_hermitian, check_positive_definite, normest, solve_shifted_system!, ldiv!
 import LinearAlgebra.ldiv!
+import LinearAlgebra: opnorm       
+using GenericLinearAlgebra
+using TSVD
+
+export opnorm
 
 """
   normest(S) estimates the matrix 2-norm of S.
@@ -287,3 +292,120 @@ function ldiv!(
   solve_shifted_system!(x, B, b, T(0.0))
   return x
 end
+
+"""
+  opnorm(B; kwargs...)
+
+Compute the operator 2-norm (largest singular value) of a matrix or linear operator `B`.
+This method dispatches to efficient algorithms depending on the type and size of `B`:
+for small dense matrices, it uses direct LAPACK routines; for larger or structured operators,
+it uses iterative methods (ARPACK or TSVD) to estimate the norm efficiently.
+
+# Arguments
+- `B`: A matrix or linear operator.
+- `kwargs...`: Optional keyword arguments passed to the underlying norm estimation routines.
+
+# Returns
+- The estimated operator 2-norm of `B` (largest singular value or eigenvalue in absolute value).
+"""
+
+function LinearAlgebra.opnorm(B; kwargs...)
+  _opnorm(B, eltype(B); kwargs...)
+end
+
+# This method will be picked if eltype is one of the four types Arpack supports
+# (Float32, Float64, ComplexF32, ComplexF64).
+function _opnorm(
+  B,
+  ::Type{T};
+  kwargs...,
+) where {T <: Union{Float32, Float64, ComplexF32, ComplexF64}}
+  m, n = size(B)
+  return (m == n ? opnorm_eig : opnorm_svd)(B; kwargs...)
+end
+
+# Fallback for everything else
+function _opnorm(B, ::Type{T}; kwargs...) where {T}
+  _, s, _ = tsvd(B)
+  return s[1], true  # return largest singular value
+end
+
+function opnorm_eig(B; max_attempts::Int = 3)
+  n = size(B, 1)
+  # 1) tiny dense Float64: direct LAPACK
+  if n ≤ 5
+    return maximum(abs, eigen(Matrix(B)).values), true
+  end
+
+  # 2) iterative ARPACK
+  nev, ncv = 1, max(20, 2*1 + 1)
+  attempt, λ, have_eig = 0, zero(eltype(B)), false
+
+  while !(have_eig || attempt >= max_attempts)
+    attempt += 1
+    try
+      # Estimate largest eigenvalue in absolute value
+      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)) && ncv < n
+        @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)
+  m, n = size(J)
+  # 1) tiny dense Float64: direct LAPACK
+  if min(m, n) ≤ 5
+    return maximum(svd(Matrix(J)).S), true
+  end
+
+  # 2) iterative ARPACK‐SVD
+  nsv, ncv = 1, 10
+  attempt, σ, have_svd = 0, zero(eltype(J)), false
+  n = min(m, n)
+
+  while !(have_svd || attempt >= max_attempts)
+    attempt += 1
+    try
+      # Estimate largest singular value
+      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)) && ncv < n
+        @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
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 2c8e6d1..196aecb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -13,6 +13,7 @@ include("test_kron.jl")
 include("test_callable.jl")
 include("test_deprecated.jl")
 include("test_normest.jl")
+include("test_opnorm.jl")
 include("test_diag.jl")
 include("test_chainrules.jl")
 include("test_solve_shifted_system.jl")
diff --git a/test/test_opnorm.jl b/test/test_opnorm.jl
new file mode 100644
index 0000000..f610f04
--- /dev/null
+++ b/test/test_opnorm.jl
@@ -0,0 +1,31 @@
+function test_opnorm()
+    @testset "opnorm and TSVD tests (LinearOperators)" begin
+        # 1) Square Float64 via direct LAPACK or ARPACK
+        A_mat = [2.0 0.0; 0.0 -1.0]
+        A_op = LinearOperator(A_mat)
+        λ, ok = opnorm(A_op)
+        @test ok == true
+        @test isapprox(λ, 2.0; rtol=1e-12)
+
+        # 2) Rectangular Float64 via direct LAPACK or ARPACK SVD
+        J_mat = [3.0 0.0 0.0; 0.0 1.0 0.0]
+        J_op = LinearOperator(J_mat)
+        σ, ok_sv = opnorm(J_op)
+        @test ok_sv == true
+        @test isapprox(σ, 3.0; rtol=1e-12)
+
+        # 3) Square BigFloat via TSVD
+        B_mat = Matrix{BigFloat}([2.0 0.0; 0.0 -1.0])
+        B_op = LinearOperator(B_mat)
+        λ_bf, ok_bf = opnorm(B_op)
+        @test ok_bf == true
+        @test isapprox(λ_bf, BigFloat(2); rtol=1e-12)
+
+        # 4) Rectangular BigFloat via rectangular TSVD
+        JR_mat = Matrix{BigFloat}([3.0 0.0 0.0; 0.0 1.0 0.0])
+        JR_op = LinearOperator(JR_mat)
+        σ_bf, ok_bf2 = opnorm(JR_op)
+        @test ok_bf2 == true
+        @test isapprox(σ_bf, BigFloat(3); rtol=1e-12)
+    end
+end 
\ No newline at end of file

From a67ac201806101a4676524fe5d3e90ff283e00da Mon Sep 17 00:00:00 2001
From: farhadrclass <31899325+farhadrclass@users.noreply.github.com>
Date: Thu, 21 Aug 2025 18:34:53 -0400
Subject: [PATCH 2/2] Update test_opnorm.jl

---
 test/test_opnorm.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/test/test_opnorm.jl b/test/test_opnorm.jl
index f610f04..c178ab6 100644
--- a/test/test_opnorm.jl
+++ b/test/test_opnorm.jl
@@ -28,4 +28,6 @@ function test_opnorm()
         @test ok_bf2 == true
         @test isapprox(σ_bf, BigFloat(3); rtol=1e-12)
     end
-end 
\ No newline at end of file
+end
+
+test_opnorm()
\ No newline at end of file