running out of memory while executing matlab script?

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(:))

Answer 1

[Riccardo Scorretti]

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.

[asim asrar]

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]

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);                                        % ***