c **** SCD21:*S*ign-*C*hanging *D*escent, *1*-Hump Solution, *2*-D region
c ****                                                        (Unit Square).
C **** Modified from SCD22.for to seek positive (or negative) solutions.
C ****    
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 One-Sign 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  onto $S$, 
c ****       Sobolev descent on member of $S$.,
c ****       convergence yields $w_1$ (or $w_2$)  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) = u^3 + lambda * 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 ****    * 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**********************************************************************


c **** Number of divisions of Square $\Omega = I \times I$      
      integer n, m
      parameter (m = 50)
      
c **** fixed bifurcation parameter, used in defs of $f$, $F$, and  $f'$
      real lambda 
      common lambda              
      
      real dx, dx2, tol /.001/                  
      
c **** Stepsizes for rAY projecting onto $S$ and sOBOLEV gradient step
      real rstep /.05/, sstep /1/
      
c **** $J(u)$, sup and $L^2$ norms of residuals, $L^2$ norm of $\nabla J(u)$      
      real bju, de1, de2, gnorm, convg

c **** functions in $H=H_0^{1,2}$.  
c **** Iterate, pos/neg parts, $\nabla J(u)$, and a tempORARY variable
      real u, u1, u2, grdj, v
      
c **** counters, maxIMUM number iterations, convergence flag      
      integer i, j, num /0/, max /150/, go /1/ 

      dimension u(0:m,0:m)
      dimension grdj(0:m,0:m),v(0:m,0:m)
      
      n        = m
      dx       = 1.0 / n
      dx2      = dx ** 2
      lambda   = 9
      
c **** WRITE OUTPUT to file
      open(unit=1,file='d22a.txt')
      
c **** initialize u, parse, project parts onto $S$, reassemble in $S_1$
      do i = 1,n-1                                 
         do j = 1,n-1
            u(i,j) = g(float(i)/n,float(j)/n)
         end do        
      end do
      call sobray(u(0,0),n,rstep)

c *** MAIN LOOP ***********************************************
      do while (go .eq. 1)
         num = num + 1 
         
c **** calculate gradient, sSTEP in the negative direction.
c **** (add normalizing code for "farout" initial estimates)
         call gj(u(0,0), grdj(0,0),gnorm, n,dx2,v(0,0))
         call sobstep(u(0,0),grdj(0,0), sstep, convg, sup, n)

c **** Project back onto $S$
         call sobray(u(0,0),n,rstep)

c *****  DIAGNOSTICS and OUTPUT ************************************
                      
c **** calculate $J(u)$ and norms of RESIDUAL
         call ju(u(0,0), bju, n, dx2)
         call test(u(0,0),de1, de2, n,dx2)
                      
c********* beginning of screen/file data out section ***************
         write(*,*) " num, sup, convg: ", num,sup,convg
         write(*,*) " lambda, n ", lambda, n
         write(*,*) "J(u)     = ",bju
         write(*,*) "max/avg  = ",de1, de2
         write(*,*) "gnorm    = ",gnorm

         WRITE (*, 110) ((u(i,j), i = n/10, 9*n/10,n/10),
     +   j = n/10, 9*n/10,n/10)
         
         write(*,*) ' '
         WRITE (*, 110) ((grdj(i,j), i = n/10, 9*n/10,n/10),
     +   j = n/10, 9*n/10,n/10)

         write (*,*) ' '
         write (1,*) ' num, sup, convg: ', num,sup,convg
         write (1,200)
         write (1,*) "J(u) = ",bju          
         write (1,200)
         write (1,*) " lambda, n ", lambda, n
         write (1,200)
         write (1,*) "gnorm    = ",gnorm
         write (1,200)
         write (1,*) "max/avg = ",de1, de2
         write (1,200)

         WRITE (1, 110) ((u(i,j), i = n/10, 9*n/10,n/10),
     +   j = n/10, 9*n/10,n/10)
         write (1,200)
         WRITE (1, 110) ((grdj(i,j), i = n/10, 9*n/10,n/10),
     +   j = n/10, 9*n/10,n/10)
      
110      format (9(f8.5))
200      format (' ',/)     
c********* end of screen/file data out section ********************
         
                        
c **** check for convergence/maxIMUM number of iterations
         if ((num .gt. max) .or. (convg .lt. tol)) then
            go = 0                              
            close(unit=1)
         end if
      end do                            
      pause 'hit any key'

      end                       
c**** END MAIN program *****************************************      
c
c**** nonlinear term function (and derivative/antiderivative)
c**** allow NON-odd functions via diff exponents for pos/neg branches
      real function f(x)
      real x,lambda, p1 /3/, p2 /3/
      common lambda
      
      if (x .gt. 0) then
         f =      x**p1 + lambda*x
      else
         f =  -(-x)**p2 + lambda*x
      end if
            
      end
      real function fp(x)
      real x,lambda, p1 /3/, p2 /3/
      common lambda
      
      if (x .gt. 0) then
         fp =     p1*x**(p1-1) + lambda
      else
         fp =  p2*(-x)**(p2-1) + lambda
      end if

      end
      real function bf(x)
      real x,lambda, p1 /3/, p2 /3/
      common lambda

      if (x .gt. 0) then
         bf =     (1/(p1+1))*x**(p1+1) + .5*lambda*x**2
      else
         bf =  (1/(p2+2))*(-x)**(p2+1) + .5*lambda*x**2
      end if

      end
c**************************
c**** initial guess function, choose 1, comment out others, recompile
      real function g(x,y)

      real x,y, pi /3.14159/, d

c **** This is initial guess:
      g= 10*(sin(pi*x)*sin(pi*y)+sin(pi*x)*sin(pi*y))
      
      end
c************************** 
c**** Solves $-\Delta u = b$ for $u$, standard Num-Anal algorithm.
c**** [Renka] SOR: w = 2/(1+sqr(1-m**2)), m = spect rad Jac It mat
c****               here m = 1-2*sin(pi/(2*n))**2
      subroutine ilap(b,u,n,dx2)
      
      real w /1.73/, total, convg, u, b, max, tol /.000001/, dx2
      dimension u(0:n,0:n),b(0:n,0:n)
      integer i,j,go, n
      
      go = 1
      do while (go .eq. 1)
         max = 0
         do i = 1,n-1
            do j = 1,n-1
               total = u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1) 
               total = (total-b(i,j)*dx2)/4-u(i,j)
               u(i,j) = u(i,j) + w*total
               convg = abs(total) / (1+abs(u(i,j)))
               if (convg .gt. max) then 
                  max = convg
               end if
            end do
         end do
         if (max .lt. tol) then
            go = 0
         end if
      end do
      end
c**************************
c **** Compute $\nabla J(u) = u + (-\Delta)^{-1} f(u)$
      subroutine gj(u,grdj,gnorm, n,dx2,v)
      
      real u, v, grdj, dx2, gnorm
      dimension u(0:n,0:n),v(0:n,0:n),grdj(0:n,0:n)
      integer i,j,n                           
      
      do i = 1,n-1
         do j = 1,n-1
            v(i,j)    =  f(u(i,j))
         end do
      end do                    
      
      call ilap(v(0,0),grdj(0,0),n,dx2)
      
      gnorm = 0
      do i = 1,n-1
         do j = 1,n-1
            grdj(i,j) = grdj(i,j) + u(i,j)
            gnorm = gnorm + grdj(i,j)**2 
         end do
      end do                            
      
c **** compute $L^2$ norm of $\nabla J(u)$      
      gnorm = (gnorm * dx2)**.5
      
      end
c**************************            
c **** ray project in ASCENT direction elements $u$ onto $S$.
      subroutine sobray(u,n,rstep)
      
      real u, rstep, tol /.00005/, gamu, uu
      dimension u(0:n,0:n)
      integer i,j,n,max /200/, go,num

      go = 1
      num = 0
      do while (go .eq. 1)
         num = num + 1
         
c **** recall gamu = $J'(u)(u)$
         gamu = gam(u(0,0), n)
         uu   = dot(u(0,0),u(0,0),n)
         do i = 1,n-1
            do j = 1,n-1

c              rSTEP by projection of $\nabla J(u)$ onto ray
               u(i,j) = u(i,j) + rstep * gamu * u(i,j) / uu

            end do
         end do
         
         if ((num .gt. max) .or. (abs(gamu) .lt. tol)) then
            go = 0
         end if
      end do                    

      end
c**************************                                    
c **** take negative gradient sSTEP
      subroutine sobstep(u,grdj, sstep, convg, sup, n)
      
      real u, grdj, sstep, convg, uold, sup
      dimension u(0:n,0:n), grdj(0:n,0:n)
      integer i,j,n

      convg = 0
      sup   = 0                                   

      do i = 1,n-1
         do j = 1,n-1
            uold = u(i,j)
            
c ****      this is the sSTEP            
            u(i,j) = u(i,j) - sstep * grdj(i,j)
            
            if (abs(u(i,j)) .gt. sup) then
               sup = abs(u(i,j))
            end if
            convg = convg +  (u(i,j) - uold) ** 2
         end do
      end do                                 
      convg = (convg / (n**2)) ** .5
      
      end
c****************************************
c**** compute Convergence Criteria
c *** max = sup norm of PDE(u), avg = $L^2$ norm of PDE(u)
      subroutine test(u,max, avg, n,dx2)
     
      real u, term, dx2, max , avg
      integer i,j,n
      dimension u(0:n,0:n)        
      avg =    0
      max =    0
            
      sum = 0
      do i=1,n-1
         do j=1,n-1                                       
         
c           SECOND DIFFerencing, use grid points.  Improvable?
            term = u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)-4*u(i,j)
            term = (term*n*n + f(u(i,j)))**2
            avg = avg + term
            
            if (term .gt. max) then
               max = term
            end if
         end do
      end do                           
      max = max ** .5
      avg = (avg * dx2)**.5

      end
c****************************************
c**** compute J(u), dx2/dx2 = 1 in term
c**** $J(u) = ||u||^2/2 - \int F(u) dx$
      subroutine ju(u,bju,n,dx2)
     
      real u, dx2, bju, ui
      integer i,j,n
      dimension u(0:n,0:n)
      
      bju = (dot(u(0,0),u(0,0),n) / 2)
      do i=1,n
         do j=1,n

c           average 4 grid corners to do calculation at centers
            ui    =  (u(i,j)+u(i-1,j)+u(i-1,j-1)+u(i,j-1))/4
            bju   =  bju - bf(ui)*dx2

          end do
      end do       

      end
c**** compute gamma = $||u||^2 - \int u f(u) dx$
      real function gam(u,n)
     
      real u,sum
      dimension u(0:n,0:n)
      integer i,j,n
      
      sum = 0
      do i=1,n
         do j=1,n                
            
c           average 4 grid corners to do calculation at centers
            ui    =  (u(i,j)+u(i-1,j)+u(i-1,j-1)+u(i,j-1))/4
            sum   =  sum - ui * f(ui)

         end do
      end do  
      gam = (sum / n ** 2) + dot(u(0,0),u(0,0),n)
      
      end                                        
c COMPUTE $rangle u, v \langle$ = $\int \nabla u \cdot \nabla v dx$
      real function dot(u,v,n)
            
      real u,v,sum, term
      dimension u(0:n,0:n), v(0:n,0:n)
      integer i,j,n
            
      sum = 0
      do i=1,n
         do j =1,n
                        
c           average adjacent grid corners to do differencing at centers
            u1 = (u(i-1,j)+u(i-1,j-1))/2.
            u2 = (u(i,j)+u(i,j-1))/2.
            u3 = (u(i-1,j-1)+u(i,j-1))/2.
            u4 = (u(i-1,j)+u(i,j))/2.
            
            v1 = (v(i-1,j)+v(i-1,j-1))/2.
            v2 = (v(i,j)+v(i,j-1))/2.
            v3 = (v(i-1,j-1)+v(i,j-1))/2.
            v4 = (v(i-1,j)+v(i,j))/2.

c           (u_x v_x + u_y v_y) dx dy
            term = (u2 - u1)*(v2-v1)+(u4 - u3)*(v4-v3)
         
            sum = sum + term
            
         end do  
      end do  
      
      dot = sum
      end