EDIT based on comments below:
I add the mathematical formulation of my problem below. I am trying to solve an equation of the form
$$
\partial_t f(x,y,t)= (\partial^2_x +\partial^2_y) f(x,y,t) \equiv G(x,y,t),
$$
discretizing this equation we have
$$
f^{k+1}_{i,j}= f^k_{i,j} + \Delta t G^k_{i,j}
$$
where $i,j$ refer to discretized spatial coordinates $x,y$ and $k$ corresponds to the iteration step.
However here $\Delta t $ is a fixed step size, I want to use a line search to find an optimal step size. Defining $\Delta t \equiv \alpha_k$, I want to find $\alpha_k $ such that $f^{k+1}(i,j) < f^{k}(i,j)-c\alpha_k G^\top G$ which is a backtracking Armijo line search. So the equation I am trying to solve is :
$$
f^{k+1}_{i,j}= f^k_{i,j}+\alpha^k G^k (i,j)
$$
Below is a back tracking line search algorithm to find $\alpha_k$ but it is not being computed correctly I realize.
I updated my algorithm based on the comment below however it still seems my stepsize at each iteration, $\alpha_k$ is not being updated properly. When I print it out it just prints the initial value I inputted for it. Is this algorithm not updating $\alpha_k$ correctly? I thought the point of the backtracking line search was to find me an optimal value $\alpha_k$ such that I get to the minimum. How can I fix this? thanks!
I am trying to code the backtracking-Armijo line search algorithm on page 10 here https://people.maths.ox.ac.uk/hauser/hauser_lecture2.pdf.
Below is a sample code for a back tracking line search algorithm . I can check that the algorithm is not correct but I am not sure where I am going wrong. A few errors i've realized are possibly the if statement condition and updating alphak at the very end ; lastly I'm not sure what else is wrong (I am unsure how to fix these problems). I tried to follow the algorithm in the book but it is not too detailed. It is very clear there is a problem with the algorithm but I am not sure what.
Note, the search direction I choose given by Pk is in the direction of the negative gradient, which I call -g. Assume below that g(i,j) and fk(i,j) are given at the first iteration, and are 2D arrays since they depend on spatial positions i,j.
integer, parameter :: nx=10,ny=10, k=10
real, dimension(-nx:nx,-ny:ny) :: fk,fk1,g,gt,Pk
integer :: i,j,m
real :: alphak,c,rho !step size at iteration k
c=0.0001
rho=0.5
do m=1,k
alphak = 1.0
Pk(i,j) = -g(i,j) !search direction = -gradient
fk1(i,j) = fk(i,j) + alphak*Pk(i,j)
gt(i,j) = g(j,i) !transpose of g
if (fk1(i,j) > fk(i,j)-c*alphak*gt*g) then
alphak = rho*alphak
do j=-ny+1,ny-1
do i=-nx+1,nx-1
fk1(i,j) = fk1(i,j)-alphak*g(i,j)
end do
end do
end if
print*, "print out alphak for m=', m
print "(//(5(5x,e22.14)))", alphak
end do
fk1
. $\endgroup$