Nearly identical code, hugely different runtimes

Hi all,

I've got a weird situation going on. I have the following code:

if true
for ii = 2:natom 
  tcheck = 1; 
         while tcheck ~= 0
             tcheck = 0; 
             lx =  rand*lxmax;
             ly =  rand*lymax;
             lz =  rand*lzmax;

               for jj = 1:ii-1

                   x12 = lx-x(jj);
                   y12 = ly-y(jj);
                   z12 = lz-z(jj);
            %
                   x12 = x12-lxmax* round(x12/lxmax);
                   y12 = y12-lymax* round(y12/lymax);
                   z12 = z12-lzmax* round(y12/lzmax);
            %
                   dist = sqrt(x12^2 +y12^2 + z12^2);
                   if dist <= rmax
                        tcheck =tcheck +1; 
                        break
                   end
               end
           end
           x(ii) = lx;
           y(ii) = ly;
           z(ii) = lz;  
    ii
  endif true

natom=4000 and x,y,z are vectors of length 4000. lxmax, lymax, and lzmax are all equal at 40.5 which is what I have "cubesize" set to in the next code. In general you can assume the parameters controlling the iteration process are the same for both.

I am trying to incorporate this code into another program and rewrote it as such:

if true
for i=2:natoms
        while flag ~= 0
            flag=0;
            xtemp=rand*cubesize;
            ytemp=rand*cubesize;
            ztemp=rand*cubesize;

              for j=1:i-1
                  xvec=xtemp-store(j,3);
                  yvec=ytemp-store(j,4);
                  zvec=ztemp-store(j,5);

                  xvec = xvec-cubesize* round(xvec/cubesize);
                  yvec = zvec-cubesize* round(yvec/cubesize);
                  zvec = zvec-cubesize* round(zvec/cubesize);

                  dist=sqrt(xvec^2+yvec^2+zvec^2);

                  if dist <= moldia
                      flag=1;
                      break;
                  end
              end
          end
      store(i,2)=i;
      store(i,3)=xtemp;
      store(i,4)=ytemp;
      store(i,5)=ztemp;
      flag=1;
      i
  end

endif true

This code takes a HUGE amount more time to run. What gives? One speeds its way to i=3500 or so in a matter of seconds. The other craws to i=700 or so and then bogs. I even tried swapping out and storing everything in single vectors first with no effect.

The first code was inside a function, so I tried putting my entire new script inside a dummy function to no effect.

Thanks, Nathan

[per isakson]

See Profile to Improve Performance

[Nathan Phelps]

Thanks Per. I'll look into that.

[Nathan Phelps]

For clarity, I have included the entirely of both scripts below:

Script/Function 1:

if true
 % This program generate the 
function initialdist4
clear all; 
close all;

natom = 4000; % number of atoms 
%lxmax = 54.9521; % box length of a cubic box
lxmax = 41.5; % box length of a cubic box
%
mbin = 80; 
rmax = 2.3; 
lymax = lxmax;
lzmax = lxmax;
dist = 0; 
tcheck = 0; 

fl1 = fopen ('loc1.txt','w');
fl2 = fopen ('rdf1.txt','w');
%
for ii = 1:natom
    x(ii) = 0.0;
    y(ii) = 0.0;
    z(ii) = 0.0; 
end
%
x(1) = rand*lxmax;
y(1) = rand*lymax;
z(1) = rand*lzmax; 
fprintf(fl1,'%d %d %f %f %f\n',1,1,x(1),y(1),z(1));
% 
tic
for ii = 2:natom 
    tcheck = 1; 
           while tcheck > 0
               tcheck = 0; 
               lx =  rand*lxmax;
               ly =  rand*lymax;
               lz =  rand*lzmax;

                 for jj = 1:ii-1

                     x12 = lx-x(jj);
                     y12 = ly-y(jj);
                     z12 = lz-z(jj);
              %
                     x12 = x12-lxmax* round(x12/lxmax);
                     y12 = y12-lymax* round(y12/lymax);
                     z12 = z12-lzmax* round(y12/lzmax);
              %
                     dist = sqrt(x12^2 +y12^2 + z12^2);
                     if dist <= rmax
                          tcheck =tcheck +1; 
                          break
                     end
                 end
             end
             x(ii) = lx;
             y(ii) = ly;
             z(ii) = lz;  
      fprintf(fl1,'%d %d %f %f %f\n',ii,1,x(ii),y(ii),z(ii));
      ii
  end 
  toc
  end

Second script (I have tried making this a function to no effect):

if true
%%Nathan Phelps
% Molecule Space Filler

% This may not be necessary in the final but simplifies debugging
clc;      
clear all;
close all;

avogadros = 6.022e23; % natoms/mol
density = 2.7; % g/cm^3
molarmass = 26.9815386; % g/mol;
molperang = (1/molarmass)*(density)*(1e-24)*avogadros; % moles/angstrom^3

natoms = 4000;
atomdia=2.3;
cubesize = (natoms/molperang)^(1/3);

store=zeros(natoms,5);
store(:,1)=1;
store(:,2)=1;
store(1,3)=rand*cubesize;
store(1,4)=rand*cubesize;
store(1,5)=rand*cubesize;

flag=1;

tic
for i=2:natoms
        while flag ~= 0
            flag=0;
            xtemp=rand*cubesize;
            ytemp=rand*cubesize;
            ztemp=rand*cubesize;

              for j=1:i-1
                  xvec=xtemp-store(j,3);
                  yvec=ytemp-store(j,4);
                  zvec=ztemp-store(j,5);

                  xvec = xvec-cubesize* round(xvec/cubesize);
                  yvec = zvec-cubesize* round(yvec/cubesize);
                  zvec = zvec-cubesize* round(zvec/cubesize);

                  dist=sqrt(xvec^2+yvec^2+zvec^2);

                  if dist <= atomdia
                      flag=1;
                      break;
                  end
              end
          end
      store(i,2)=i;
      store(i,3)=xtemp;
      store(i,4)=ytemp;
      store(i,5)=ztemp;
      flag=1;
      i
  end
  toc
  end

My intent was to utilize the intent of the code (Fill a 3d space with random x,y,z points with a given spacing) while removing some of the less than ideal practices (dynamic allocation, preallocating using a loop instead of zeros() or similar) but instead I managed to bog the code down somehow. Note that I removed the fprintf stuff from the second script and instead did output in a second section of the code later on. I figured it was better to do one "big" file write at the end instead of continually accessing the file and writing to it as the loop progresses.

[Nathan Phelps]

Welp. That's embarrassing.

if true
                yvec = zvec-cubesize* round(yvec/cubesize);
end

Answer 1

[Walter Roberson]

You should be pre-allocating the storage of the variables you are assigning in to. You might not use them all, but you can calculate a maximum possible occupancy and you can shrink the variables afterwards.

[Nathan Phelps]

Sorry, I didn't copy paste the whole program. I could have been more clear there. I'm preallocating in both. Each uses a different storage scheme. The first uses separate vectors for each x,y,z coordinate generated, the second stores them in an natomsx5 matrix.

[Nathan Phelps]

I have posted the whole program as a comment to my original question if you'd like to take a look. Thanks!

Answer 2

[Roger Stafford]

The progress of your routine will be extremely sensitive to the value you give to ‘rmax’ in the one code or ‘moldia’ in the other. I would hope these are set equal. There will also be a dependence on the initial values in x(1), y(1), and z(1), or in store(1,3:5), and on the values you happen to receive from rand.

I do notice one difference in coding. In the first code you ensure that on the first value with ii = 2 you will definitely go through the while loop, whereas in the second code I see no initial value given to flag to ensure that that happens on the first trip with i = 2. If you happened to start with flag equal to 0, it could have a lasting effect on subsequent processing.

However, as Per states, you should use ‘profile’ to analyze how your code is performing.

An interesting investigation would be to previously store a large number of values from rand in an array and change the two codes to access this array sequentially instead of rand. Either that or use rng to set the same seed in advance each time. In that case if you have set all parameters equal, the results should be identical. If they are not, you have the means of finding out why they are not by checking where store begins to differ from x, y, and z.

[Nathan Phelps]

I'll check those things and report back. In the second script flag is initialized to 1. I really should have posted the whole thing. I'll do that tonight or tomorrow morning.

[Nathan Phelps]

I posted the scripts from the beginning of the code to the trouble making sections as a comment on my original question if you'd like to take a look. Thanks again for your help!