! Metoda BROYDEN pt. sistem neliniar.
! Jacobianul inlocuit cu aproxomatia B_k;
! H_k =B_k^(-1)
! H_k+1 actualizat din H_k:
! The "Bad" update (Broyden,2000; Griewank,2012)

subroutine Broyden_2(n, x, xtol, delta, lnit, f, kod, n_div, wnorm)

real(4) x(n), f(n), x0(n), f0(n)
real(4) A(n,n), b(n), w(n)

real(4) H(n,n), H_old(n,n), A0(n,n) !, ID(n,n)
real(4) s(n), y(n)  !;s =dx; y =df;

logical print_it, xtol_met
real(4) norm$
real(4) prag/1E-6/		! pt. elim


print_it = kod /=0

iter=0
xtol_met=.false.

190 format(/, a, t14, a, t<n*16+13>, a)

199 format(i6, t10, <n>G16.7, t<n*16+10>, <n>(1pG16.7), :, <n>(1pG20.7))

	if(print_it)then
	print 190, 'Iteratia', 'X', 'F(X)'
	write (1,190) 'Iteratia', 'X', 'F(X)'
	endif

! INITIALIZEAZA:
! Calculeaza A=F(w0) - jacobianul lui 'f';
! introdu valorile X in X0 :
90 x0 =x

! introdu valorile f(X0) in f0 :
call fcn(n, x0, f0)

! Calculeaza jacobianul lui f in X0:
! Variaza pe rand, cate un x: x(1) = x0(1)+d; x(k#1) =x0(1); etc.
! Calculeaza A pe coloane: 1; 2; etc. 
jcol_do: do j =1,n
!x =x0
x(j) =x0(j) +delta

! Calculeaza f =f(... x0(j-1), x0(j) +delta, x0(j+1)...)
call fcn(n, x, f)

! Coloana "j" din A:
A(1:n,j) =(f(1:n) -f0(1:n))/delta

! restore x(j)
x(j) =x0(j)
enddo jcol_do

! INVERSEAZA pentru H_0 =[A]^(-1)
! ---------------------
! Initializ. pt. ELIM:
nb =1; kodPrint =0;
!**
A0 =A


do icol =1,n
!initializeaza b:
b =0; b(icol) =1;
! rezolva A*x= b:
A=A0
call ELIM (n, nb, A, b, prag, kod, kodPrint)
H(:,icol) =b
end do

!**
!ID =matmul(A0,H)

!do i=1,n
!write(10,*) ID(i,:)
!enddo

iter =0;
x =x0

n_div =0
wnorm_ant =100.

!  Iteratie:
100 iter=iter+1
call fcn(n, x, f)

! Stocheaza f:
f0 =f

! Calculeza x:
! rezolva A*d_iter= -f(w_iter):
w =-matmul(H,f)
x =x +w

	if(print_it)then
	call fcn(n, x, f)
	print 199, iter, x, f
	write(1,199) iter, x, f		!, w
	endif

!Test xtol:
wnorm =norm$(n,w)
if(wnorm <= xtol) then
xtol_met =.true.
kod =0
endif

! Test Divergenta locala:
IF(abs(wnorm) > abs(wnorm_ant)) THEN
!kod =-1 
n_div =n_div+1

if(print_it) then

write(1, '(a)') ' ** Divergenta locala'
write(*, '(a)') ' ** Divergenta locala'
endif

ENDIF

IF(xtol_met) THEN
print*; print*,' Solutia:',x;
if(.not. print_it) call fcn(n, x, f)
print*,' Functia:',f
print*,' Nr. iteratii:',iter

	if(print_it) then
	print*,' Norma corectiei:', wnorm
	write(1,*) ' Norma corectiei:', wnorm
	endif
lnit =iter
return

ELSEIF(iter > lnit)THEN
print*; print*,' LNIT depasit !';
print*,' Ultima aproximatie - X: ',x
if(.not. print_it) call fcn(n, x, f)
print*,' F(X):',f

	if(print_it) then
	print*,' Norma corectiei:', wnorm
	write(1,*) ' Norma corectiei:', wnorm
	endif
kod =1
return

ELSE
!   goto 90
wnorm_ant =wnorm
call fcn(n, x, f)
s =w; y =f -f0;
H_old =H
call Broyden_Update_2(n, H_old, s, y, H)
goto 100
ENDIF

END

