- 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?
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
asim asrar
le 31 Mar 2022
Commenté : Riccardo Scorretti
le 1 Avr 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 commentaires
Réponse acceptée
Riccardo Scorretti
le 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 commentaires
Riccardo Scorretti
le 1 Avr 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); % ***
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Spatial Search dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!