Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
preallocation has no effect - or just slow integration?

Subject: preallocation has no effect - or just slow integration?

From: Reimund Loshka

Date: 30 Dec, 2008 02:07:01

Message: 1 of 5

Hello,

I'm currently working on speeding up a matlab code including a for loop. I know that it is probably most effective to vectorize when using matlab, however this could be a bit complicated in my case. Nonetheless I am using the traditional preallocation of memory for the nested for loop - and this is the strange thing. Using tic/toc on my PC Matlab says, that the code takes about 15 sec to execute, with as well as without preallocation. So I'd like to know, whether in fact the preallocation works perfectly and it is the quad integration that inevitably is that time consuming or whether I can change preallocation to speed up the programme. Here is the code:

clear all

phi = 90; p = phi*pi/180;

t = [0.1:0.1:pi/2];
tmin = 76*pi/180;
tmax = 78*pi/180;

ymax = 2;
ystep = 0.05;
y = [-ymax:ystep:ymax];

xmax = 2;
xstep = 0.05;
x = [-xmax:xstep:xmax];

nx=2*(xmax/xstep)+1;
ny=2*(ymax/ystep)+1;

ER=zeros(nx,ny);
EZ=zeros(nx,ny);
EP=zeros(nx,ny);

for k = 1:nx;
    for m = 1:ny;
        a=y(k).^2 + x(m).^2;
        ER(k,m) = 2*cos(p)*quad(@(t) (cos(t).^(1/2)).*sin(t).*cos(t).*besselj(1,sin(t)*abs(sqrt(a))*2*pi),tmin,tmax);
        EZ(k,m) = 2*i*cos(p)*quad(@(t) (cos(t).^(1/2)).*(sin(t).^2).*besselj(0,sin(t)*abs(sqrt(a))*2*pi),tmin,tmax);
        EP(k,m) = 2*sin(p)*quad(@(t) (cos(t).^(1/2)).*sin(t).*besselj(1,sin(t)*abs(sqrt(a))*2*pi),tmin,tmax);
    end
end

I would be grateful for some comments.

regards

Reimund

Subject: preallocation has no effect - or just slow integration?

From: Roger Stafford

Date: 30 Dec, 2008 02:33:01

Message: 2 of 5

"Reimund Loshka" <spamfilterer@arcor.de> wrote in message <gjbvo5$38t$1@fred.mathworks.com>...
> Hello,
>
> I'm currently working on speeding up a matlab code including a for loop. I know that it is probably most effective to vectorize when using matlab, however this could be a bit complicated in my case. Nonetheless I am using the traditional preallocation of memory for the nested for loop - and this is the strange thing. Using tic/toc on my PC Matlab says, that the code takes about 15 sec to execute, with as well as without preallocation. So I'd like to know, whether in fact the preallocation works perfectly and it is the quad integration that inevitably is that time consuming or whether I can change preallocation to speed up the programme. Here is the code:
>
> clear all
>
> phi = 90; p = phi*pi/180;
>
> t = [0.1:0.1:pi/2];
> tmin = 76*pi/180;
> tmax = 78*pi/180;
>
> ymax = 2;
> ystep = 0.05;
> y = [-ymax:ystep:ymax];
>
> xmax = 2;
> xstep = 0.05;
> x = [-xmax:xstep:xmax];
>
> nx=2*(xmax/xstep)+1;
> ny=2*(ymax/ystep)+1;
>
> ER=zeros(nx,ny);
> EZ=zeros(nx,ny);
> EP=zeros(nx,ny);
>
> for k = 1:nx;
> for m = 1:ny;
> a=y(k).^2 + x(m).^2;
> ER(k,m) = 2*cos(p)*quad(@(t) (cos(t).^(1/2)).*sin(t).*cos(t).*besselj(1,sin(t)*abs(sqrt(a))*2*pi),tmin,tmax);
> EZ(k,m) = 2*i*cos(p)*quad(@(t) (cos(t).^(1/2)).*(sin(t).^2).*besselj(0,sin(t)*abs(sqrt(a))*2*pi),tmin,tmax);
> EP(k,m) = 2*sin(p)*quad(@(t) (cos(t).^(1/2)).*sin(t).*besselj(1,sin(t)*abs(sqrt(a))*2*pi),tmin,tmax);
> end
> end
>
> I would be grateful for some comments.
>
> regards
>
> Reimund

  A couple of questions to start with

1. As you have defined p, the value of cos(p) is zero and this nullifies the entries to ER and EZ

2. The only thing changing in your for-loops is 'a'. Why not just vary 'a' in a single for-loop and make a reasonable approximation to extend it to your square mesh points.

Roger Stafford

Subject: preallocation has no effect - or just slow integration?

From: Roger Stafford

Date: 30 Dec, 2008 03:36:02

Message: 3 of 5

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gjc18t$ahv$1@fred.mathworks.com>...
> "Reimund Loshka" <spamfilterer@arcor.de> wrote in message <gjbvo5$38t$1@fred.mathworks.com>...
> > Hello,
> >
> > I'm currently working on speeding up a matlab code including a for loop. I know that it is probably most effective to vectorize when using matlab, however this could be a bit complicated in my case. Nonetheless I am using the traditional preallocation of memory for the nested for loop - and this is the strange thing. Using tic/toc on my PC Matlab says, that the code takes about 15 sec to execute, with as well as without preallocation. So I'd like to know, whether in fact the preallocation works perfectly and it is the quad integration that inevitably is that time consuming or whether I can change preallocation to speed up the programme. Here is the code:
> >
> > clear all
> >
> > phi = 90; p = phi*pi/180;
> >
> > t = [0.1:0.1:pi/2];
> > tmin = 76*pi/180;
> > tmax = 78*pi/180;
> >
> > ymax = 2;
> > ystep = 0.05;
> > y = [-ymax:ystep:ymax];
> >
> > xmax = 2;
> > xstep = 0.05;
> > x = [-xmax:xstep:xmax];
> >
> > nx=2*(xmax/xstep)+1;
> > ny=2*(ymax/ystep)+1;
> >
> > ER=zeros(nx,ny);
> > EZ=zeros(nx,ny);
> > EP=zeros(nx,ny);
> >
> > for k = 1:nx;
> > for m = 1:ny;
> > a=y(k).^2 + x(m).^2;
> > ER(k,m) = 2*cos(p)*quad(@(t) (cos(t).^(1/2)).*sin(t).*cos(t).*besselj(1,sin(t)*abs(sqrt(a))*2*pi),tmin,tmax);
> > EZ(k,m) = 2*i*cos(p)*quad(@(t) (cos(t).^(1/2)).*(sin(t).^2).*besselj(0,sin(t)*abs(sqrt(a))*2*pi),tmin,tmax);
> > EP(k,m) = 2*sin(p)*quad(@(t) (cos(t).^(1/2)).*sin(t).*besselj(1,sin(t)*abs(sqrt(a))*2*pi),tmin,tmax);
> > end
> > end
> >
> > I would be grateful for some comments.
> >
> > regards
> >
> > Reimund
>
> A couple of questions to start with
>
> 1. As you have defined p, the value of cos(p) is zero and this nullifies the entries to ER and EZ
>
> 2. The only thing changing in your for-loops is 'a'. Why not just vary 'a' in a single for-loop and make a reasonable approximation to extend it to your square mesh points.
>
> Roger Stafford

  As to that second point that I hastily mentioned in my earlier response, what I had in mind is this. Clearly the numerous (3*nx^2=19,683) calls on 'quad' are the major time usage in your code. Allocation delays would pale in comparison.

  For each of the three integrations, the quantity 'a' is the only thing changing from one call to another. You could construct a single for-loop that varied 'a' from 0 in steps of, say, .05 up to 2*sqrt(2). Then for each point (x((m),y(k)) you could approximate the values of ER, EZ, and EP at

 a = sqrt(x(m)^2+y(k)^2)

using interpolation for the nearby values of 'a' that have been already computed.

  This way you have only 3*sqrt(2)*nx (=344) instead of 3*nx^2 calls on 'quad' to make.

Roger Stafford

Subject: preallocation has no effect - or just slow integration?

From: Roger Stafford

Date: 30 Dec, 2008 05:03:03

Message: 4 of 5

> > "Reimund Loshka" <spamfilterer@arcor.de> wrote in message <gjbvo5$38t$1@fred.mathworks.com>...
> > > ......
> > > t = [0.1:0.1:pi/2];
> > > ......
> > > xmax = 2;
> > > xstep = 0.05;
> > > x = [-xmax:xstep:xmax];
> > > nx=2*(xmax/xstep)+1;
> > > ......

  I thought of a couple more items. When you wrote nx=2*(xmax/xstep)+1, it is clear you expect that nx will exactly equal the integer 81 (which in fact it does in this case.) However, you should realize that matlab cannot exactly represent such decimal fractions as .05, .1, etc. and the above kind of calculation might fool you. For example 7/.07 does not yield an exact 100 for matlab. It is a healthier practice to write such things as nx = round(2*(xmax/xstep)+1) and x = linspace(-xmax,xmax,nx) and you will be assured of no unwelcome surprises.

  As a second item I noticed that you wrote t = [0.1:0.1:pi/2]. The question is why. This t will not be used in the quadrature taking place in 'quad'. It uses its own strategically selected variable values for t to suit the needs of its computations. If you want quadrature using the particular values of t you indicated, you should be calling on something like 'trapz' which carries out trapezoidal integration on uniformly spaced variables of integration.

Roger Stafford

Subject: preallocation has no effect - or just slow integration?

From: Reimund Loshka

Date: 30 Dec, 2008 14:33:02

Message: 5 of 5

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gjca27$q2c$1@fred.mathworks.com>...
> > > "Reimund Loshka" <spamfilterer@arcor.de> wrote in message <gjbvo5$38t$1@fred.mathworks.com>...
> > > > ......
> > > > t = [0.1:0.1:pi/2];
> > > > ......
> > > > xmax = 2;
> > > > xstep = 0.05;
> > > > x = [-xmax:xstep:xmax];
> > > > nx=2*(xmax/xstep)+1;
> > > > ......
>
> I thought of a couple more items. When you wrote nx=2*(xmax/xstep)+1, it is clear you expect that nx will exactly equal the integer 81 (which in fact it does in this case.) However, you should realize that matlab cannot exactly represent such decimal fractions as .05, .1, etc. and the above kind of calculation might fool you. For example 7/.07 does not yield an exact 100 for matlab. It is a healthier practice to write such things as nx = round(2*(xmax/xstep)+1) and x = linspace(-xmax,xmax,nx) and you will be assured of no unwelcome surprises.
>
> As a second item I noticed that you wrote t = [0.1:0.1:pi/2]. The question is why. This t will not be used in the quadrature taking place in 'quad'. It uses its own strategically selected variable values for t to suit the needs of its computations. If you want quadrature using the particular values of t you indicated, you should be calling on something like 'trapz' which carries out trapezoidal integration on uniformly spaced variables of integration.
>
> Roger Stafford

Thank's for your advises; I am not an expert in Matlab (basically approached it from my C++ knowledge some time ago) so these are the information I searched for... it seems that I really have to bear in mind, that the calculations and representations are meant to be numerical;

the idea for treating 'a' in the for loop sounds quite reasonable! though I am still a bit undecided regarding using trapz, as the efficieny of th simpson quadrature approach should be far better than trapezoidal integration - so first I will workout whether the quad function in Matlab gives sufficient good results;

regards

RL

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us