Incorporating diagonal scaling into partial Cholesky fixes cheb2legTH

Author: MikaelSlevinsky

docs/api/FastTransforms.md

@@ -13,7 +13,7 @@ See also [`icjt`](#method__icjt.1) and [`jjt`](#method__jjt.1).
 
 
 *source:*
-[FastTransforms/src/cjt.jl:121](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/cjt.jl#L121)
+[FastTransforms/src/cjt.jl:121](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/cjt.jl#L121)
 
 ---
 
@@ -23,7 +23,7 @@ Calculates the Gaunt coefficients in 64-bit floating-point arithmetic.
 
 
 *source:*
-[FastTransforms/src/gaunt.jl:24](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/gaunt.jl#L24)
+[FastTransforms/src/gaunt.jl:24](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/gaunt.jl#L24)
 
 ---
 
@@ -43,7 +43,7 @@ This is a Julia implementation of the stable recurrence described in:
 
 
 *source:*
-[FastTransforms/src/gaunt.jl:14](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/gaunt.jl#L14)
+[FastTransforms/src/gaunt.jl:14](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/gaunt.jl#L14)
 
 ---
 
@@ -56,7 +56,7 @@ See also [`cjt`](#method__cjt.1) and [`jjt`](#method__jjt.1).
 
 
 *source:*
-[FastTransforms/src/cjt.jl:129](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/cjt.jl#L129)
+[FastTransforms/src/cjt.jl:129](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/cjt.jl#L129)
 
 ---
 
@@ -69,7 +69,7 @@ See also [`cjt`](#method__cjt.1) and [`icjt`](#method__icjt.1).
 
 
 *source:*
-[FastTransforms/src/cjt.jl:137](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/cjt.jl#L137)
+[FastTransforms/src/cjt.jl:137](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/cjt.jl#L137)
 
 ---
 
@@ -88,7 +88,7 @@ Optionally:
 
 
 *source:*
-[FastTransforms/src/cjt.jl:158](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/cjt.jl#L158)
+[FastTransforms/src/cjt.jl:158](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/cjt.jl#L158)
 
 ---
 
@@ -107,7 +107,7 @@ Optionally:
 
 
 *source:*
-[FastTransforms/src/cjt.jl:172](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/cjt.jl#L172)
+[FastTransforms/src/cjt.jl:172](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/cjt.jl#L172)
 
 ## Internal
 
@@ -122,7 +122,7 @@ Modified Chebyshev moments of the first kind with respect to the Jacobi weight:
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:399](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L399)
+[FastTransforms/src/specialfunctions.jl:399](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L399)
 
 ---
 
@@ -135,7 +135,7 @@ Modified Chebyshev moments of the second kind with respect to the Jacobi weight:
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:417](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L417)
+[FastTransforms/src/specialfunctions.jl:417](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L417)
 
 ---
 
@@ -145,7 +145,7 @@ Compute weights of the Clenshaw—Curtis quadrature rule with a Jacobi weight.
 
 
 *source:*
-[FastTransforms/src/clenshawcurtis.jl:12](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/clenshawcurtis.jl#L12)
+[FastTransforms/src/clenshawcurtis.jl:12](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/clenshawcurtis.jl#L12)
 
 ---
 
@@ -155,7 +155,7 @@ Compute nodes and weights of the Clenshaw—Curtis quadrature rule with a Jacobi
 
 
 *source:*
-[FastTransforms/src/clenshawcurtis.jl:6](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/clenshawcurtis.jl#L6)
+[FastTransforms/src/clenshawcurtis.jl:6](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/clenshawcurtis.jl#L6)
 
 ---
 
@@ -165,7 +165,7 @@ Compute Jacobi expansion coefficients in Pₙ^(α-1,β) given Jacobi expansion c
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:467](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L467)
+[FastTransforms/src/specialfunctions.jl:467](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L467)
 
 ---
 
@@ -175,7 +175,7 @@ Compute Jacobi expansion coefficients in Pₙ^(α-1,α-1) given Jacobi expansion
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:489](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L489)
+[FastTransforms/src/specialfunctions.jl:489](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L489)
 
 ---
 
@@ -185,7 +185,7 @@ Compute Jacobi expansion coefficients in Pₙ^(α,β-1) given Jacobi expansion c
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:478](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L478)
+[FastTransforms/src/specialfunctions.jl:478](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L478)
 
 ---
 
@@ -195,7 +195,7 @@ Compute nodes and weights of Fejer's first quadrature rule with a Jacobi weight.
 
 
 *source:*
-[FastTransforms/src/fejer.jl:7](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/fejer.jl#L7)
+[FastTransforms/src/fejer.jl:7](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/fejer.jl#L7)
 
 ---
 
@@ -205,7 +205,7 @@ Compute nodes and weights of Fejer's second quadrature rule with a Jacobi weight
 
 
 *source:*
-[FastTransforms/src/fejer.jl:12](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/fejer.jl#L12)
+[FastTransforms/src/fejer.jl:12](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/fejer.jl#L12)
 
 ---
 
@@ -215,7 +215,7 @@ Compute weights of Fejer's first quadrature rule with a Jacobi weight.
 
 
 *source:*
-[FastTransforms/src/fejer.jl:21](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/fejer.jl#L21)
+[FastTransforms/src/fejer.jl:21](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/fejer.jl#L21)
 
 ---
 
@@ -225,7 +225,7 @@ Compute weights of Fejer's second quadrature rule with a Jacobi weight.
 
 
 *source:*
-[FastTransforms/src/fejer.jl:26](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/fejer.jl#L26)
+[FastTransforms/src/fejer.jl:26](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/fejer.jl#L26)
 
 ---
 
@@ -235,7 +235,7 @@ Compute a typed 0.5.
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:12](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L12)
+[FastTransforms/src/specialfunctions.jl:12](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L12)
 
 ---
 
@@ -245,7 +245,7 @@ Compute Jacobi expansion coefficients in Pₙ^(α+1,β) given Jacobi expansion c
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:432](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L432)
+[FastTransforms/src/specialfunctions.jl:432](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L432)
 
 ---
 
@@ -255,7 +255,7 @@ Compute Jacobi expansion coefficients in Pₙ^(α+1,α+1) given Jacobi expansion
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:454](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L454)
+[FastTransforms/src/specialfunctions.jl:454](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L454)
 
 ---
 
@@ -265,7 +265,7 @@ Compute Jacobi expansion coefficients in Pₙ^(α,β+1) given Jacobi expansion c
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:443](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L443)
+[FastTransforms/src/specialfunctions.jl:443](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L443)
 
 ---
 
@@ -275,7 +275,7 @@ Pochhammer symbol (x)_n = Γ(x+n)/Γ(x) for the rising factorial.
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:32](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L32)
+[FastTransforms/src/specialfunctions.jl:32](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L32)
 
 ---
 
@@ -285,7 +285,7 @@ Stirling series for Γ(z).
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:63](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L63)
+[FastTransforms/src/specialfunctions.jl:63](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L63)
 
 ---
 
@@ -295,7 +295,7 @@ Compute a typed 2.
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:20](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L20)
+[FastTransforms/src/specialfunctions.jl:20](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L20)
 
 ---
 
@@ -307,7 +307,7 @@ For 64-bit floating-point arithmetic, the Lambda function uses the asymptotic se
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:147](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L147)
+[FastTransforms/src/specialfunctions.jl:147](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L147)
 
 ---
 
@@ -317,7 +317,7 @@ The Lambda function Λ(z) = Γ(z+½)/Γ(z+1) for the ratio of gamma functions.
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:141](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L141)
+[FastTransforms/src/specialfunctions.jl:141](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L141)
 
 ---
 
@@ -327,7 +327,7 @@ The Lambda function Λ(z,λ₁,λ₂) = Γ(z+λ₁)/Γ(z+λ₂) for the ratio of
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:160](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L160)
+[FastTransforms/src/specialfunctions.jl:160](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L160)
 
 ---
 
@@ -337,7 +337,7 @@ The Kronecker δ function.
 
 
 *source:*
-[FastTransforms/src/specialfunctions.jl:26](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/specialfunctions.jl#L26)
+[FastTransforms/src/specialfunctions.jl:26](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/specialfunctions.jl#L26)
 
 ---
 
@@ -351,5 +351,5 @@ where the Hankel matrix `H` is non-negative definite. This allows a Cholesky dec
 
 
 *source:*
-[FastTransforms/src/toeplitzhankel.jl:8](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/9a0fff4d389746d63718883170537e95ad849e76/src/toeplitzhankel.jl#L8)
+[FastTransforms/src/toeplitzhankel.jl:8](https://github.com/MikaelSlevinsky/FastTransforms.jl/tree/b6e4fdf35a2160f0d9948e6565fd69c2177f4d82/src/toeplitzhankel.jl#L8)
 

src/toeplitzhankel.jl

@@ -19,8 +19,10 @@ function ToeplitzHankelPlan(T::TriangularToeplitz,C::Vector,DL::AbstractVector,D
 end
 ToeplitzHankelPlan(T::TriangularToeplitz,C::Matrix) =
     ToeplitzHankelPlan(T,C,ones(size(T,1)),ones(size(T,2)))
-ToeplitzHankelPlan(T::TriangularToeplitz,H::Hankel,D...) =
-    ToeplitzHankelPlan(T,partialchol(H),D...)
+ToeplitzHankelPlan(T::TriangularToeplitz,H::Hankel,DL::AbstractVector,DR::AbstractVector) =
+    ToeplitzHankelPlan(T,partialchol(H),DL,DR)
+ToeplitzHankelPlan(T::TriangularToeplitz,H::Hankel,D::AbstractVector,DL::AbstractVector,DR::AbstractVector) =
+    ToeplitzHankelPlan(T,partialchol(H,D),DL,DR)
 
 *(P::ToeplitzHankelPlan,v::AbstractVector)=P.DL.*toeplitzcholmult(P.T,P.C,P.DR.*v)
 
@@ -50,6 +52,33 @@ function partialchol(H::Hankel)
     C
 end
 
+function partialchol(H::Hankel,D::AbstractVector)
+    # Assumes positive definite
+    T = promote_type(eltype(H),eltype(D))
+    σ=T[]
+    n=size(H,1)
+    C=Vector{T}[]
+    v=[H[:,1];vec(H[end,2:end])]
+    d=diag(H).*D.^2
+    @assert length(v) ≥ 2n-1
+    reltol=maxabs(d)*eps(T)*log(n)
+    for k=1:n
+        mx,idx=findmax(d)
+        if mx ≤ reltol break end
+        push!(σ,inv(mx))
+        push!(C,v[idx:n+idx-1].*D.*D[idx])
+        for j=1:k-1
+            nCjidxσj = -C[j][idx]*σ[j]
+            Base.axpy!(nCjidxσj, C[j], C[k])
+        end
+        @simd for p=1:n
+            @inbounds d[p]-=C[k][p]^2/mx
+        end
+    end
+    for k=1:length(σ) scale!(C[k],sqrt(σ[k])) end
+    C
+end
+
 function toeplitzcholmult(T,C,v)
     n,K = length(v),length(C)
     ret,temp1,temp2 = zero(v),zero(v),zero(v)
@@ -86,9 +115,10 @@ function cheb2legTH{S}(::Type{S},n)
     T = TriangularToeplitz(t,:U)
     h = Λ(1:half(S):n-1,zero(S),3half(S))
     H = Hankel(h[1:n-1],h[n-1:end])
-    DL = 3half(S):n-half(S)
-    DR = -(one(S):n-one(S))/4
-    T,H,DL,DR
+    D = 1:one(S):n-1
+    DL = (3half(S):n-half(S))./D
+    DR = -(one(S):n-one(S))./4D
+    T,H,D,DL,DR
 end
 
 function leg2chebuTH{S}(::Type{S},n)