fasterfcsummation_matlab.m

RankNGG 算法实现. 含有dll文件。源码为matlab

function [E]=FastErfcSummation_MATLAB(N,M,q,x,y,p,h,rx,r,nn)
% Code for fast erfc summation.
%
% MATLAB  Implementation.
%
% Parameters for the algorithm. Use the MATLAB function [rx,r,h,p,n]=choose_parameters(epsil)
%
%% Input
%
% * N ...number of source points.
% * M ... number of target points.
% * q ... 1 x N matrix of N source weights.
% * x ... 1 x N matrix of N source points.
% * y ... 1 x M matrix of M target points.
%
% Algorithm parameters--
% 
% * rx ... the interval length
% * r ... the cutoff radius
% * h ... the series parameter
% * p ... the truncation number
% * n ... number of influential clusters
%
%% Ouput
%
% *E ... 1 x M vector of the evaluated sum .
%
%
%% Signature
%
% Author: Vikas Chandrakant Raykar
% E-Mail: [e-mail]
% Date: September 27, 2006
%
%% See also
%
%  choose_parameters

min_x=min(min(x),min(y));
max_x=max(max(x),max(y));
K=ceil((max_x-min_x)/(2*rx));
%max_x=(2*rx*K)+min_x;
cluster_width=2*rx;
cluster_radius=rx;

cluster_center=[ ];
for k=1:K
    cluster_center(k)=min_x+((k-1)*cluster_width)+cluster_radius;
end

h_square=h*h;
upper_limit=(2*p)-1;
const=4/pi;

coeffs=zeros(1,upper_limit);
A=zeros(K,upper_limit);
B=A;
cluster_Q=zeros(1,K);

for n=1:2:upper_limit
    coeffs(n)=exp(-n*n*h_square)/n;
end

for i=1:N
    cluster_index=max(ceil((x(i)-min_x)/cluster_width),1);
    x_star=cluster_center(cluster_index);
    for n=1:2:upper_limit
        A(cluster_index,n)=A(cluster_index,n)+(q(i)*cos(2*n*h*(x(i)-x_star)));
        B(cluster_index,n)=B(cluster_index,n)+(q(i)*sin(2*n*h*(x(i)-x_star)));
    end
    cluster_Q(cluster_index)=cluster_Q(cluster_index)+q(i);
end

for k=1:K
    for n=1:2:upper_limit
        A(k,n)=coeffs(n)*A(k,n);
        B(k,n)=coeffs(n)*B(k,n);
    end
end

G=zeros(1,M);
E=G;

for j=1:M
    G(j)=0;
    Q=0;
    cluster_index=max(ceil((y(j)-min_x)/cluster_width),1);
    low_limit=max(1,cluster_index-nn);
    upp_limit=min(K,cluster_index+nn);
    for k=low_limit:upp_limit
        x_star=cluster_center(k);    
       for n=1:2:upper_limit
            G(j)=G(j)+((A(k,n)*sin(2*n*h*(y(j)-x_star)))-(B(k,n)*cos(2*n*h*(y(j)-x_star))));
        end
        Q=Q+cluster_Q(k);
    end
    E(j)=Q-(const*G(j));
    for k=upp_limit+1:K
          E(j)=E(j)+(2.0*cluster_Q(k));
    end  
end