running out of memory while executing matlab script?
Show older comments
clc
clear all;
close all
tic();
um = 1e-6;
c = 3e8; % m/s
dx = 4.8 * um;
dy = 4.8 * um;
dz = 4.8 * um;
dt = 1/4 * dx /c;
lamb = 74.9*um;
k = 2*pi / lamb;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
shape = [250,250,250];
N1=shape(1);
N2=shape(2);
N3=shape(3);
R = lamb*3/dx;
eps=10;
mu=1;
xx = linspace(0,N1-1,N1);
yy = linspace(0,N2-1,N2);
zz = linspace(0,N3-1,N3);
[XX,YY,ZZ] = meshgrid(xx,yy,zz);
center = [125, 125,125];
rr = sqrt((center(2)-XX).^2 + (center(1)-YY).^2 + (center(3)-ZZ).^2);
epsr = ones(shape);
mur = ones(shape);
for ii=1:1:length(rr)
for jj = 1:1:length(rr)
for kk=1:1:length(rr)
if rr(ii,jj,kk) <= R
epsr(ii,jj,kk) = epsr(ii,jj,kk)+(eps-1);
end
end
end
end
eps = epsr(:,:,125);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% omega region configuration
N=94;
x = linspace(0,N-1,N);
y=x;
x=x*dx;
y=y*dy;
[X,Y] = meshgrid(x,y);
X=(X(:))';
Y=(Y(:))';
g=zeros(length(X),length(Y));
for ii=1:length(X)
a=X-X(ii);
g(ii,:)=a;
end
gg=zeros(length(X),length(Y));
for ii=1:length(X)
aa=Y-Y(ii);
gg(ii,:)=aa;
end
r = sqrt(g.^2 + gg.^2);
lamb=409; % wavelength
nb=1.33; % refractive index background
siz=[N^2,N^2]; % size of the region of interest (containing the support of f)
dz=16*lamb/siz(1); % axial discretization step (o have a ROI of 16*lamb)
kdz=2*pi/lamb*nb*dz; % wavenumber
%kg=9.9260
kg=kdz;
G = (exp(1i*(kg*r*nb)))/(4*pi*r);
G1=(1/4)*besselh(0,1,kg*nb*r);
figure(1)
subplot 121
imagesc(abs(G))
subplot 122
imagesc(abs(G1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gamma region configuration
N1=250;
x1 = linspace(0,N1-1,N1);
y1=x1;
x1=x1*dx;
y1=y1*dy;
[X1,Y1] = meshgrid(x1,y1);
X1=(X1(:))';
Y1=(Y1(:))';
g1=zeros(length(X1),length(Y));
for ii=1:length(X)
a1=X1-(X(ii)+78*dx);
g1(:,ii)=a1;
end
gg1=zeros(length(X1),length(Y));
for ii=1:length(Y)
aa1=Y1-(Y(ii)+78*dy);
gg1(:,ii)=aa1;
end
gg_sq1 = g1.^2;
gg1_sq1 =gg1.^2;
r1 = sqrt(gg_sq1 + gg1_sq1);
lamb=409; % wavelength
nb=1.33; % refractive index background
siz=[N1^2,N1^2]; % size of the region of interest (containing the support of f)
dz=16*lamb/siz(1); % axial discretization step (o have a ROI of 16*lamb)
kdz=2*pi/lamb*nb*dz; % wavenumber
%kg=9.9260
kg=kdz;
%H = (exp(1i*(kg*r1*nb)))/(4*pi*r1);
H1=(1/4)*besselh(0,1,kg*nb*r1);
figure(2)
subplot 121
imagesc(abs(H1))
subplot 122
imagesc(abs(H1))
toc();
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% incident field calculaion
%u_in = np.exp(1J * k*(np.sqrt((X_g-(1000/4.8+125)*dx)**2+(Y_g-125*dx)**2))).reshape(250,250)[125-47:125+47,125-47:125+47].flatten()
u_in = exp(1i*k*(sqrt((X1-(1000/4.8+125)*dx).^2+(Y1-125*dx).^2)));
u_in=reshape(u_in,[250,250]);
u_in=u_in(125-46:125+47,125-46:125+47);
%u_in = reshape(250,250)[125-47:125+47,125-47:125+47]
u_in=(u_in(:))
Accepted Answer
Hi Asim,
basically because in your code you use (in the wrong way) length instead of size:
- the function length returns the total number of elements of a vector or matrix,
- the function size returns all or some of the dimensions of a vector or matrix
In particular, you should modify your script like that (pay attention to the triple for loops marked by % ***):
clc
clear all;
close all
tic();
um = 1e-6;
c = 3e8; % m/s
dx = 4.8 * um;
dy = 4.8 * um;
dz = 4.8 * um;
dt = 1/4 * dx /c;
lamb = 74.9*um;
k = 2*pi / lamb;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
shape = [250,250,250];
N1=shape(1);
N2=shape(2);
N3=shape(3);
R = lamb*3/dx;
eps=10;
mu=1;
xx = linspace(0,N1-1,N1);
yy = linspace(0,N2-1,N2);
zz = linspace(0,N3-1,N3);
[XX,YY,ZZ] = meshgrid(xx,yy,zz);
center = [125, 125,125];
rr = sqrt((center(2)-XX).^2 + (center(1)-YY).^2 + (center(3)-ZZ).^2);
epsr = ones(shape);
mur = ones(shape);
for ii=1:size(rr,1) % ***
for jj = 1:size(rr,2) % ***
for kk=1:size(rr,3) % ***
if rr(ii,jj,kk) <= R
epsr(ii,jj,kk) = epsr(ii,jj,kk)+(eps-1);
end
end
end
end
eps = epsr(:,:,125);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% omega region configuration
N=94;
x = linspace(0,N-1,N);
y=x;
x=x*dx;
y=y*dy;
[X,Y] = meshgrid(x,y);
X=(X(:))';
Y=(Y(:))';
g=zeros(length(X),length(Y));
for ii=1:length(X)
a=X-X(ii);
g(ii,:)=a;
end
gg=zeros(length(X),length(Y));
for ii=1:length(X)
aa=Y-Y(ii);
gg(ii,:)=aa;
end
r = sqrt(g.^2 + gg.^2);
lamb=409; % wavelength
nb=1.33; % refractive index background
siz=[N^2,N^2]; % size of the region of interest (containing the support of f)
dz=16*lamb/siz(1); % axial discretization step (o have a ROI of 16*lamb)
kdz=2*pi/lamb*nb*dz; % wavenumber
%kg=9.9260
kg=kdz;
G = (exp(1i*(kg*r*nb)))/(4*pi*r);
G1=(1/4)*besselh(0,1,kg*nb*r);
figure(1)
subplot 121
imagesc(abs(G))
subplot 122
imagesc(abs(G1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gamma region configuration
N1=250;
x1 = linspace(0,N1-1,N1);
y1=x1;
x1=x1*dx;
y1=y1*dy;
[X1,Y1] = meshgrid(x1,y1);
X1=(X1(:))';
Y1=(Y1(:))';
g1=zeros(length(X1),length(Y));
for ii=1:length(X)
a1=X1-(X(ii)+78*dx);
g1(:,ii)=a1;
end
gg1=zeros(length(X1),length(Y));
for ii=1:length(Y)
aa1=Y1-(Y(ii)+78*dy);
gg1(:,ii)=aa1;
end
gg_sq1 = g1.^2;
gg1_sq1 =gg1.^2;
r1 = sqrt(gg_sq1 + gg1_sq1);
lamb=409; % wavelength
nb=1.33; % refractive index background
siz=[N1^2,N1^2]; % size of the region of interest (containing the support of f)
dz=16*lamb/siz(1); % axial discretization step (o have a ROI of 16*lamb)
kdz=2*pi/lamb*nb*dz; % wavenumber
%kg=9.9260
kg=kdz;
%H = (exp(1i*(kg*r1*nb)))/(4*pi*r1);
H1=(1/4)*besselh(0,1,kg*nb*r1);
figure(2)
subplot 121
imagesc(abs(H1))
subplot 122
imagesc(abs(H1))
toc();
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% incident field calculaion
%u_in = np.exp(1J * k*(np.sqrt((X_g-(1000/4.8+125)*dx)**2+(Y_g-125*dx)**2))).reshape(250,250)[125-47:125+47,125-47:125+47].flatten()
u_in = exp(1i*k*(sqrt((X1-(1000/4.8+125)*dx).^2+(Y1-125*dx).^2)));
u_in=reshape(u_in,[250,250]);
u_in=u_in(125-46:125+47,125-46:125+47);
%u_in = reshape(250,250)[125-47:125+47,125-47:125+47]
u_in=(u_in(:))
That's being said, your script requires a lot of memory, notably due to variables G, G1 and H1:
G 8836x8836 1249198336 double complex
G1 8836x8836 1249198336 double complex
H1 62500x8836 8836000000 double complex
On my PC I didn't got an out of memory error (but I have 128 Gb of RAM...) and the script run in a few seconds.
2 Comments
asim asrar
on 31 Mar 2022
Riccardo Scorretti
on 1 Apr 2022
I don't know what to tell. On my PC your program executed in 391 seconds, with a maximum RAM occupation of 71Gb of RAM:
Rather, the problem is that I got NaN (= wrong) values, but this depends on your algorithm / formulation of the problem. More precisely, I obtain these NaN values after executing for the first time line 149:
g =A'*(A*s -u_in);
In order to avoid misunderstanding, I report hereafter the code that I run. Modified line are marked by % ***. The only modification I added are:
- the bounds of the three nested for loops
- the last line: z1 = H1*diag(u)*f; --> z1 = H1*(diag(u)*f);
I hope it helps.
clc
clear all;
close all
tic();
um = 1e-6;
c = 3e8; % m/s
dx = 4.8 * um;
dy = 4.8 * um;
dz = 4.8 * um;
dt = 1/4 * dx /c;
lamb = 74.9*um;
k = 2*pi / lamb;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
shape = [250,250,250];
N1=shape(1);
N2=shape(2);
N3=shape(3);
R = lamb*3/dx;
eps=10;
mu=1;
xx = linspace(0,N1-1,N1);
yy = linspace(0,N2-1,N2);
zz = linspace(0,N3-1,N3);
[XX,YY,ZZ] = meshgrid(xx,yy,zz);
center = [125, 125,125];
rr = sqrt((center(2)-XX).^2 + (center(1)-YY).^2 + (center(3)-ZZ).^2);
epsr = ones(shape);
mur = ones(shape);
for ii=1:1:size(rr,1) % ***
for jj = 1:1:size(rr,2) % ***
for kk=1:1:size(rr,3) % ***
if rr(ii,jj,kk) <= R
epsr(ii,jj,kk) = epsr(ii,jj,kk)+(eps-1);
end
end
end
end
eps = epsr(:,:,125);
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% omega region configuration
N=94;
x = linspace(0,N-1,N);
y=x;
x=x*dx;
y=y*dy;
[X,Y] = meshgrid(x,y);
X=(X(:))';
Y=(Y(:))';
g=zeros(length(X),length(Y));
for ii=1:length(X)
a=X-X(ii);
g(ii,:)=a;
end
gg=zeros(length(X),length(Y));
for ii=1:length(X)
aa=Y-Y(ii);
gg(ii,:)=aa;
end
r = sqrt(g.^2 + gg.^2);
lamb=409; % wavelength
nb=1.33; % refractive index background
siz=[N^2,N^2]; % size of the region of interest (containing the support of f)
dz=16*lamb/siz(1); % axial discretization step (o have a ROI of 16*lamb)
kdz=2*pi/lamb*nb*dz; % wavenumber
%kg=9.9260
kg=kdz;
%G = (exp(1i*(kg*r*nb)))/(4*pi*r);
G1=(1/4)*besselh(0,1,kg*nb*r);
figure(1)
% subplot 121
% imagesc(abs(G))
subplot 122
imagesc(abs(G1))
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gamma region configuration
N1=250;
x1 = linspace(0,N1-1,N1);
y1=x1;
x1=x1*dx;
y1=y1*dy;
[X1,Y1] = meshgrid(x1,y1);
X1=(X1(:))';
Y1=(Y1(:))';
g1=zeros(length(X1),length(Y));
for ii=1:length(X)
a1=X1-(X(ii)+78*dx);
g1(:,ii)=a1;
end
gg1=zeros(length(X1),length(Y));
for ii=1:length(Y)
aa1=Y1-(Y(ii)+78*dy);
gg1(:,ii)=aa1;
end
gg_sq1 = g1.^2;
gg1_sq1 =gg1.^2;
r1 = sqrt(gg_sq1 + gg1_sq1);
lamb=409; % wavelength
nb=1.33; % refractive index background
siz=[N1^2,N1^2]; % size of the region of interest (containing the support of f)
dz=16*lamb/siz(1); % axial discretization step (o have a ROI of 16*lamb)
kdz=2*pi/lamb*nb*dz; % wavenumber
%kg=9.9260
kg=kdz;
%H = (exp(1i*(kg*r1*nb)))/(4*pi*r1);
H1=(1/4)*besselh(0,1,kg*nb*r1);
figure(2)
subplot 121
imagesc(abs(H1))
%
% subplot 122
% imagesc(abs(H1))
% toc();
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% incident field calculaion
%u_in = np.exp(1J * k*(np.sqrt((X_g-(1000/4.8+125)*dx)**2+(Y_g-125*dx)**2))).reshape(250,250)[125-47:125+47,125-47:125+47].flatten()
u_in = exp(1i*k*(sqrt((X1-(1000/4.8+125)*dx).^2+(Y1-125*dx).^2)));
u_in=reshape(u_in,[250,250]);
u_in=u_in(125-46:125+47,125-46:125+47);
%u_in = reshape(250,250)[125-47:125+47,125-47:125+47]
u_in=(u_in(:));
%%
% omega region scatterer
omega = eps(125-46:125+47 , 125-46: 125+47);
omega=(omega(:))';
% omega region function
f = k^2 * diag(omega.^2 -1);
A = eye(length(f)) - G1*f;
%% iterative parameter configuration
delta = 5*10e-7*(norm(u_in,2));
u_prev = u_in;
u_prevprev = u_in;
t_prev =0;
iter =1;
%%
while iter<10
t = sqrt(1 + 4*t_prev^2 )/2;
mu =(1-t_prev)/t;
s = (1-mu)*u_prev + mu*u_prevprev;
g =A'*(A*s -u_in);
gamma = (norm(g,2))^2 / (norm(A*g ,2))^2;
if gamma*(norm(g,2)) <= delta
break
end
u = s - gamma*g;
u_prev = u;
u_prevprev = u_prev;
t_prev = t;
iter =iter+1;
end
X11 = H1*diag(u);
z1 = H1*(diag(u)*f); % ***
More Answers (0)
Categories
Find more on Execution Speed in Help Center and File Exchange
on 31 Mar 2022
on 1 Apr 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
