Do I implement this algorithm correctly?

I'm trying to implement the a block conjugate gradient algorithm that is not subject to breakdown from non invertible residual matrices; But I'm getting nonsensical results (in each iteration, the rank of Rcurrent should be getting smaller, not increasing). It is presented in the paper "A breakdown-free block conjugate gradient method" by Hao Ji and Yaohang Li.

Here is the algorithm:

This is my implementation in Julia:

function orth(M::Matrix)
  matrixRank = rank(M)
  Ufactor = svdfact(M)[:U]
  return Ufactor[:,1:matrixRank]
end

function BFBCG(A::Matrix, Xcurrent::Matrix, M::Matrix, tol::Number, maxit::Number, Rcurrent::Matrix)
    # initialization
    #Rcurrent = B - A*Xcurrent;
    Zcurrent = M*Rcurrent;
    Pcurrent = orth(Zcurrent);

    Xnext::Matrix = ones(size(Xcurrent))
    # iterative method
    for i = 0:maxit
        Qcurrent = A*Pcurrent
        acurrent =  (Pcurrent' * Qcurrent)\(Pcurrent'*Rcurrent)
        Xnext = Xcurrent+Pcurrent*acurrent
        Rnext = Rcurrent-Qcurrent*acurrent
        # if Residual norm of columns in Rcurrent < tol, stop
        Znext = M*Rnext
        bcurrent = -(Pcurrent' * Qcurrent)\ (Qcurrent'*Znext)
        Pnext = orth(Znext+Pcurrent*bcurrent)

        Xcurrent = Xnext
        Zcurrent = Znext
        Rcurrent = Rnext
        Pcurrent = Pnext
        @printf("\nRANK:\t%d",rank(Rcurrent))
        @printf("\nNORM column1:\t%1.8f",vecnorm(Rcurrent[:,1]))
        @printf("\nNORM column2:\t%1.8f\n=============",vecnorm(Rcurrent[:,2]))
    end
    return Xnext
end

The results of the paper for those inputs:

A = [15 5 4 3 2 1; 5 35 9 8 7 6; 4 9 46 12 11 10; 3 8 12 50 14 13; 2 7 11 14 19 15; 1 6 10 13 15 45]
M = eye(6)
guess = rand(6,2)
R0 = [1 0.537266261211281;2 0.043775211060964;3 0.964458562037146;4 0.622317517840541;5 0.552735938776748;6 0.023323943544997]
X = BFBCG(A,guess,M,tol,9,R0)

are a rank that reaches zero in the third iteration.