How to convert MATLAB code to C++ code?

How do I convert my MATLAB code to a C++ code?

Answer 1

[Ameer Hamza]

You will need MATLAB Coder to perform the conversion. If it is already installed, you can use MATLAB coder app under Apps tab in MATLAB. It will guide you through steps to generate C or C++ code.

[PatrizioGraziosi]

Hi Ameer,

thanks for the answer.

There are mainly two reasons, as far as I see, to port a code from MATLAB to C++:

1) to speed it up by using it inside MATLAB environment,

2) to use it outside MATLAB, without having MATLAB installed.

MATLAB Coder should work for the first, does it work also for the second?

I don' t have any experience in this task but I'm going to face this challenge with a big code (tens of thousands of MATLAB code lines in many files) in the forthcoming months.

Thanks in advance

Patrizio

[Ameer Hamza]

Yes, you can use the generated C/C++ code in external projects without the need to have MATLAB installed. See here: https://www.mathworks.com/help/coder/deployment.html?s_tid=CRUX_lftnav, also check this example: https://www.mathworks.com/help/coder/ug/generate-and-modify-an-example-cc-main-function.html.

[PatrizioGraziosi]

Hi Ameer,

thank you!

the C code generated in this way can be used in any platform, while the mex files called in matlab are sensitive to the operating system, isn't it?

Thank

Patrizio

[Ameer Hamza]

Yes, mex files are platform specific. I think that the generated C code should work on other platforms too, but I am not entirely sure.

Answer 2

[Raghu Boggavarapu]

Use MATLAB Coder to convert your MATLAB Code into C/C++

[Harshal Suresh]

Does it come free with the standard student license?

[Walter Roberson]

@Harshal Suresh

MATLAB Coder is an optional extra-cost product for the Student license.

Answer 3

[Ahmad Nadeem]

I want to Change MATLAB code to C++ can anyone tell me how can I convert this? I am working on phase field simulation

[Raghu Boggavarapu]

As mentioned above, please use MATLAB Coder - MATLAB (mathworks.com) for the purpose.

[Ankit Singh Jaswal]

%MAC algorithm based code developed by Rohit Kanchi
%Contact: [phone]
%All units are in SI System kg,m,s,K,Pa
clc
clear
close all
savepath='D:\CFD_MEE4006\Project\2D MAC Algorithm FDM Code\Simulation_Files\';
fprintf("This is an FDM MAC algorithm based code to predict 2-D unsteady incompressible viscous flows. \n\n");
fprintf("Please make sure to enter the input parameters correctly");
nx=100;
ny=40;
deltax=1;
deltay=1;
a=deltax;
b=deltay;
deltaT=0.001;
timesteps=1000;
mu=10^(-3);
rho=1000;
%Grid Point Generator
[Pgrid,Ugrid,Vgrid]=gridpointgenerator(nx,ny,deltax,deltay);

%Boundary Conditions:
%Top: Slip Wall [Axis of rotation]
%Bottom: No Slip Wall, Adiabatic
%Left: Velocity Inlet
Inletvel=0.1;
%Right: Pressure Boundary Condition
ExitP=101325;

%Grid Points Initialization at t=0s
%Pressure Grid Points:
l=length(Pgrid(:,1));
for i=1:l
    Pgrid(i,3)=ExitP;
end

%u vel grid points:
l=length(Ugrid(:,1));
for i=1:l
    if Ugrid(i,1)==(-(deltax/2))
        Ugrid(i,3)=Inletvel;
    else
        Ugrid(i,3)=0;
    end
end

%v vel grid points:
l=length(Vgrid(:,1));
for i=1:l
    Vgrid(i,3)=0;
end

%%%% START MAC ALGORITHM
utild=Ugrid+0; %Initializing Provisional Velocity matrices
vtild=Vgrid+0;
Pprime=Pgrid+0;
Ugridtemp=Ugrid+0;
Vgridtemp=Vgrid+0;

for timestep=1:timesteps
    
    
    
    flag=0; %This var, if zero, means residual are not sufficiently low.
    %When the residual comes down in the while loop below, a break-
    %-statement is encountered and the soln proceeds to next time-
    %-step after updating values
    iter=1;
    while flag==0
        %Calculate provisional velocities:
        
        %Calculate U~
        for i=1:length(Ugrid(:,1))
            
            if Ugrid(i,1)==(-(deltax/2)) %Left Boundary
                continue
            elseif Ugrid(i,1)==((nx*deltax)-(deltax/2)) %Right Boundary
                
                continue
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            elseif Ugrid(i,2)==0 && Ugrid(i,1)~=((nx*deltax)-(deltax/2)) %Bottom Boundary
                u1=Ugrid(i,3); %i,j
                u2=Ugrid(i+1,3); %i+1,j
                u3=Ugrid(i-1,3); %i-1,j
                
                %i,j+1
                for k=i:length(Ugrid(:,1))
                    if(Ugrid(k,1))==(Ugrid(i,1))
                        u4=Ugrid(k,3);
                        break
                    end
                end
                %i,j-1
                u5=0-u1;
                %v Velocity points for averaging
                for k=1:length(Vgrid(:,1))
                    if Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)
                        vforavg1=Vgrid(k,3); %i,j
                    elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)
                        vforavg2=Vgrid(k,3); %i+1,j
                    end
                end
                vforavg3=0;
                vforavg4=0;
                %Calculate Pressures:
                for k=1:length(Pgrid(:,1))
                    if Pgrid(k,1)==Ugrid(i,1)+(deltax/2) && Pgrid(k,2)==Ugrid(i,2)
                        P1up=Pgrid(k,3); %P i+1,j
                    end
                    if Pgrid(k,1)==Ugrid(i,1)-(deltax/2) && Pgrid(k,2)==Ugrid(i,2)
                        P2up=Pgrid(k,3); %P i,j
                    end
                end
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            elseif Ugrid(i,2)==((ny*deltay)-deltay) && Ugrid(i,1)~=((nx*deltax)-(deltax/2)) %Top Boundary
                
                u1=Ugrid(i,3); %i,j
                u2=Ugrid(i+1,3); %i+1,j
                u3=Ugrid(i-1,3); %i-1,j
                
                %i,j+1
                u4=u1+0;
                %i,j-1
                for k=i:-1:1
                    if(Ugrid(k,1))==(Ugrid(i,1))
                        u5=Ugrid(k,3);
                        break
                    end
                end
                %v Velocity points for averaging
                for k=1:length(Vgrid(:,1))
                    if Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)
                        vforavg1=Vgrid(k,3); %i,j
                    elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)
                        vforavg2=Vgrid(k,3); %i+1,j
                    elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)
                        vforavg3=Vgrid(k,3); %i+1,j-1
                    elseif Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)
                        vforavg4=Vgrid(k,3); %i-1,j-1
                    end
                end
                
                %Calculate Pressures:
                for k=1:length(Pgrid(:,1))
                    if Pgrid(k,1)==Ugrid(i,1)+(deltax/2) && Pgrid(k,2)==Ugrid(i,2)
                        P1up=Pgrid(k,3); %P i+1,j
                    end
                    if Pgrid(k,1)==Ugrid(i,1)-(deltax/2) && Pgrid(k,2)==Ugrid(i,2)
                        P2up=Pgrid(k,3); %P i,j
                    end
                end
                
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            else %Interior Grid Points
                
                u1=Ugrid(i,3); %i,j
                u2=Ugrid(i+1,3); %i+1,j
                u3=Ugrid(i-1,3); %i-1,j
                
                %i,j+1
                for k=i:length(Ugrid(:,1))
                    if(Ugrid(k,1))==(Ugrid(i,1))
                        u4=Ugrid(k,3);
                        break
                    end
                end
                %i,j-1
                for k=i:-1:1
                    if(Ugrid(k,1))==(Ugrid(i,1))
                        u5=Ugrid(k,3);
                        break
                    end
                end
                %v Velocity points for averaging
                for k=1:length(Vgrid(:,1))
                    if Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)
                        vforavg1=Vgrid(k,3); %i,j
                    elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)
                        vforavg2=Vgrid(k,3); %i+1,j
                    elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)
                        vforavg3=Vgrid(k,3); %i+1,j-1
                    elseif Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)
                        vforavg4=Vgrid(k,3); %i-1,j-1
                    end
                end
                
                %Calculate Pressures:
                for k=1:length(Pgrid(:,1))
                    if Pgrid(k,1)==Ugrid(i,1)+(deltax/2) && Pgrid(k,2)==Ugrid(i,2)
                        P1up=Pgrid(k,3); %P i+1,j
                    end
                    if Pgrid(k,1)==Ugrid(i,1)-(deltax/2) && Pgrid(k,2)==Ugrid(i,2)
                        P2up=Pgrid(k,3); %P i,j
                    end
                end
            end
            
            delcU= rho*((u1*((u2-u3)/(2*a))) + (((vforavg1+vforavg2+vforavg3+vforavg4)/4)*((u4-u5)/(2*b))));
            deldU= mu*( ((u2-(2*u1)+u3)/(a^2)) + ((u4-(2*u1)+u5)/(b^2)));
            %Calc provisional Up:
            utild(i,3)= u1 + ((deltaT/rho)* (deldU-delcU + ((P2up-P1up)/a)));
        end
        
        %Calculate V~
        for i=1:length(Vgrid(:,1))
            if Vgrid(i,1)==0 || Vgrid(i,1)==((nx-1)*deltax) || Vgrid(i,2)==(0-(deltay/2)) || Vgrid(i,2)==((ny*deltay)-(deltay/2))  %left, right, bottom and top boundaries
                continue
            else %Interior grid points
                v1=Vgrid(1,3); %i,j
                for k=1:length(Vgrid(:,1)) %i,j+1
                    if Vgrid(k,1)==Vgrid(i,1) && Vgrid(k,2)==(Vgrid(i,2)+deltay)
                        v2=Vgrid(k,3);
                        break;
                    end
                end
                for k=1:length(Vgrid(:,1)) %i,j-1
                    if Vgrid(k,1)==Vgrid(i,1) && Vgrid(k,2)==(Vgrid(i,2)-deltay)
                        v3=Vgrid(k,3);
                        break;
                    end
                end
                v4=Vgrid(i+1,3); %i+1,j
                v5=Vgrid(i-1,3); %i-1,j
                %calc u vel avg
                for k=1:length(Ugrid(:,1))
                    if Ugrid(k,1)==Vgrid(i,1)-(deltax/2) && Ugrid(k,2)==Vgrid(i,2)+(deltay/2)
                        uforavg1=Ugrid(k,3); %i,j
                    elseif Ugrid(k,1)==Vgrid(i,1)+(deltax/2) && Ugrid(k,2)==Vgrid(i,2)+(deltay/2)
                        uforavg2=Ugrid(k,3); %i+1,j
                    elseif Ugrid(k,1)==Vgrid(i,1)+(deltax/2) && Ugrid(k,2)==Vgrid(i,2)-(deltay/2)
                        uforavg3=Ugrid(k,3); %i+1,j-1
                    elseif Ugrid(k,1)==Vgrid(i,1)-(deltax/2) && Ugrid(k,2)==Vgrid(i,2)-(deltay/2)
                        uforavg4=Ugrid(k,3); %i-1,j-1
                    end
                end
                %Calc Pressures
                for k=1:length(Pgrid(:,1))
                    if Pgrid(k,1)==Vgrid(i,1) && Pgrid(k,2)==(Vgrid(i,2)-(deltay/2))
                        P1vp=Pgrid(k,3); %P i,j
                    end
                    if Pgrid(k,1)==Vgrid(i,1) && Pgrid(k,2)==(Vgrid(i,2)+(deltay/2))
                        P2vp=Pgrid(k,3); %P i,j+1
                    end
                end
            end
            delcV=rho*( (v1*((v2-v3)/(2*b))) + (((uforavg1+uforavg2+uforavg3+uforavg4)/4)*((v4-v5)/(2*a))) );
            deldV=mu*( ((v4-(2*v1)+v5)/(a^2)) + ((v2-(2*v1)+v3)/(b^2)) );
            vtild(i,3)=v1 + ((deltaT/rho)* (deldV-delcV + ((P1vp-P2vp)/b)));
        end
        
        %Calculate P' field
        for i=1:length(Pgrid(:,1))
            px1=Pgrid(i,1)+(deltax/2);
            px2=Pgrid(i,1)-(deltax/2);
            py1=Pgrid(i,2)+(deltay/2);
            py2=Pgrid(i,2)-(deltay/2);
            
            for k=1:length(utild(:,1))
                if utild(k,1)==px1 && utild(k,2)==Pgrid(i,2)
                    up=utild(k,3);
                end
                if utild(k,1)==px2 && utild(k,2)==Pgrid(i,2)
                    uw=utild(k,3);
                end
            end
            
            for k=1:length(vtild(:,1))
                if vtild(k,1)==Pgrid(i,1) && vtild(k,2)==py1
                    vp=vtild(k,3);
                end
                if vtild(k,1)==Pgrid(i,1) && vtild(k,2)==py2
                    vs=vtild(k,3);
                end
            end
            Pprime(i,3)= -(( ((up-uw)/a) +((vp-vs)/b) )/( ((2*deltaT)/rho)*((1/(a^2)) + (1/(b^2))) ));
        end
        
        %Updating utild and vtild to n+1 temporarily to check for-
        %-convergence
        
        for i=1:length(Ugrid(:,1))
            if Ugrid(i,1)==(-(deltax/2)) || Ugrid(i,1)==((nx*deltax)-(deltax/2)) %Left Boundary and right boundary
                continue
            end
            for k=1:length(Pprime(:,1))
                if Pprime(k,1)==Ugrid(i,1)-((deltax)/2) && Pprime(k,2)==Ugrid(i,2)
                    Ugridtemp(i,3)=utild(i,3) - ( (deltaT/(rho*a))*(Pprime(k+1,3)-Pprime(k,3)));
                end
            end
        end
        %For right boundary we need special treatment for uvel:
        for i=1:length(Ugrid(:,1))
            if Ugrid(i,1)==((nx*deltax)-(deltax/2))
                Ugridtemp(i,3)=Ugridtemp(i-1,3);
            end
        end
        
        for i=1:length(Vgrid(:,1))
            if Vgrid(i,1)==0 || Vgrid(i,1)==((nx-1)*deltax) || Vgrid(i,2)==(0-(deltay/2)) || Vgrid(i,2)==((ny*deltay)-(deltay/2))  %left, right, bottom and top boundaries
                continue
            end
            for k=1:length(Pprime(:,1))
                if Pprime(k,1)==Vgrid(i,1) && Pprime(k,2)==(Vgrid(i,2)-(deltay/2))
                    paa=Pprime(k,3);
                elseif Pprime(k,1)==Vgrid(i,1) && Pprime(k,2)==(Vgrid(i,2)+(deltay/2))
                    pbb=Pprime(k,3);
                end
            end
            Vgridtemp(i,3)=vtild(i,3) - ( (deltaT/(rho*b))*(pbb-paa));
        end
        
        %For right and top boundary we need special treatment for vvel:
        for i=1:length(Vgrid(:,1))
            if Vgrid(i,1)==((nx-1)*deltax) %right
                Vgridtemp(i,3)=Vgrid(i-1,3);
            elseif Vgrid(i,1)==((ny*deltay)-(deltay/2)) %top
                for k=1:length(Vgrid(:,1))
                    if Vgrid(k,1)==Vgrid(i,1) && Vgrid(k,2)==(Vgrid(i,2)-deltay)
                        Vgridtemp(i,3)=Vgrid(k,3);
                    end
                end
            end
        end
        %Pressure correction
        for i=1:length(Pgrid(:,1))
            
            if Pgrid(i,1)==((nx-1)*deltax)
                continue
            else
                Pgrid(i,3)=Pgrid(i,3)+Pprime(i,3);
            end
        end
        %To ensure dp/dy at boundary is zero:
        
        %Check for convergence using temporary values:
        %[Convergence Check: L1 norm]
        sum=0;
        for i=1:length(Pgrid(:,1))
            px1=Pgrid(i,1)+(deltax/2);
            px2=Pgrid(i,1)-(deltax/2);
            py1=Pgrid(i,2)+(deltay/2);
            py2=Pgrid(i,2)-(deltay/2);
            
            for k=1:length(Ugridtemp(:,1))
                if Ugridtemp(k,1)==px1 && Ugridtemp(k,2)==Pgrid(i,2)
                    up=Ugridtemp(k,3);
                end
                if Ugridtemp(k,1)==px2 && Ugridtemp(k,2)==Pgrid(i,2)
                    uw=Ugridtemp(k,3);
                end
            end
            
            for k=1:length(Vgridtemp(:,1))
                if Vgridtemp(k,1)==Pgrid(i,1) && Vgridtemp(k,2)==py1
                    vp=Vgridtemp(k,3);
                end
                if Vgridtemp(k,1)==Pgrid(i,1) && Vgridtemp(k,2)==py2
                    vs=Vgridtemp(k,3);
                end
            end
            residual=( ((up-uw)/a) +((vp-vs)/b) );
            sum=sum+abs(residual);
        end
        
        m=sum/length(Pgrid(:,1));
        disp(m);
        %disp(iter);
        
        if m<=0.001 %Converged
            flag=1;
            Vgrid=Vgridtemp; %Proceeding to n+1
            Ugrid=Ugridtemp;
        else %Not converged
            iter=iter+1;
        end
        
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Proceeding to next time step n -> n+1:
    %The points in the u and v grids are first updated.
    %While updation, points such as the inlet grid points for uvel are-
    %-not touched because they are constant vel points and is a -
    %-dirichlet bc.
    
    %Similarly the Pgrid is also updated and the exit P points are not-
    %-touched because it is a constant pressure Dirichlet BC
    
    
    %Save files at n+1 for the current time-step
    if mod((timestep*deltaT),0.01)==0
        file_name = sprintf('U_vel_timestep_%d.txt',timestep);
        id=fopen([savepath file_name],'w');
        for i=1:length(Ugrid(:,1))
            fprintf(id,'%f %f %f\n',Ugrid(i,1),Ugrid(i,2),Ugrid(i,3));
        end
        fclose(id);
        
        file_name = sprintf('V_vel_timestep_%d.txt',timestep);
        id=fopen([savepath file_name],'w');
        for i=1:length(Vgrid(:,1))
            fprintf(id,'%f %f %f\n',Vgrid(i,1),Vgrid(i,2),Vgrid(i,3));
        end
        fclose(id);
        
        file_name = sprintf('P_timestep_%d.txt',timestep);
        id=fopen([savepath file_name],'w');
        for i=1:length(Pgrid(:,1))
            fprintf(id,'%f %f %f\n',Pgrid(i,1),Pgrid(i,2),Pgrid(i,3));
        end
        fclose(id);
        %Saving complete and proceed to next time-step
    end
end
function [Pgrid,Ugrid,Vgrid] = gridpointgenerator(nx,ny,deltax,deltay)
% Grid point generator
%Row by row, upward
%Pressure collocated grid points
Pgridx(1)=0;
Pgridy(1)=0;
for i=2:nx
    Pgridx(i)=Pgridx(i-1)+deltax;
end

for i=2:ny
    Pgridy(i)=Pgridy(i-1)+deltay;
end
count=1;
for i=1:ny
    for j=1:nx
        Pgridxp(count)=Pgridx(j);
        Pgridyp(count)=Pgridy(i);
        count=count+1;
    end
end

%x vel collocated grid points
nx2=nx+1;
ugridx(1)=-(deltax/2);
ugridy(1)=0;
for i=2:nx2
    ugridx(i)=ugridx(i-1)+deltax;
end
for i=2:ny
    ugridy(i)=ugridy(i-1)+deltay;
end

count=1;
for i=1:ny
    for j=1:nx2
        ugridxp(count)=ugridx(j);
        ugridyp(count)=ugridy(i);
        count=count+1;
    end
end
%y vel collocated grid points
ny2=ny+1;
vgridy(1)=-(deltay/2);
vgridx(1)=0;
for i=2:nx
    vgridx(i)=vgridx(i-1)+deltax;
end
for i=2:ny2
    vgridy(i)=vgridy(i-1)+deltay;
end

count=1;
for i=1:ny2
    for j=1:nx
        vgridxp(count)=vgridx(j);
        vgridyp(count)=vgridy(i);
        count=count+1;
    end
end
Pgridxp=Pgridxp';
Pgridyp=Pgridyp';
ugridxp=ugridxp';
ugridyp=ugridyp';
vgridxp=vgridxp';
vgridyp=vgridyp';
figure
plot(Pgridxp,Pgridyp,'o','MarkerSize',5,'Markeredgecolor','k','markerfacecolor','k');
hold on
plot(ugridxp,ugridyp,'x','MarkerSize',5,'Markeredgecolor','b','markerfacecolor','b');
hold on
plot(vgridxp,vgridyp,'*','MarkerSize',5,'Markeredgecolor','g','markerfacecolor','g');
title('FDM Grid Points')
legend('Pressure Points','u vel points','v vel points');
Pgrid=[Pgridxp,Pgridyp];
Ugrid=[ugridxp,ugridyp];
Vgrid=[vgridxp,vgridyp];
end

[Steven Lord]

If you're trying to ask how to convert that large section of MATLAB code to C code, use MATLAB Coder as @Raghu Boggavarapu said in their answer. If while doing so you encounter difficulties, post as a new question the specific error and/or warning messages you receive and ask for help in that new question.

Answer 4

[Israa Alhuniti]

function [Top_predator_fit,Top_predator_pos,Convergence_curve]=MPA(SearchAgents_no,Max_iter,lb,ub,dim,fobj)





Top_predator_pos=zeros(1,dim);

Top_predator_fit=inf; 



Convergence_curve=zeros(1,Max_iter);

stepsize=zeros(SearchAgents_no,dim);

fitness=inf(SearchAgents_no,1);





Prey=initialization(SearchAgents_no,dim,ub,lb);



Xmin=repmat(ones(1,dim).*lb,SearchAgents_no,1);

Xmax=repmat(ones(1,dim).*ub,SearchAgents_no,1);





Iter=0;

FADs=0.2;

P=0.5;



while Iter<Max_iter    
    
    %------------------- Detecting top predator -----------------    
    
    for i=1:size(Prey,1)  
        
        
        
        Flag4ub=Prey(i,:)>ub;
        
        Flag4lb=Prey(i,:)<lb;    
        
        Prey(i,:)=(Prey(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;                    
        
        
        
        fitness(i,1)=fobj(Prey(i,:));
        
        
        
        if fitness(i,1)<Top_predator_fit 
            
            Top_predator_fit=fitness(i,1); 
            
            Top_predator_pos=Prey(i,:);
            
        end          
        
    end
    
    
    
    %------------------- Marine Memory saving ------------------- 
    
    
    
    if Iter==0
        
        fit_old=fitness;    Prey_old=Prey;
        
    end
    
    
    
    Inx=(fit_old<fitness);
    
    Indx=repmat(Inx,1,dim);
    
    Prey=Indx.*Prey_old+~Indx.*Prey;
    
    fitness=Inx.*fit_old+~Inx.*fitness;
    
    
    
    fit_old=fitness;    Prey_old=Prey;
    
    
    
    %------------------------------------------------------------   
    
    
    
    Elite=repmat(Top_predator_pos,SearchAgents_no,1);  %(Eq. 10) 
    
    CF=(1-Iter/Max_iter)^(2*Iter/Max_iter);
    
    
    
    RL=0.05*levy(SearchAgents_no,dim,1.5);   %Levy random number vector
    
    RB=randn(SearchAgents_no,dim);          %Brownian random number vector
    
    
    
    for i=1:size(Prey,1)
        
        for j=1:size(Prey,2)        
            
            R=rand();
            
            %------------------ Phase 1 (Eq.12) ------------------- 
            
            if Iter<Max_iter/3 
                
                stepsize(i,j)=RB(i,j)*(Elite(i,j)-RB(i,j)*Prey(i,j));                    
                
                Prey(i,j)=Prey(i,j)+P*R*stepsize(i,j); 
                
                
                
                %--------------- Phase 2 (Eqs. 13 & 14)----------------
                
            elseif Iter>Max_iter/3 && Iter<2*Max_iter/3 
                
                
                
                if i>size(Prey,1)/2
                    
                    stepsize(i,j)=RB(i,j)*(RB(i,j)*Elite(i,j)-Prey(i,j));
                    
                    Prey(i,j)=Elite(i,j)+P*CF*stepsize(i,j); 
                    
                else
                    
                    stepsize(i,j)=RL(i,j)*(Elite(i,j)-RL(i,j)*Prey(i,j));                     
                    
                    Prey(i,j)=Prey(i,j)+P*R*stepsize(i,j);  
                    
                end  
                
                
                
            else 
                
                
                
                stepsize(i,j)=RL(i,j)*(RL(i,j)*Elite(i,j)-Prey(i,j)); 
                
                Prey(i,j)=Elite(i,j)+P*CF*stepsize(i,j);  
                
                
                
            end  
            
        end                                         
        
    end         
    
    for i=1:size(Prey,1)  
        
        
        
        Flag4ub=Prey(i,:)>ub;  
        
        Flag4lb=Prey(i,:)<lb;  
        
        Prey(i,:)=(Prey(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
        
        
        
        fitness(i,1)=fobj(Prey(i,:));
        
        
        
        if fitness(i,1)<Top_predator_fit 
            
            Top_predator_fit=fitness(i,1);
            
            Top_predator_pos=Prey(i,:);
            
        end     
        
    end
    
    
    
    
    if Iter==0
        
        fit_old=fitness;    Prey_old=Prey;
        
    end
    
    
    
    Inx=(fit_old<fitness);
    
    Indx=repmat(Inx,1,dim);
    
    Prey=Indx.*Prey_old+~Indx.*Prey;
    
    fitness=Inx.*fit_old+~Inx.*fitness;
    
    
    
    fit_old=fitness;    Prey_old=Prey;
    
    
    
    if rand()<FADs
        
        U=rand(SearchAgents_no,dim)<FADs;                                                                                              
        
        Prey=Prey+CF*((Xmin+rand(SearchAgents_no,dim).*(Xmax-Xmin)).*U);
        
        
        
    else
        
        r=rand();  Rs=size(Prey,1);
        
        stepsize=(FADs*(1-r)+r)*(Prey(randperm(Rs),:)-Prey(randperm(Rs),:));
        
        Prey=Prey+stepsize;
        
    end
    
    
    
    Iter=Iter+1;  
    
    Convergence_curve(Iter)=Top_predator_fit; 
    
    
    
end