!
!     main (nonlineq)
!     Sept 24, 2016
!
!      implicit real *8 (a-h,o-z)

	use msflib
	logical chamada

      parameter(neqmax=15)   
	dimension u(neqmax),g(neqmax),h(neqmax)
	common /pert/ epsilon
	common /method/ iflag
	common /parpicard/ alfa
!
!     input file
!
      open (unit=1,file='inp1.txt',status='old')
!
!     output file
!
      open (unit=2,file='out1.txt',status='unknown') 

	  open (10,file='dados1.txt')
	  open (11,file='dados2.txt')

!****************************************
!
!	Decide what method to use	(iflag = 1 - exact Newton; 2 - approximate Newton; 3 - Picard)
!
	read(1,*)iflag
      write(*,*)'iflag=',iflag
	read(1,*)alfa
      write(*,*)'alfa=',alfa

      read (1,*)maxit,neq,tol1,tol2,epsilon
      write(*,*)'maxit=',maxit,'neq=',neq
	write(*,*)'tol1=',tol1,'tol2=',tol2,'epsilon=',epsilon

!
!	dados para o grafico das duas equacoes
!
	xmax=1.
	xmin=0.
	dx=(xmax-xmin)/100
	xx=xmin
	do i=1,101
		y1=cos(xx)
		write(10,*)xx,y1
	!	write(*,*)xx,y1
		y2=asin(xx)
		write(11,*)xx,y2
	!	write(*,*)xx,y2
		xx=xx+dx
	enddo

	close(10)
	close(11)

	chamada=systemqq('wgnuplot dados.gnu')

!	stop

!	read initial guesses

	read(1,*)(u(i),i=1,neq)
	write(*,*)'the initial values are u(i):'
	write(*,*)(i,u(i),i=1,neq)
 !     call fcn(neq,u,g)
!	  gnorm=rnorm2(neq,g,neqmax) 
!	  write(2,*)'--------------------------------------------------- i = ',i
 !     write(*,*)'gnorm=',gnorm
!	  write(2,*)'gnorm=',gnorm
!	stop


!******************************
!
!     solution of nonlinear equations
!
	if(iflag.eq.1.or.iflag.eq.2) then
      call mnewt(maxit,u,neq,tol1,tol2)
	endif
	if(iflag.eq.3) then

	do i=1,maxit
      call fcn(neq,u,g)
	  gnorm=rnorm2(neq,g,neqmax) 
      write(*,*)'gnorm=',gnorm
	  write(2,*)'gnorm=',gnorm
	  if(gnorm.le.tol1) then
         write(*,*)'converged in ',i-1,' iterations'
		 write(2,*)'converged in ',i-1,' iterations'
         stop
      endif
		write(*,*) 'Picard method'
		write(2,*) 'Picard method'
      write(*,*)'--------------------------------------------------- i = ',i
	  write(2,*)'--------------------------------------------------- i = ',i
      call fcn2(neq,u,h)
	  do j=1,neq
		u(j)=alfa*u(j)+(1.-alfa)*h(j)
        write(*,*)'u(',j,')=',u(j)
		write(2,*)'u(',j,')=',u(j)
	  enddo	
	enddo	
	endif


      stop
      end
!      
!----------------------------------------------------------------
      function rnorm2(n,x,nd)
!
!	compute euclidean norm of a vector
!
!      implicit real *8 (a-h,o-z)
      dimension x(nd)
      sum=0.d0
      do i=1,n
      	sum=sum+x(i)*x(i)
      enddo
      aux=sqrt(sum)
      rnorm2=aux
      return
      end
!------------------------------------------------------------
      function rninf(n,x,nd)
!
!	compute infinity norm of a vector
!
!      implicit real *8 (a-h,o-z)
      dimension x(nd)
      rmax=0.d0
      do i=1,n
      	if(abs(x(i)).gt.rmax) then
           rmax=abs(x(i))
        endif
      enddo
      rninf=rmax
      return
      end
!------------------------------------------------------------------ 
      subroutine usrfun(u,n,np,g,a)
!
!	compute the values of g (through function fcn) 
!	for a given u and the jacobean
!
!      implicit real *8 (a-h,o-z)
      parameter (np1=15)
      dimension u(n),g(np),a(np,np),ge(np1)
      common /pert/ epsilon
      call fcn(n,u,g)
!      do it=1,n
!         write(*,*)'g(',it,')=',g(it)
!      enddo
!
!     build the jacobean numerically
!
!      write(*,*)'epsilon=',epsilon
      do j=1,n
         u(j)=u(j)+epsilon
         call fcn(n,u,ge)
         u(j)=u(j)-epsilon
         do i=1,n
            a(i,j)=(ge(i)-g(i))/epsilon
!            write(*,*)'ge(',i,')=',ge(i)
         enddo
      enddo
      return
      end
!234567890123456789012345678901234567890123456789012345678901234567890 
!--------------------------------------------
