Program NEWT_SYS  ! main pt. Newton_Sys
real, allocatable :: x(:), f(:)
!external fcn, jfcn
character(260) fname, nume*100
character*1 w
character*80 text(-2:2)
integer width
character(6) str_n

character(132), allocatable :: f_text(:), x0_text(:)*6
character(132) f_source
character(2) i_char
real(4) norm$
integer status_result


text(-2)=' ** Matrice jacobian singulara'
text(-1)=' ** Divergenta'
text(0) =' * XTOL atins'
text(1) =' ** LNIT depasit'
text(2) =' ** Proces stationar'

print '(a, \)', ' Introdu numarul de ecuatii: '
read *, n

allocate(f_text(n), x0_text(n))
! Initialize with "f(i) = ?":
do i=1,n
iwi =width(i)
write(i_char, '(I<iwi>)') i
f_text(i) ='f('//trim(i_char)//') = ?'
x0_text(i) ='x0('//trim(i_char)//')'
enddo

! Open "f_source" and read:
f_source ='fcn.f90'
call F_Expr_Sys(f_source, n, f_text)

do i=1,n
write(*, '(1x, a)') trim(f_text(i))
enddo

print '(/, a, \)', ' Nume fisier de iesire (fara extensie; <=100 char.): '
read(*,'(a)') nume
fname =trim(nume)//'.rez'
open(1, file= fname, status ='unknown')


allocate (x(n), f(n))

write(1, '(a)') ' Metoda NEWTON - pentru Sisteme'

do i=1,n
if(i==1) write(1,*)

write(1, '(1x, a)') trim(f_text(i))
enddo

deallocate(f_text)

nrad =1

iwn =width(n)
write(str_n, '(I<iwn>)') n


!100 print '(/, a, \)',' Aproximatia initiala x0(1:'//trim(str_n)//'): '
100 print '(/, a, <n-1>(a), a)',' Aproximatia initiala: ', &
( trim(x0_text(i))//', ', i=1,n-1 ), trim(x0_text(n))//':'
read *,x

print '(/,a, \)',' Toleranta (xtol), Numar limita de iteratii (lnit): '
read *, xtol, lnit

print '(/,a, \)',' Cod de tiparire iteratii (0/#0): '
read *, kod

iwnrad =width(nrad); iwlnit =width(lnit);
write(1, '(/, a, I<iwnrad>, a1)') ' Radacina nr.', nrad

write(1,*) 'Aproximatia initiala:', x
!write(1,*) ' lnit, xtol:', lnit, xtol
write(1, '(a, 1pG10.3, a, I<iwlnit>, /)') ' xtol:', xtol, '; lnit: ', lnit

! Check XTOL:
wnorm =norm$(n, x)
xtol_min =ulp(wnorm)
IF(xtol <xtol_min) THEN
write(1,'(/,a, 1pe10.3, a)') 'Atentie: XTOL mai mic decat XTOL_Min =', xtol_min, ' !'
write(*,'(/,a, 1pe10.3, a)') 'Atentie: XTOL mai mic decat XTOL_Min =', xtol_min, ' !'

write(*, '(/,a,\)') 'Vrei sa corectezi XTOL ? (Y,y/char): '
read(*,'(a1)') w
if(w =='Y' .or. w =='y') then
xtol = xtol_min
write(1,'(a,1pe10.3,/)') 'XTOL corectat la:', xtol_min
else
write(1,*)
endif
ENDIF

!call Newton_Sys(fcn, jfcn, n, x, lnit, xtol, f, kod, n_div)
call Newton_Sys(n, x, lnit, xtol, f, kod, n_div)

! Explore divergence:
IF(n_div /=0) THEN
iwn_div =width(n_div)
write(1,'(/, a, I<iwn_div>)') '  * Nr. de divergente locale: ', n_div

ENDIF

write(1,'(/, a)') trim(text(kod))
write(*,'(/, a)') trim(text(kod))

iwlnit =width(lnit)
if(kod ==0) then
write(1,*) ' Solutia X: ', x
write(1,*) ' F(X):      ', f
write(1, '(a, I<iwlnit>)') '  Numar de iteratii = ', lnit

else
write(1, '(a, I<iwlnit>)') '  Iteratia numarul: ', lnit
write(1,*) ' Aproximatia curenta X: ', x
write(1,*) ' F(X): ', f
endif



print'(/, a, \)',' Alta radacina ? (Y, char): '
read'(a)',w

if(w =='Y' .or. w =='y') then
nrad =nrad+1
goto 100
endif

print '(/, a, a, /)', ' * Fisier de rezultate: ', trim(fname)

close(1)

END program NEWT_SYS

