Using Clenshaw to build multmats

[ajt60gaibb]

There is the possibility of using Clenshaw's algorithm to build multmats. For example, here is the code of ChebT (I put this inside +util for a fair comparison):

function M = multmat1d_chebT_clenshaw( n, a, Mx )
%MULTMAT1D_CHEBT_CLENSHAW   Compute the 1D N x N multiplication matrix for the function
%
%   f(x) = sum_j a(j) T_j(x)
%
%  Alex Townsend, April 2020.

m = 2*n;
Mold = 0*Mx;
Ms = 0*Mx;
twoMx = 2*Mx;
I = speye(m); 
for k = n:-1:2
    Mnew = a(k)*I + twoMx*Ms - Mold;
    Mold = Ms; Ms = Mnew;
end
M = a(1)*I + Mx*Ms - Mold; 
M = M(1:n,1:n);
end

It can be faster than the three-term recurrence in certain situations.

n = 100; m = 2*n; 
Mx = spdiags(ones(m,2)/2, [-1 1], m, m);
Mx(2,1) = 1;

r = 10; 
a = [randn(r,1); zeros(n-r,1)];

tic
for j = 1:100
    M1 = util.multmat1d_chebT(n, a, Mx);
end
toc

tic
for j = 1:100
    M2 = util.multmat1d_chebT_clenshaw(n, a, Mx);
end
toc
normest(M1 - M2)
Elapsed time is 0.278272 seconds.
Elapsed time is 0.133217 seconds.
ans =
     1.122498678062323e-15

The precise speedup depends on n and r. Do we want to use clenshaw instead of three-term recurrence?

[ajt60gaibb]

Here's the code for multmat1d_ultraS using Clenshaw's algorithm:

function M = multmat1d_ultraS_clenshaw(n, a, Mx, lambda)
%MULTMAT1D_ULTRAS   Compute the 1D N x N multiplication matrix for the function
%
%   f(x,y) = sum_j a(j) C^{(lambda)}_j(x)
%
%  Alex Townsend, April 2020.

m = 2*n;
Mold = 0*Mx;
Ms = 0*Mx;
twoMx = 2*Mx;
I = speye(m); 
for k = n:-1:2
    Mnew = a(k)*I + (k+lambda-1)/(k)*twoMx*Ms - (k+2*lambda-1)/(k+1)*Mold;
    Mold = Ms; Ms = Mnew;
end
M = a(1)*I + lambda*twoMx*Ms - lambda*Mold; 
M = M(1:n,1:n);
end