Constructing positive-definite Wendland Compact-Support RBF Matrix

Here is my scilab version. The resultingt matrix appears to be positive definite, so the problem must be with the matrix referencing in my C++ code.

###########################

M = load("C:/tmp2/PolySphereConstraints.txt");

n = length(M);

alpha = 1.0;



for i=1:n

    x1 = M(i,1);

    y1 = M(i,2);

    z1 = M(i,3);



    for j=1:n

        x2 = M(j,1);

        y2 = M(j,2);

        z2 = M(j,3);



        xdiff = x2 - x1;

        ydiff = y2 - y1;

        zdiff = z2 - z1;



          r = sqrt( xdiff^2 + ydiff^2 + zdiff^2 );          

        rr = r/alpha;



          if (rr < 1)

            phi(i,j) = (1-rr)^4*(4*rr+1);

           else

                 phi(i,j) = 0;

        endif



    end

end



eig(phi)

My matrix produced by my C++ code (without the constraint blocks) matches the result of my octave (not scilab) script i.e. both are positive definite. However, when I add the constraint blocks, the matrix is no longer positive definite. Below is my octave script:

 #################

M = load("C:/tmp2/PolySphereConstraints.txt");

n = length(M);

alpha = 3.5;



for i=1:n

    x1 = M(i,1);

    y1 = M(i,2);

    z1 = M(i,3);



    for j=1:n

        x2 = M(j,1);

        y2 = M(j,2);

        z2 = M(j,3);



        xdiff = x2 - x1;

        ydiff = y2 - y1;

        zdiff = z2 - z1;



          r = sqrt( xdiff^2 + ydiff^2 + zdiff^2 );         

        rr = r/alpha;



          if (rr < 1)

            matrix(i,j) = (1-rr)^4*(4*rr+1);

           else

                 matrix(i,j) = 0;

        endif

    end

    #constraints (refer to paper)

    matrix(n+1,i) = matrix(i,n+1) = x1;

    matrix(n+2,i) = matrix(i,n+2) = y1;

    matrix(n+3,i) = matrix(i,n+3) = z1;

    matrix(n+4,i) = matrix(i,n+4) = 1;

end



#zeros

for i=1:4

    for j=1:4

        matrix(n+i,n+j) = 0;    

    end

end



eig(matrix)

# chol(matrix) #alternatively