function[x, x_d, hat_Q1, sigma_d, ifcvg] = cVAMPa_dampling(y, A, lambda, alpha, delta, iter_max, sigma)

[M, N] = size(A);

gamma = M/N;
p = A'*y;
h1 = p/gamma;
hat_Q1 = gamma;

%Eigenvalue Decomposition
[V, D] = eig(A'*A);
d = diag(D); 

t = 0;
diff = 1;

while((diff > delta) && (t < iter_max))
    
    h1_pre = h1;
    
    % Factorized
    x1 = sft_thd(h1, lambda/hat_Q1);
    chi1 = sum((abs(x1) > 1e-4).* (2 - lambda./((hat_Q1*abs(x1) + lambda)))) / 2 / N / hat_Q1;
    % Message F to G
    hat_Q2 = 1/chi1 - hat_Q1;
    h2 = (x1/chi1 - h1*hat_Q1)/hat_Q2;
    % Gaussian
    tmp = V'*(p + h2*hat_Q2);
    tmp = tmp./(d+hat_Q2);
    x2 = V*tmp;
    chi2 = sum(1./(d+hat_Q2))/N;
    % Message G to F
    hat_Q1 = 1/chi2 - hat_Q2;
    h1 = alpha*(x2/chi2 - h2*hat_Q2)/hat_Q1+(1-alpha)*h1;
    
    diff = sum(abs(h1_pre - h1))/sum(abs(h1));
    t = t+1;
    
end

ifcvg = diff <= delta;

x = x1;
x_d = h1;

chi = chi1;
t = -hat_Q2;
t_prime = -1/mean((1./(d+hat_Q2)).^2);
G_prime = t + 1/chi;
G_wprime = t_prime + 1/chi/chi;
RSS = sum(abs(y - A*x).^2)/M;
hat_chi = gamma*G_wprime/(2*G_prime-2*chi*G_wprime)*RSS +...
    (-G_wprime*gamma+G_prime*G_prime)/(2*G_prime-2*chi*G_wprime)*sigma^2;
sigma_d = sqrt(2*hat_chi)/hat_Q1;

end