Chebyshev smoother, I cannot get this to run, I dont understand why.

dtVx = rho*resolx.^2./(4*eta_nbv*(1.+eta_b)*4.1); %dtVx is a very large vector
A = (4*eta_nbv*(1.+eta_b)*4.1); %Here I am using Ax=b for the above dtVx vector
A = diag(A);
[V,D] = eigs(A);
lam1 = min(diag(D))
lam2 = max(diag(D))
%END OF ARNOLDI
S = sparse(A);
S = S \ A;
b = rho*resolx.^2; %This b is from the dtVx vector above
r = S.*(b - A.*dtVx); %redidual
thet = 0.5*(lam2+lam1); %Theta value 
delt = 0.5*(lam2-lam1); %Delta value
sig = thet/delt; %sigma value
rho1 = 1/sig;
dk = r/thet;
xk = 0;
nbv = 44229;
for k = 1:nbv
    xk = xk + dk;
    r = r - S*A*dk;
    rhok = 1/(2*sig-rhok);
    dk = (rhok*rho1*dk + 2*rhok*r)/(delt+1);
      rho1 = rhok
      end
xk = xk + dk;

Answer 1

[Paul]

Hi Nathan,

Can't run the code because not all of the information needed to run it is provided.

dtVx = rho*resolx.^2./(4*eta_nbv*(1.+eta_b)*4.1); %dtVx is a very large vector

Unrecognized function or variable 'rho'.

Comparing the code to the Algorithm, it looks like there may be a few issues.

A = (4*eta_nbv*(1.+eta_b)*4.1); %Here I am using Ax=b for the above dtVx vector
A = diag(A);
[V,D] = eigs(A);
lam1 = min(diag(D))
lam2 = max(diag(D))
%END OF ARNOLDI
S = sparse(A);

This line should result in S being the identity matrix, to within numerical precision.

S = S \ A;

If that's the dersired result, why not just initialize S using eye.

b = rho*resolx.^2; %This b is from the dtVx vector above

This line using element-wise multiplication and not matrix multiplication. Is thar correct? Also, this value of r is used and updated in the loop, even though the Algorithm doesn't initialize it as r_1. Is that a typo in the Algorithm?

r = S.*(b - A.*dtVx); %redidual

thet = 0.5*(lam2+lam1); %Theta value 
delt = 0.5*(lam2-lam1); %Delta value
sig = thet/delt; %sigma value
rho1 = 1/sig;
dk = r/thet;
xk = 0;
nbv = 44229;
for k = 1:nbv
    xk = xk + dk;
    r = r - S*A*dk;
    rhok = 1/(2*sig-rhok);

The next line doesn't look correct. Accoding to the algorithm it should be

    dk = (rhok*rho1*dk + 2*rhok*r)/(delt+1);
    % dk = rhok*rho1*dk + 2*rhok*r/delt   % should be this?

    rho1 = rhok;
end
xk = xk + dk;