add modified chebyshev generic routine

rename variables in STpH

Author: MikaelSlevinsky

Project.toml

@@ -1,6 +1,6 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.16.1"
+version = "0.16.2"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/SymmetricToeplitzPlusHankel.jl

@@ -37,7 +37,7 @@ end
 bandwidths(A::SymmetricBandedToeplitzPlusHankel{T}) where T = (A.b, A.b)
 
 #
-# Jac*W-W*Jac' = G*J*G'
+# X'W-W*X = G*J*G'
 # This returns G and J, where J = [0 I; -I 0], respecting the skew-symmetry of the right-hand side.
 #
 function compute_skew_generators(A::SymmetricToeplitzPlusHankel{T}) where T
@@ -56,80 +56,82 @@ function cholesky(A::SymmetricToeplitzPlusHankel{T}) where T
     n = size(A, 1)
     G, J = compute_skew_generators(A)
     L = zeros(T, n, n)
-    r = A[:, 1]
-    r2 = zeros(T, n)
+    c = A[:, 1]
+    ĉ = zeros(T, n)
     l = zeros(T, n)
     v = zeros(T, n)
-    col1 = zeros(T, n)
-    STPHcholesky!(L, G, r, r2, l, v, col1, n)
+    row1 = zeros(T, n)
+    STPHcholesky!(L, G, c, ĉ, l, v, row1, n)
     return Cholesky(L, 'L', 0)
 end
 
-function STPHcholesky!(L::Matrix{T}, G, r, r2, l, v, col1, n) where T
+function STPHcholesky!(L::Matrix{T}, G, c, ĉ, l, v, row1, n) where T
     @inbounds @simd for k in 1:n-1
-        x = sqrt(r[1])
+        d = sqrt(c[1])
         for j in 1:n-k+1
-            L[j+k-1, k] = l[j] = r[j]/x
+            L[j+k-1, k] = l[j] = c[j]/d
         end
         for j in 1:n-k+1
-            v[j] = G[j, 1]*G[1,2]-G[j,2]*G[1,1]
+            v[j] = G[j, 1]*G[1, 2] - G[j, 2]*G[1, 1]
         end
-        F12 = k == 1 ? T(2) : T(1)
-        r2[1] = (r[2] - v[1])/F12
+        X21 = k == 1 ? T(2) : T(1)
+        ĉ[1] = (c[2] - v[1])/X21
         for j in 2:n-k
-            r2[j] = (r[j+1]+r[j-1] + r[1]*col1[j] - col1[1]*r[j] - v[j])/F12
+            ĉ[j] = (c[j+1] + c[j-1] + c[1]*row1[j] - row1[1]*c[j] - v[j])/X21
         end
-        r2[n-k+1] = (r[n-k] + r[1]*col1[n-k+1] - col1[1]*r[n-k+1] - v[n-k+1])/F12
-        cst = r[2]/x
+        ĉ[n-k+1] = (c[n-k] + c[1]*row1[n-k+1] - row1[1]*c[n-k+1] - v[n-k+1])/X21
+        cst = c[2]/d
         for j in 1:n-k
-            r[j] = r2[j+1] - cst*l[j+1]
+            c[j] = ĉ[j+1] - cst*l[j+1]
         end
+        cst = X21/d
         for j in 1:n-k
-            col1[j] = -F12/x*l[j+1]
+            row1[j] = -cst*l[j+1]
         end
-        c1 = G[1, 1]
-        c2 = G[1, 2]
+        gd1 = G[1, 1]/d
+        gd2 = G[1, 2]/d
         for j in 1:n-k
-            G[j, 1] = G[j+1, 1] - l[j+1]*c1/x
-            G[j, 2] = G[j+1, 2] - l[j+1]*c2/x
+            G[j, 1] = G[j+1, 1] - l[j+1]*gd1
+            G[j, 2] = G[j+1, 2] - l[j+1]*gd2
         end
     end
-    L[n, n] = sqrt(r[1])
+    L[n, n] = sqrt(c[1])
 end
 
 function cholesky(A::SymmetricBandedToeplitzPlusHankel{T}) where T
     n = size(A, 1)
     b = A.b
     R = BandedMatrix{T}(undef, (n, n), (0, bandwidth(A, 2)))
-    r = A[1:b+2, 1]
-    r2 = zeros(T, b+3)
+    c = A[1:b+2, 1]
+    ĉ = zeros(T, b+3)
     l = zeros(T, b+3)
-    col1 = zeros(T, b+2)
-    SBTPHcholesky!(R, r, r2, l, col1, n, b)
+    row1 = zeros(T, b+2)
+    SBTPHcholesky!(R, c, ĉ, l, row1, n, b)
     return Cholesky(R, 'U', 0)
 end
 
-function SBTPHcholesky!(R::BandedMatrix{T}, r, r2, l, col1, n, b) where T
+function SBTPHcholesky!(R::BandedMatrix{T}, c, ĉ, l, row1, n, b) where T
     @inbounds @simd for k in 1:n
-        x = sqrt(r[1])
+        d = sqrt(c[1])
         for j in 1:b+1
-            l[j] = r[j]/x
+            l[j] = c[j]/d
         end
         for j in 1:min(n-k+1, b+1)
             R[k, j+k-1] = l[j]
         end
-        F12 = k == 1 ? T(2) : T(1)
-        r2[1] = r[2]/F12
+        X21 = k == 1 ? T(2) : T(1)
+        ĉ[1] = c[2]/X21
         for j in 2:b+1
-            r2[j] = (r[j+1]+r[j-1] + r[1]*col1[j] - col1[1]*r[j])/F12
+            ĉ[j] = (c[j+1] + c[j-1] + c[1]*row1[j] - row1[1]*c[j])/X21
         end
-        r2[b+2] = (r[b+1] + r[1]*col1[b+2] - col1[1]*r[b+2])/F12
-        cst = r[2]/x
+        ĉ[b+2] = (c[b+1] + c[1]*row1[b+2] - row1[1]*c[b+2])/X21
+        cst = c[2]/d
         for j in 1:b+2
-            r[j] = r2[j+1] - cst*l[j+1]
+            c[j] = ĉ[j+1] - cst*l[j+1]
         end
+        cst = X21/d
         for j in 1:b+2
-            col1[j] = -F12/x*l[j+1]
+            row1[j] = -cst*l[j+1]
         end
     end
 end

src/libfasttransforms.jl

@@ -1250,3 +1250,16 @@ for (fC, T) in ((:execute_jacobi_similarityf, Float32), (:execute_jacobi_similar
         end
     end
 end
+
+function modified_jacobi_matrix(R, XP)
+    n = size(R, 1) - 1
+    XQ = SymTridiagonal(zeros(n), zeros(n-1))
+    XQ.dv[1] = (R[1, 1]*XP[1, 1] + R[1, 2]*XP[2, 1])/R[1, 1]
+    for i in 1:n-1
+        XQ.ev[i] = R[i+1, i+1]*XP[i+1, i]/R[i, i]
+    end
+    for i in 2:n
+        XQ.dv[i] = (R[i, i]*XP[i,i] + R[i, i+1]*XP[i+1, i] - XQ[i, i-1]*R[i-1, i])/R[i, i]
+    end
+    return XQ
+end