function [x] = conjgrad_trig(A, b)
    [Nx,Nv,Nz] = size(A);
    x = zeros(Nx,Nv,Nz);
    r = b - Trig_Multi(A,x);
    p = r;
    rsold = sum(r .* r,2);

    for i = 1:length(b)
        Ap = Trig_Multi(A,p);
        alpha = rsold ./ sum(p .* Ap,2);
        x = x + repmat(alpha,[1,Nv,1]) .* p;
        r = r - repmat(alpha,[1,Nv,1]) .* Ap;
        rsnew = sum(r .* r,2);
        temp = max(reshape(rsnew,Nx*Nz,1));
        if sqrt(temp) < 1e-10
              break;
        end
        p = r + repmat(rsnew ./ rsold,[1,Nv,1]) .* p;
        rsold = rsnew;
    end
end