function [C,Cinv,C_handle,Csecondpart,Cinv_gam_handle,Csecondpart_fast_handle,...
    C_handle_neg] = compute_Cmat_LRRS(Phi,W_1_small_c,W_1_diag_c,lambda,weight_op,weight_vec_c,d1,d2,fac)
Wvec = ((lambda).*(fac.*weight_vec_c-fac/weight_op.S_c_eps)).^(-1);
Bvec = (fac.*(lambda./weight_op.S_c_eps+W_1_diag_c)).^(-1);
Bvec_sm = ((lambda./(weight_op.S_c_eps/fac)+fac.*W_1_small_c)).^(-1);

USigU = weight_op.U'*(Bvec_sm.*weight_op.U);
USigU_inv = inv(USigU);
r = size(weight_op.U,2);


% Csecondpart_hdl = @(gam) matrixspace_to_tangspace(reshape(Bvec.*reshape(tangspace_to_matrixspace(gam,weight_op),d1*d2,1),d1,d2),weight_op);
% Csecondpart_action = @(Y) cell2mat(cellfun(Csecondpart_hdl,num2cell(Y,1),'UniformOutput',false));
% Csecondpart = Csecondpart_action(eye(length(weight_vec_c)));
Csecondpart = [];
%cell2mat(cellfun(A,num2cell(X,1),'UniformOutput',false));
% Cinv_gam_handle = @(gam) action_Cinv(gam,weight_op,USigU,USigU_inv,Bvec_sm,d1,d2,r);
Cinv_gam_handle  = [];
Csecondpart_fast_handle = @(gam) action_Csecondpart(gam,weight_op,USigU,USigU_inv,Bvec_sm,d1,d2,r);
% C_handle_slow = @(gam) Wvec.*gam...
 %  +matrixspace_to_tangspace(reshape(Bvec.*reshape(tangspace_to_matrixspace(gam,weight_op),d1*d2,1),d1,d2),weight_op);
C_handle = @(gam) Wvec.*gam+Csecondpart_fast_handle(gam);
C_handle_neg = [];%@(gam) -(Wvec.*gam...
  % +Csecondpart_fast_handle(gam));
handle_matrixaction = @(Y) cell2mat(cellfun(C_handle,num2cell(Y,1),'UniformOutput',false));
C = handle_matrixaction(eye(length(weight_vec_c)));
% Cinv = inv(C);
Cinv = [];

end


function Cinv_gam = action_Cinv(gam,weight_op,USigU,USigU_inv,Bvec_sm,d1,d2,r)

[~,Gam_1_tilde,Gam_3_tilde,Gam_2_tilde] = get_Tk_matrices(gam,weight_op);

Gam_2 = USigU_inv*Gam_2_tilde;

Tmp = Bvec_sm.^(-1/2).*Gam_3_tilde;
ProjMat = eye(d1)- Bvec_sm.^(1/2).*weight_op.U*(USigU_inv*(weight_op.U'.* Bvec_sm'.^(1/2)));
Tmp2 = ProjMat\Tmp;
Gam_3 =  Bvec_sm.^(-1/2).*Tmp2;
Gam_1 = USigU_inv*(Gam_1_tilde-weight_op.U'*(Bvec_sm.*Gam_3))+weight_op.U'*Gam_3;

Cinv_gam = zeros(length(gam),1);
    Cinv_gam(1:r^2)              = reshape(Gam_1,[r^2,1]);
    Cinv_gam((r^2+1):(r*(r+d1))) = reshape(Gam_3,[r*d1,1]);
    Cinv_gam((r*(r+d1)+1):end)   = reshape(Gam_2,[r*d2,1]);

    
end

function Csec_gam = action_Csecondpart(gam,weight_op,USigU,USigU_inv,Bvec_sm,d1,d2,r)

[~,Gam_1,Gam_3,Gam_2] = get_Tk_matrices(gam,weight_op);

Gam_2_tilde = USigU*Gam_2;
Tmp = weight_op.U'*Gam_3;
Gam_3_tilde = Bvec_sm.*(weight_op.U*Gam_1+Gam_3)-(Bvec_sm.*weight_op.U)*Tmp;
Gam_1_tilde = USigU*(Gam_1-Tmp)+weight_op.U'*(Bvec_sm.*Gam_3);
% Tmp = Bvec_sm.^(-1/2).*Gam_3_tilde;
% ProjMat = eye(d1)- Bvec_sm.^(1/2).*weight_op.U*(USigU_inv*(weight_op.U'.* Bvec_sm'.^(1/2)));
% Tmp2 = ProjMat\Tmp;
% Gam_3 =  Bvec_sm.^(-1/2).*Tmp2;
% Gam_1 = USigU_inv*(Gam_1_tilde-weight_op.U'*(Bvec_sm.*Gam_3))+weight_op.U'*Gam_3;

Csec_gam = zeros(length(gam),1);
Csec_gam(1:r^2)              = reshape(Gam_1_tilde,[r^2,1]);
Csec_gam((r^2+1):(r*(r+d1))) = reshape(Gam_3_tilde,[r*d1,1]);
Csec_gam((r*(r+d1)+1):end)   = reshape(Gam_2_tilde,[r*d2,1]);

    
end

%function C = construct_Cmat