running out of memory while executing matlab script?

6 views (last 30 days)
asim asrar
asim asrar on 31 Mar 2022
Commented: Riccardo Scorretti on 1 Apr 2022
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(:))

Sign in to answer this question.

Accepted Answer

Riccardo Scorretti
Riccardo Scorretti on 31 Mar 2022
Ran in:
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));
Out of memory.
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
asim asrar on 31 Mar 2022
this is th complete script which is running out of memory. whihc even after making the corrections
as prescribed is running out of memory , i too am using a PC of 128 GB ram
please suggest necessary corrections.
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;
Riccardo Scorretti
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:
  1. the bounds of the three nested for loops
  2. 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); % ***

Sign in to comment.

More Answers (0)

Categories

Find more on Spatial Search in Help Center and File Exchange

Tags

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!