- 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
running out of memory while executing matlab script?
1 view (last 30 days)
Show older comments
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(:))
0 Comments
Accepted Answer
Riccardo Scorretti
on 31 Mar 2022
Hi Asim,
basically because in your code you use (in the wrong way) length instead of size:
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
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/948779/image.png)
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)
See Also
Categories
Find more on Spatial Search in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!