Why built-in conjugated gradient solver for sparse linear matrix runs much faster than normal script, no matter GPU or CPU?

I am trying to profile the conjugated gradient method for Ax = b. A is a large sparse matrix. The conjugated gradient method could be code in a .m file. The method is from wiki. But the speed of this two ways differs in a stable manner. The matlab build in function, pcg, without any preconditioner, is always faster than the .m file. Could anyone tell me the reason of it? Even when matrix A and vector b are gpuArray. The pcg still runs faster. Does Matlab use some way to speed up implicitly? My m file is as below. I had checked the answers x1~x4 are equal under tolerance and the iter number k1~k4 are close. Is there any thing wrong with my code? Any suggestion would be appreicated.

function [x,dist,k] = cg1(A,b,tol,maxit)
N = length(b);
norm_b = norm(b);
%% initial value
x = zeros(N,1);
r = b;
dist = norm(r)/norm_b;
if dist < tol
    return;
end
p = b;
k = 0;
while (dist>tol && k<maxit)
%     disp(['loop ' num2str(k)]);
    %% update variables
    r0 = r;
    x0 = x;
    p0 = p;
    %% evaluate step k+1
    Ap0 = A*p0;
    r02 = r0'*r0;
    alpha = r02/(p0'*Ap0);
    x = x0+alpha*p0;
    r = r0-alpha*Ap0;
    dist = norm(r)/norm_b;
    beta = r'*r/r02;
    p = r+beta*p0;
    k = k+1;
end

The profile code is

tol = 1e-9;
maxit = length(b);
%% CPU version 1 : build in pcg 
tic;[x1,~,dist1,k1]  = pcg(A,b,tol,maxit);toc; % 187s
%% CPU version 2 : self_made
tic;[x2,dist2,k2] = cg1(A,b,tol,maxit);toc; % 207.26s
%% gpuArray
Ag = gpuArray(A);
bg = gpuArray(b);
%% GPU version 1 : build in pcg
tic;[x3,dist3,k3] = cg1(Ag,bg,tol,maxit);toc; % 16.6s
%% GPU version 2 : self_made
tic;[x4,~,dist4,k4] = pcg(Ag,bg,tol,maxit);toc; % 14.6s

[Bjorn Gustavsson]

Try to use the profiler to compare/investigate the execution of your implementation and pcg. That should help you figure out what goes on where and for how long. Check the help and documentation for profile.

HTH