MATLAB: How can vectorize the below codes

for loopMATLABvectorization

Nx=512;
Ny=Nx
    x0 = Nx/2;
    y0 = Ny/2;
isolve = 2;
tic
    for i=1:512
        for j=1:512
            
            ii=(i-1)*Nx+j;
            
            if(isolve == 2)
                etas(ii,1)=1.0;
                phis(ii,1) =0;
                phis(ii,2)=0;
                etas(ii,2)=0.0;
            else
                etas(i,j,1) =1.0;
                etas(i,j,2) =0.0;
                phis(ii,1) =0;
                phis(ii,2)=0;
            end
            
            xlength =sqrt((i-x0)^2+(j-y0)^2);
            if(xlength <= 14)
                
                if(isolve == 2)
                    etas(ii,1)=0.0;
                    etas(ii,2)=1.0;
                    phis(ii,1) =0;
                    phis(ii,2)=0;
                else
                    phis(ii,1) =0;
                    phis(ii,2)=0;
                end
  
            end %if
        end %j
    end %i

Best Answer

Well, not vectorizing it, but if you preallocate etas and phis with these line:

etas = zeros(262144,2);
phis = zeros(262144,2);

you will speed up the code by a huge factor. (For me, it was a factor of over 1000.)

Related Solutions

MATLAB: How can calcuate the lines of code efficiently

Yes, you can if you avoid using a loop and vectorize your code - see here:

tic
%%%code

nx=512;
ny=nx;
dx=1;
dy=dx;
format long;
nxny=nx*ny;
r=zeros(1,nx);
r(1:2)=[2,-1];
T=toeplitz(r);
E=speye(nx);
grad=-(kron(T,E)+kron(E,T));%Tensor


%-- for periodic boundaries


ii=(0:nx-1).*nx+1;
jj = ii+nx-1;
grad(ii,jj)=1.0;
grad(jj,ii)=1.0;
kk=(nxny-nx+1):nxny;
grad(1:nx,kk)=1.0;
grad(kk,1:nx)=1.0;
grad = grad /(dx*dy);
toc
Elapsed time is 0.145899 seconds.

Your code:

tic
%%%code
nx=512;
ny=nx;
dx=1;
dy=dx;
format long;
nxny=nx*ny;
r=zeros(1,nx);
r(1:2)=[2,-1];
T=toeplitz(r);
E=speye(nx);
grad=-(kron(T,E)+kron(E,T));%Tensor
%-- for periodic boundaries
for i=1:nx
    ii=(0:nx-1).*nx+1;
    jj=ii+nx-1;
    grad(ii,jj)=1.0;
    grad(jj,ii)=1.0;
    kk=nxny-nx+i;
    grad(i,kk)=1.0;
    grad(kk,i)=1.0;
end 
toc
Elapsed time is 11.333450 seconds.

Are the results correct? - to check this we call one of the results grad1 and subtract grad1 from grad:

nx=512;
ny=nx;
dx=1;
dy=dx;
format long;
nxny=nx*ny;
r=zeros(1,nx);
r(1:2)=[2,-1];
T=toeplitz(r);
E=speye(nx);
grad=-(kron(T,E)+kron(E,T));%Tensor
%-- for periodic boundaries
for i=1:nx
    ii=(0:nx-1).*nx+1;
    jj=ii+nx-1;
    grad(ii,jj)=1.0;
    grad(jj,ii)=1.0;
    kk=nxny-nx+i;
    grad(i,kk)=1.0;
    grad(kk,i)=1.0;
end 
ii=(0:nx-1).*nx+1;
jj = ii+nx-1;
grad(ii,jj)=1.0;
grad(jj,ii)=1.0;
kk=(nxny-nx+1):nxny;
grad(1:nx,kk)=1.0;
grad(kk,1:nx)=1.0;
grad1 = grad /(dx*dy);
expect_all_zeros = grad1-grad
expect_all_zeros =
   All zero sparse: 262144×262144

Appears to work fine...

MATLAB: Is the contour plot flipped i.e. the value that should be at the bottom right is plotted on the top left

YOu need not to use loops....it is very slow. Check the below code:

Nx = 100 ; Ny = 100 ; 
x = linspace(0,1,Nx) ; 
y = linspace(0,1,Ny) ; 
[X,Y] = meshgrid(x,y) ; 
phi = 500*exp(-50*((1-X).^2+Y.^2))+100*X.*(1-Y) ; 
contour(X,Y,phi)

Best Answer

Related Solutions

MATLAB: How can calcuate the lines of code efficiently

MATLAB: Is the contour plot flipped i.e. the value that should be at the bottom right is plotted on the top left

Related Question