running out of memory while executing matlab script?

1 次查看（过去 30 天）
asim asrar
asim asrar 2022-3-31
评论： Riccardo Scorretti 2022-4-1
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(:))

请先登录，再回答此问题。

采纳的回答

Riccardo Scorretti
Riccardo Scorretti 2022-3-31
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 个评论
asim asrar
asim asrar 2022-3-31
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 2022-4-1
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); % ***

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Classification Learner App 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by