super slow loop

[smateria]

Hi everyone,
I have been trying to use a function to calculate standardized precipitation
index (SPI) that works fine in Matlab, but seems having huge problems in
Octave. In fact, after two integrations, the loop starts proceeding
incredibly slowly. I have vectorized my matrices as much as I could but the
problem persists. Here is the code:

[...]
size (obsmean)
ans =
   204   350

size (pre_mod)
ans =
   74   204   350

spi_obs = NaN(192,350); 
spi_mod_res = NaN(192,25900);  
pre_mod_res = reshape(pre_mod,204,25900); % reshaped from a 74ens x 204time
x 350latlon

for i = 1:350
    if isfinite(obsmean(:,i))==1;
        spi_obs(:,i) = SPI(obsmean(:,i),3,12);     
    end
end

for ens = 1:25900
    if isfinite(obsmean(:,i))==1;
        spi03mod_res(:,ens) = SPI(pre_mod_res(:,ens),3,12);
    end
end

The first loop takes around 50 seconds to complete, or about 30 if I
suppress a few printouts from the function nmsmax.m. 
The second loop slows down at the 2nd iteration, and then starts proceeding
superslowly (if I ^C after five minutes, I see we are at the third iteration
:O).

Here below, the SPI Matlab function by T.Lee (2009), perfectly working on
Matlab

/function [Z]=SPI(Data,scale,nseas)

% Standardized Precipitation Index 
% Input Data
% Data : Monthly Data vector not matrix (monthly or seasonal precipitation)
% scale : 1,3,12,48
% nseas : number of season (monthly=12)
% Example
% Z=SPI(gamrnd(1,1,1000,1),3,12); 3-monthly scale, 
% Notice that the rest of the months of the first year are removed.
% eg. if scale =3, fist year data 3-12 SPI values are not estimated.

% if row vector then make coloumn vector
% if (sz==1); Data(:,1) = Data;end

erase_yr = ceil(scale/12);

% Data setting to scaled dataset

A1 = [];
for is = 1:scale, A1=[A1,Data(is:length(Data)-scale+is)]; end
XS = sum(A1,2);

if (scale>1), XS(1:nseas*erase_yr-scale+1)=[]; end

for is = 1:nseas
    tind = is:nseas:length(XS);
    Xn = XS(tind);
    [zeroa] = find(Xn==0);
    Xn_nozero = Xn; Xn_nozero(zeroa)=[];
    q = length(zeroa)/length(Xn);
    parm = gamfit(Xn_nozero);
    Gam_xs = q+(1-q)*gamcdf(Xn,parm(1),parm(2));
    Z(tind) = norminv(Gam_xs);
end/



--
View this message in context: 
http://octave.1599824.n4.nabble.com/super-slow-loop-tp4683040.html
Sent from the Octave - General mailing list archive at Nabble.com.