How to reduce time in a matlab code?

function [f , t] = jacobi1(n)
tic;
if n < 4 
    error ('n not in the range')
end
if rem(n,1)~=0
    error ('n has to be a whole mumber')
end
n=n^2;
b=zeros(n,1);
b(1,1)=1;
b(n,1)=1;
A = zeros(n,n);
mymatrix=[-1 0 -1 4 -1 0 -1];
for i=1:3
    A(i,1:i+3)=mymatrix(5-i:7);
    A(n-(3-i),n-(6-i):n)=mymatrix(1:7-i);   
end
for i=4:n-3
    A(i,i-3:i+3) = mymatrix;
end
epsilon = 1e-3;
f = zeros(n,1);
counter = 0;
flag = 0;
L = tril(A,-1);
D = diag(A);
U = triu(A,1);
B = -1./D.*(L+U);
C = (1./D).*b;
while flag == 0 
    counter = counter+1;
    if counter > 10000
        error ('Too many iteration')
    end 
f_n = (B*f) + C;
    if max(abs(f_n-f)/(f_n))<epsilon
        flag = 1;
    else 
        f = f_n;
    end 
end 

t=toc;
end

[Rik]

This time I edited your question for you. Next time, please use the tools explained on this page to make your question more readable.

Did you try running the profiler to see where your code spends the most time?

[Tzach Berlinsky]

no, i dont really know how to do that. can you guide me through please?

[Tzach Berlinsky]

i ran the profiler, and it says that line 3 - if n > 4 takes the most time, is it possible to decrease that time by changing something?

[Rik]

What input did you use when running the profiler? Because I see a completely different profile. For me line 37 (if max(abs(f_n-f)/(f_n))<epsilon) is taking up 25% of the run time.

[Walter Roberson]

It makes sense for abs(f_n-f)/(f_n) to take a long time, as that is a matrix right division, similar to

abs(f_n-f) * pinv(f_n)

By the way, I suspect you should be taking

max(abs((f_n-f)/f_n))

[Tzach Berlinsky]

abs(f_n-f) * pinv(f_n)

yes!! this one reduces it for 8 sec to 0.8 sec.

what exactly in pinv? what does it do?

[Walter Roberson]

https://www.mathworks.com/help/matlab/ref/pinv.html

If you are going to use pinv() then I believe you need

max(abs( (f_n-f) * pinv(f_n) ))

[Tzach Berlinsky]

what is the difference between the / that i put and the *pinv that you put? are the results similar? does it change the outcome of the results?

[Tzach Berlinsky]

becuase i also see that ther is if max(abs(f_n-f) * rank(f_n)) and it runs even faster.

[Walter Roberson]

https://www.mathworks.com/help/matlab/ref/mrdivide.html

x = B/A solves the system of linear equations x*A = B for x. The matrices A and B must contain the same number of columns.

Mathematically speaking, if A were an invertable matrix and you were using calculations that did not suffer from floating point round-off, then x*A = B could be solved by multiplying both sides on the right by inv(A), giving you x*A*inv(A) = B*inv(A) which is x = B*inv(A) .

And generally speaking, when A is not fully invertable, x = B*pinv(A) can produce a solution -- pinv(A) is, in a way, a generalized inverse that extends to singular matrices as well. For invertable matrices, pinv(A) = inv(A) to within round-off error.

When pinv() is used, then

x_pinv = B*pinv(A)
Breconstructed_pinv = x_pinv * A
sserr_pinv = sum((Breconstructed_pinv - B).^2)

then except for round-off error problems, the sum-squared error between B and the reconstructed B will mathematically be as small as possible.

When / is used:

x_mrd = B/A
Breconstructed_mrd = x_mrd * A
sserr_pinv = sum((Breconstructed_mrd - B).^2)

is not mathematically as small as possible.

However, if I recall correctly, B/A has the property that the reconstructed B will be all 0 except for as many leading entries as there are columns in A -- so in your case with A having one column, the reconstructed B would have one leading non-zero entry and then the rest would be 0.

Using rank() in the situation is just plain WRONG for the situation. It is irrelevant and has nothing to do with the mathematical operation you are trying to perform.

What you are trying to calculate in that step is a measure of how bad the current location is. The calculation (f_n-f)/f_n gives a relative slope, derivatives for each of the components. When the ratio is near 0 in all components, then the location is a good one. If the absolute value of any of the components is too large, then you should keep continuing to look for a better position. But you should not be fitting abs(f_n-f) to f_n as that could give you a large negative component where some element of f_n is negative, and max(abs(f_n-f)/f_n) would pay not attention to the large negative value. Instead, as I have been saying, you should be taking max(abs((f_n-f)/f_n))) as that takes the absolute value of each of the components and looks for the maximum of those. Your code would ignore a component that came out as -1000 if the other components were small enough, whereas my suggestion certainly would not.