Need help/advise constructing positive-definite Wendland Compact-Support RBF Matrix

Hello,

I am trying to construct a matrix using Wendland's Compact-Support RBF as described in the following the paper, http://www.cs.umbc.edu/~rheingan/pubs/smi2001.pdf . In particular, I'm using the C2, d = 3 RBF: (1 - r)^4*(4r + 1). It seems quite simple enough, but none of the matrices I've constructed are positive definite, although Wendland's RBFs guarantee positive definite matrices -- the sparse solver I intend to use requires positive definite matrices.

Below are the relevant portions of my code in which I initialize my matrix. Is there anything I'm doing wrong?

Thanks,

- Olumide

//////////////////////////// code ////////////////////////////

for( unsigned i = 0; i < constraintCount; ++i )
{
    for( unsigned j = 0; j < i; ++j )
    {
        rbfMatrix[i][j] = rbfMatrix[j][i] = computeRBF( j , i );
    }

    rbfMatrix[i][i] = 1.0;
    rbfMatrix[constraintCount+0][i]= rbfMatrix[i][constraintCount+0] = constraint[i].x;
    rbfMatrix[constraintCount+1][i]= rbfMatrix[i][constraintCount+1] = constraint[i].y;
    rbfMatrix[constraintCount+2][i]= rbfMatrix[i][constraintCount+2] = constraint[i].z;
    rbfMatrix[constraintCount+3][i]= rbfMatrix[i][constraintCount+3] = 1.0;
}

double computeRBF( unsigned i, unsigned j )
{
    double radius2 = ( constraint[i].x - constraint[j].x )*( constraint[i].x - constraint[j].x )  + ( constraint[i].y - constraint[j].y )*( constraint[i].y - constraint[j].y ) + ( constraint[i].z - constraint[j].z )*( constraint[i].z - constraint[j].z );
    double radius = sqrt(radius2);

    double r = radius/alpha;
    return ( r < 1.0 ) ? pow( (1 - r) , 4 ) * (4*r + 1) : 0;
}