1; # octave script file
N = 5;
D = 2;
EPS = 0.02;
X = rand(N, D);
Pi = zeros(N,1);
qEPS = EPS^2 / 4;
# pre-compute the values Pi:
for i = 1:rows(X)
# BLAS level 1 dsdot:
Pi(i) = qEPS - ( dot( X(1,:), X(1,:) ) / 2);
endfor
# BLAS level 3 dgemm:
scalarp = X * X';
# loop over lower triangle of scalarp
for i = 2:rows(scalarp)
for j = 1:(i-1)
# join partner i and j
scalarp(i,j) += Pi(i) + Pi(j);
endfor
endfor
# if scalarp(i,j) greater than 0, then we have a join partner
scalarp
