Optimization of matrix calculations for loop & if statements

Hi,

I´m fairly new to matlab (have taken one course) and I now need to make some scripts in a separate, material science, course. I have made this script below which works fine but takes a really long time. I have done what I can/know how to do.

The problem is that the values on row 1 & 100 aswell as col 1 & 100 needs values from the first/last rows/col as the 'if' statements show.

Could someone give me a few pointers on how to improve it for efficiency?

Thanks in advance!

kappa = -9e-16;
dz = 2.7e-8/100

fileID = fopen('xf0.txt','r');
formatSpec = '%f';
A = fscanf(fileID,formatSpec);
dxrow = A(1:2:end,:);
dxtran = reshape(dxrow,[100,100]);
dx = dxtran';

% Row/Col 2,2 -> 99,99 excl col 1 & 100 %
grad = zeros(100,100);
for i = 1:100
    for j = 1:100
        if i == 1 && j == 1
            z = kappa*(0.5*(dx(end,j)/dz)^2 + 0.5*(dx(i,j+1)/dz)^2 + 0.5*(dx(i+1,j)/dz)^2 + 0.5*(dx(i,end)/dz)^2);
        elseif i == 1 && j == 100
            z = kappa*(0.5*(dx(end,j)/dz)^2 + 0.5*(dx(i,1)/dz)^2 + 0.5*(dx(i+1,j)/dz)^2 + 0.5*(dx(i,j-1)/dz)^2);
        elseif i == 100 && j == 1
            z = kappa*(0.5*(dx(i-1,j)/dz)^2 + 0.5*(dx(i,j+1)/dz)^2 + 0.5*(dx(1,j)/dz)^2 + 0.5*(dx(i,end)/dz)^2);    
        elseif i == 100 && j == 100
            z = kappa*(0.5*(dx(i-1,j)/dz)^2 + 0.5*(dx(i,1)/dz)^2 + 0.5*(dx(1,j)/dz)^2 + 0.5*(dx(i,j-1)/dz)^2);
        elseif i-1 < 1
            z = kappa*(0.5*(dx(end,j)/dz)^2 + 0.5*(dx(i,j+1)/dz)^2 + 0.5*(dx(i+1,j)/dz)^2 + 0.5*(dx(i,j-1)/dz)^2);
        elseif j+1 > 100
            z = kappa*(0.5*(dx(i-1,j)/dz)^2 + 0.5*(dx(i,1)/dz)^2 + 0.5*(dx(i+1,j)/dz)^2 + 0.5*(dx(i,j-1)/dz)^2);
        elseif i+1 > 100
            z = kappa*(0.5*(dx(i-1,j)/dz)^2 + 0.5*(dx(i,j+1)/dz)^2 + 0.5*(dx(1,j)/dz)^2 + 0.5*(dx(i,j-1)/dz)^2);
        elseif j-1 < 1
            z = kappa*(0.5*(dx(i-1,j)/dz)^2 + 0.5*(dx(i,j+1)/dz)^2 + 0.5*(dx(i+1,j)/dz)^2 + 0.5*(dx(i,end)/dz)^2);
        else
            z = kappa*(0.5*(dx(i-1,j)/dz)^2 + 0.5*(dx(i,j+1)/dz)^2 + 0.5*(dx(i+1,j)/dz)^2 + 0.5*(dx(i,j-1)/dz)^2);
        end
        grad(i,j) = z;
    end
end

[Dyuman Joshi]

If I am understanding correct, you are trying to get the values of 4 elements around the the element (i,j) and perform the calculation as per the mentioned formula?

[Alexander]

Yes, that is correct. The 4 adjacent values are needed for the calculation of each point. Where there are none, the ones on "opposite" side has to be used.

I also discovered that if the last "else" row is placed at the begining with it's own conditions it runs a bit faster.

Answer 1

[Matt J]

[m,n]=size(dx);

dX=(  dx(   [m,1:m,1]  ,   [n,1:n,1]  )     ).^2;

k=0.25*kappa/dz^2*[0 1 0; 
                   1 0 1;
                   0 1 0];

grad=conv2(dX,k,'valid');

[Alexander]

Thank you! That seems to work like a charm :)

[Alexander]

Hi again,

Maybe you or someone is able to help me with this. It's along the same lines as my last question, about optimization. There is an additional optimization within the same project that absolutly needs doing. it's along line 90-105 but I will provide the whole code so that its easier to understand. Where the problem lies is that the composition needs to be interpolated linearly to find the right value. I do not know much about optimization...

Thankful for any help or pointers on where the right and applicable information kan be found.

clc, clear, clf

kappa = 9e-16;

dz = 2.7e-8/100;

prompt = "Original aton % Chromium? ";

comp_Cr = 0.48; % Komposition

comp_Fe = 1-comp_Cr;

atom_massCr=(8.6341548e-26); % kg

atom_massFe=(9.2732796e-26);+

range=[1 36]; % FYLL I ANTALET FILER

prefix = 'yapfisim';

n = range(2) ;

M = cell(n, 1) ;

for k = 1 : n

M{k} = zeros(100,100) ;

end

%===================================== Gm

Gmtxt = readmatrix('Gm.txt'); %puts the txt file in a matrix

B = sortrows(Gmtxt); %sorts the values based on composition

%plot(B(1:end,1), B(1:end,2)) %plots the Molar-Gibbs energy

Gfe = B(1,2);

Gcr = B(102,2);

x = B(:,1);

x1 = B(1,1);

x2 = B(102,1);

k = (Gcr-Gfe)/(x2-x1);

Gm = [B(1:end,2)]; %Gibbs molar energy values

EG = [];

for i = 1:length(x)

y = Gfe+x(i)*k;

EG = [EG; (Gm(i)-y)];

end

f = [x,EG];

%=====================================

glist = zeros(range(2),100,100);

g = zeros(range(2),1,1);

gradlist = zeros(range(2),100,100);

S = zeros(range(2),1,1);

%tid = zeros(1,range(2))

tid=[]

for i=range(1):range(2)

if i <= 9

filename=strcat(prefix,'_000',num2str(i),'.txt');

fid = fopen(filename,'r');

line_num = 3; % specify the line number you want to read

formatSpec = '%s'; % specify the format specifier for reading a string

C = textscan(fid, formatSpec, 1, 'HeaderLines', line_num-1); % read the desired line

my_number = C{1};

time=my_number{1}

tid = [tid,my_number];

for k=1:4

tline = fgets(fid);

end

formatSpec = '%f';

A = fscanf(fid,formatSpec);

dxrow = A(1:2:end,:);

dxtran = reshape(dxrow,[100,100]);

dx = dxtran';

else i <= 99

filename=strcat(prefix,'_00',num2str(i),'.txt');

fid = fopen(filename,'r');

line_num = 3; % specify the line number you want to read

formatSpec = '%s'; % specify the format specifier for reading a string

C = textscan(fid, formatSpec, 1, 'HeaderLines', line_num-1); % read the desired line

my_number = C{1};

time=my_number{1}

tid = [tid,my_number];

for k=1:4

tline = fgets(fid);

end

formatSpec = '%f';

A = fscanf(fid,formatSpec);

dxrow = A(1:2:end,:);

dxtran = reshape(dxrow,[100,100]);

dx = dxtran';

end

M{i}=M{i}(i,:)+dx;

for j = 1:10000

target = dxrow(j)

so1 = min(x(x>target)); % nearest value after threshold

so2 = max(x(x<target)); % nearest value before threshold

% nearest value regardless of being bigger or smaller than threshold

[~,u] = sort(abs(x-target));

so1 = x(u(1)) % 1st nearest value

intp1 = find(x==so1)

so2 = x(u(2)) % 2nd nearest value

intp2 = find(x==so2)

index = sort([intp1,intp2])

c = polyfit(x(index(1):index(2)), EG(index(1):index(2)), 1)

P = @(x) polyval(c,x)

glist(i,:,:) = P(target);

g(i,:,:) = sum(glist(i,:,:),'all') %[J/mol]

end

[m,n]=size(dx);

dX=( M{i}( [m,1:m,1] , [n,1:n,1] ) ).^2;

k=0.5*kappa/dz^2*[0 1 0;

1 0 1;

0 1 0];

grad=conv2(dX,k,'valid');

gradlist(i,:,:) = grad;

S(i,:,:) = sum(gradlist(i,:,:),'all') %[J/m3]

end

%====================== MOL Cr/Fe

times=[]

for i = 1:length(tid)

times(i) = str2num(cell2mat(tid(i)));

end

volume= (27e-9)*(27e-9)*1; %m3

% MoleCr = (atom_massCr*1000*(10000*comp))*0.019232211646643; % Number of moles Cr

% MoleFe = (atom_massFe*1000*(10000*comp_Fe))*0.01790670606142; % Number of moles Fe

dens_Cr = 7140; % kg/m3

dens_Fe = 7873; % kg/m3

mCr = dens_Cr*volume*comp_Cr; %kg

mFe = dens_Fe*volume*comp_Fe; %kg

nCr = (mCr*1000)/51.99610; %[g]/[g/mol] --> [mol]

nFe = (mFe*1000)/55.8450; %[g]/[g/mol] --> [mol]

Gmjoule = g*(nCr+nFe);

figure(1)

plot(times,Gmjoule)

xlabel('Time [s]'); ylabel('Chemical energy [J]'); title('Chemical energy')

nablax=S*volume;

figure(2)

plot(times,nablax)

xlabel('Time [s]'); ylabel('Gradient energy [J]'); title('Gradient energy')

G = nablax + Gmjoule;

figure(3)

plot(times,G)

xlabel('Time [s]'); ylabel('Gibbs energy [J]'); title('Gibbs energy')

[Matt J]

Since this is a different question, I suggest you post it in its own separate thread, where it will be more visible.