c **** sub12.for: sublinear, 1-dim (ODE), 2-hump (sign-changing sols)
c **** Existence has not yet been proven!
c **** Don't completely undertand my own algorithm...
c **** Funky ascent/descent projections of u and parts...
c **** Iterates not monotone..oscillate towards solution...
c **** A little unstable?
c **** 
c **** Might NOT be most current/accurate differencing scheme?
c ****                                                        
c ****                                                 (Unit Square).
C **** John M. Neuberger
C ****    * Department of Mathematics and Statistics
C ****    * Mississippi State University
C ****    * Msstate, MS 39762
C ****    
C ****    * email:  neuberge@math.msstate.edu
C ****    * web:    http://www2.msstate.edu/~neuberge
C ****    * FAX:    (601) 325-0005
C ****    * office: (601) 325-7146
C ****    * home:   (601) 338-0007
c ****
c **** SOLVES: $-\Delta u = f(u)$ on $\Omega$, zero-D conditions.
C ****    * Finds Sign Changing Solutions by Sobolev Gradient Descent
C ****      and Projections on to Submanifolds of H=H_0^{1,2}.
c ****                                                                      
C **** Based on alogithm by John M. Neuberger
C ****    * Analysis: "A sign-changing solution to a superlinear 
C ****      Dirichlet problem", with A. Castro and J. Cossio.
C ****      To appear: Rocky Mnt. J. Math., 1995. ([CCN])
C ****    * Numerical Algorithm: "A numerical method for finding 
C ****      sign-changing solutions of superlinear Dirichlet problems".
C ****      To appear: Nonlinear World, 1996. ([NUM]}
C ****    * Look for a related sequence of recent papers by Neuberger,
C ****      Castro, and/or  Cossio.
c ****
c **** LOOP: Project pos/neg parts onto $S$, 
c ****       Sobolev descent on assembled $S_1$ member,
c ****       convergence yields $w_3$ from [CCN] and [NUM].                      
c ****
c **** More current than Jan 1, 1997.
C ****
C **** Can be easily modified to:  
C ****    * Use different nonlinearites (initially f(u) = tanh(u))
C ****      For more, see the 3 function routines f, fp, and bf and the
C ****      assignment of a value to lambda in the initial code fragment.
C ****    * Find one sign solutions...or other solutions which should exist 
C ****      when symmetry is present (we conjecture much...).
c ****      NOTE: trivial sol is ascent + Proj_S.
c ****            one sign sols require descent (No Projection Necessary).
c ****            Witness below the combination requied for sign-changing!
C ****    * Use more divisions (more accuracy)
C ****      and tougher convergence criteria.
C ****    * Use different initial conditions.  See the function routine g.
C ****    * Run on different regions..I have versions for disk and annulus.
C ****    * Can find semipositone and sublinear solutions with slight to
C ****      somewhat drastic modifications.  In general, this involves 
C ****      different combinations of gradient ascent/descent and projections
C ****      onto S and S_1.
C ****
c**********************************************************************
      integer n, m, num /0/, max /100/, go /1/, pskip /1/
      parameter (m = 10000)
      real tol /.0001/
c      real rstep /-.5/, sstep /.15/
c ****
c **** Note that Ray Projections are now DESCENT...J|_S \leq 0.
c ****
      real rstep /-.75/, sstep /.5/
      real u, u1, u2, grdj
      real gnorm, pnorm, sum, sum2, bju
      integer i
      dimension u(0:m), u1(0:m), u2(0:m), grdj(0:m)
      
      n   = m                        
      
      open(unit=1,file='d12.txt')
                  
c **** initialize u

      do i = 1,n-1                                 
         u(i) = g(float(i)/n)
      end do

c **** Project parts onto S...do NOT (?) project whole onto S_1!
      call posneg(u(0),u1(0),u2(0),n)
      call sobray(u1(0),n,rstep)
      call sobray(u2(0),n,rstep)
c      do i=1,n-1
c         u(i) = u1(i) + u2(i)
c      end do                 

                         
c *** MAIN LOOP ***********************************************
      do while (go .eq. 1)
         num = num + 1 

c **** take ASCENT steps on parts individually
         sstep = -sstep
         call gj(u1(0),  grdj(0), n)
         call sobstep(u1(0), grdj(0), sstep, convg, sup, n)
         call gj(u2(0),  grdj(0), n)
         call sobstep(u2(0), grdj(0), sstep, convg, sup, n)
         sstep = -sstep

c **** Reassemble to sign-changing sol...near S_1.         
         do i=1,n-1
            u(i) = u1(i) + u2(i)
         end do                                   
         
c **** Take DESCENT step on whole, project onto S (not S_1), and split
c **** into pos/neg parts for next iteration.
         call gj(u(0),  grdj(0), n)
         call sobstep(u(0), grdj(0), sstep, convg, sup, n)
         call sobray(u(0),n,rstep)
         call posneg(u(0),u1(0),u2(0),n)
      
c **** ALGORITHM RECAP:  See J(r(t)) graph in [Semipos], [CCN],..........      
c **** It's a tug-of-war...the parts on the end go up...the sign-changing
c **** in the middle goes down...minimizes while remaining...????????????
c ****
c*********** start diagnostics and output *************************


         sum = 0
         sum2 = 0
         do i=1,n
            term = (u(i)-u(i-1))**2        
            ui   = (u(i)+u(i-1))/2
            sum  = sum + term * n / 2 - bf(ui) / n
            sum2 = sum2 + grdj(i)**2
         end do  
         bju   = sum                  
         gnorm = (sum2 / n) ** .5
         
         write(*,*) "J(u) = ",bju       
         write(*,*) "|g|  = ",gnorm
         write(1,*) "J(u) = ",bju       
         write(1,*) "|g|  = ",gnorm
         write(1,*) 'num = ',num,'    sup = ',sup, '    convg = ',convg
         write(*,*) ' num, sup, convg: ', num,sup,convg
         
         WRITE (*, 110) (u(i), i = n/10, 9*n/10,n/10)
         write(*,*)
         WRITE (*, 110) (grdj(i), i = n/10, 9*n/10,n/10)

         WRITE (1, 110) (u(i), i = n/10, 9*n/10,n/10)
         write (1, 200)
         WRITE (1, 110) (grdj(i), i = n/10, 9*n/10,n/10)

110      format (9(f8.5))
200      format (' ',/)     

c*********** end diagnostics and output *************************

         if ((num .gt. max) .or. (10*convg .lt. tol)) then
            go = 0
            close(unit=1)
         end if
         
      end do                            
      pause 'hit any key'
      
      end                       

c**** END MAIN program *****************************************      

      subroutine posneg(u,u1,u2,n)
            
      real u,u1,u2
      dimension u(0:n),u1(0:n),u2(0:n)
      integer i, n
      
      do i = 1,n-1   
         if (u(i) .gt. 0) then
            u1(i) = u(i)
            u2(i) = 0
         else
            u2(i) = u(i)
            u1(i) = 0                  
         end if
      end do      
           
      end                                       
c**************************  
c**** nonlinear term function (and derivative/antiderivative)
C****
c**** SUBLINEAR: f(u) = tanh(u)
c****
      real function f(x)
      real x,lambda /200.0/

c tanh(x)      
         f = lambda*(exp(x)-exp(-x))/(exp(x)+exp(-x))
            
      end
      real function fp(x)
      real x,lambda /200.0/           
      
c (sech(x))^2
         fp = 4*lambda/(exp(x)+exp(-x))**2
            
      end
      real function bf(x)
      real x,lambda /200.0/      
      
c ln(cosh(x))
         bf = lambda*log((exp(x)+exp(-x))/2)
            
      end
c**************************
c**** initial guess function 
      real function g(x)

      real x, pi /3.14159/

      g = sin(2*pi*x)
            
c      if (x .le. .66666) then
c         g = g + 3 * sin(3*pi*x/2)
c      else  
c         g = g - 3 * sin(3*pi*x)
c      end if
           
      end
c
c**************************
c** USE separate compiled versions of dot, gam, fixed integration scheme **
c****************************************
      subroutine ilap(y,n)

      dimension y(0:n)
      integer   i,n

c      inverse laplacian via Gaussian Elimination: solves -lap(u)=f

      do i = 1, n - 1
         y(i) = y(i) / (n**2)
      end do
      do i = 2, n - 1
         y(i) = y(i) + y(i - 1) * (i - 1) / i
      end do
      y(n) = y(n) * n / (n + 1)
      do i = n - 1, 1, -1
         y(i) = (y(i + 1) + y(i)) * i / (i + 1)
      end do                                                       
      end

c******************************************************************
c NOT using SOR this version...
c**************************

c****  grad(gammma(u))
      subroutine gg(u,ggu,n)
      
      real u, ggu
      dimension u(0:n),ggu(0:n)
      integer i,n                           
      
      do i = 1,n-1
         ggu(i) = -(u(i)*fp(u(i))+f(u(i)))
      end do                    
      
      call ilap(ggu(0),n)
      
      do i = 1,n-1                 
         ggu(i) =  2*u(i) + ggu(i)
      end do     
      
      end
c**************************
      subroutine gj(u,grdj,n)
      
      real u, grdj
      dimension u(0:n),grdj(0:n)
      integer i,n                           
      
      do i = 1,n-1
         grdj(i) =   -f(u(i))
      end do                    
      
      call ilap(grdj(0), n)

      do i = 1,n-1
         grdj(i) = u(i) + grdj(i)
      end do     
      
      end
c**************************
      subroutine sobray(u,n,rstep)
      
      real u, rstep, tol /.00005/, gamu, uu
      dimension u(0:n)
      integer i,n,max /50/, go,num

      go = 1
      num = 0
      do while (go .eq. 1)
         num = num + 1
         
         gamu = gam(u(0), n)
         uu   = dot(u(0),u(0),n)
         do i = 1,n-1
            u(i) = u(i) + rstep * gamu * u(i) / uu
         end do
         
         if ((num .gt. max) .or. (abs(gamu) .lt. tol)) then
            go = 0
         end if
      end do                    

      write(*,*) 'sobray steps:',num

      end
c**************************
      subroutine sobstep(u,grd, sstep, convg, sup, n)
      
      real u, grdj, sstep, convg, uold, sup
      dimension u(0:n), grd(0:n)
      integer i,n

      convg = 0
      sup   = 0                                   

      do i = 1,n-1
         uold = u(i)

         u(i) = u(i) - sstep * grd(i)

         if (abs(u(i)) .gt. sup) then
            sup = abs(u(i))
         end if
         convg = convg +  (u(i) - uold) ** 2
      end do     

      convg = (convg ** .5)/ n

      end
c*********************************************************
c**** compute gamma
      real function gam(u,n)
     
      real u,sum,sum2, term
      dimension u(0:n)
      integer i,n
      
      sum = 0
      do i=1,n

         ui = (u(i)+u(i-1))/2
         sum = sum - ui * f(ui) / n         

      end do  
      sum = sum + dot(u(0),u(0),n)
      gam = sum                

      end                                        
c*********************************************************
      real function dot(u,v,n)

      real u,v,sum, term
      dimension u(0:n), v(0:n)
      integer i,n
            
      sum = 0
      do i=1,n
         term = (v(i)-v(i-1))*(u(i)-u(i-1))
         sum = sum + term * n
      end do  
      
      dot = sum
      end
c******************** END OF FILE *************************