function P5Grade()

clc
n = 100000;
r = 100;
n = 1000;
r = 5;
P = randn(n,r);
Q = randn(n,r);
R = randn(n,r);
S = randn(n,r);
d = randn(n,1);
x = randn(n,1);

tic
yYOU = FastProd(P,Q,R,S,d,x);
tYOU = toc;

tic
yCVL1 = FastProdCVL1(P,Q,R,S,d,x);
tCVL1 = toc;
tic
yCVL2 = FastProdCVL2(P,Q,R,S,d,x);
tCVL2 = toc;

fprintf('tYOU/tCVL1 = %6.3f    ||yYOU - yCVL1 || = %10.3e\n',...
                       tYOU/tCVL1,norm(yYOU-yCVL1))
fprintf('tYOU/tCVL2 = %6.3f    ||yYOU - yCVL2 || = %10.3e\n',...
                       tYOU/tCVL2,norm(yYOU-yCVL2))

   function y = FastProdCVL2(P,Q,R,S,d,x)
 % P,Q,R, and S are n-by-r and n is a multiple of r
 % d and x are column n-vectors
 % y = Ax where A = tril(R*S',-1) + diag(d) + triu(P*Q',1)
 [n,r] = size(P);
 m = n/r;
 yU = zeros(n,1);
 yL = zeros(n,1);
 for k=1:r:n
     idx = k:k+r-1;
     yU(idx) = Q(idx,:)'*x(idx);
     yL(idx) = S(idx,:)'*x(idx);
 end
 for k=n-2*r+1:-r:1
     idx = k:k+r-1;
     yU(idx) = yU(idx)+yU(idx+r);
 end
 for k=r+1:r:n
     idx = k:k+r-1;
     yL(idx) = yL(idx)+yL(idx-r);
 end
 for k=1:r:n
     idx = k:k+r-1;
     yU(idx) = P(idx,:)*yU(idx);
     yL(idx) = R(idx,:)*yL(idx);
 end
 
 for k=1:r:n
     idx = k:k+r-1;
     yU(idx)= yU(idx)-tril(P(idx,:)*Q(idx,:)')*x(idx);
     yL(idx)= yL(idx)-triu(R(idx,:)*S(idx,:)')*x(idx);
 end
 y = yL + d.*x + yU;

 
 
   function y = FastProdCVL0(P,Q,R,S,d,x)
 % P,Q,R, and S are n-by-r and n is a multiple of r
 % d and x are column n-vectors
 % y = Ax where A = tril(R*S',-1) + diag(d) + triu(P*Q',1)
 y = (tril(R*S',-1) + diag(d) + triu(P*Q',1))*x;
 
 
    function y = FastProdCVL1(P,Q,R,S,d,x)
 % P,Q,R, and S are n-by-r and n is a multiple of r
 % d and x are column n-vectors
 % y = Ax where A = tril(R*S',-1) + diag(d) + triu(P*Q',1)
 [n,r] = size(P);
 y = d.*x;
 for k=1:r
     y = y + FastProdv(P(:,k),Q(:,k),R(:,k),S(:,k),x);
 end

 
   function y = FastProdv(p,q,r,s,x)
 % p,q,r,s, and x are column n-vectors.
 % y = Ax where A = tril(r*s',-1)  + triu(p*q',1)
 
 n = length(x);
 % Compute u = triu(pq')*x
 u = q.*x; u0 = u;
 for k=n-1:-1:1
     u(k) = u(k)+u(k+1);
 end
 % Remove the diagonal contribution...
 u = p.*u - p.*u0;
 
 % Compute v = tril(rs')*x
 v = s.*x; v0 = v;
 for k=2:n
     v(k) = v(k)+v(k-1);
 end
 % Remove the diagonal contribution...
 v = r.*v - r.*v0;
 % Combine
 y = u + v;
 