c========================================================================
c program to solve the model nonlinear elliptic equation
c
c f'' + f' * 2/r - f^m = 0
c
c for f(r)
c
c with f'=0 at f=0 and f = 1 at r = r_max
c
c to within a tolerance "tolerance"
c========================================================================
        program model_elliptic

        implicit none
        integer iargc, i4arg
        real*8 r8arg

        integer Nr_max,Nr,l_iter,max_iter,i,j
        real*8 correction_weight,r_iter
        logical trace,vtrace
        parameter (Nr_max=100000,trace=.true.,vtrace=.true.)
        real*8 f(Nr_max)
        real*8 dr,r_max,m,dr2,drsqrd,tolerance,res_nrm,r(Nr_max)

        !---------------------------------------------------------------
        ! variables needed for DGBSV 
        !
        ! N - rank of matrix (Nr here)
        ! NRHS - number of right-hand-sides (1 here)
        ! AB - We are solving the matrix system A.x = B
        !      AB holds the storage for the matrix A(i,j),
        !      which is the Jacobian of our differential operator,
        !      i.e. the dependence of the linearized
        !      equation for the unknown at location r_i=(i-1)*dr
        !      w.r.t the variable at location r_j=(j-1)*dr. 
        !      The mapping to LAPACK storage is AB(KL+KU+1+i-j,j) = A(i,j)
        ! B -  On input, the RHS in the matrix equation, which for
        !      here is minus the residual of the differential operator.
        !      On output, holds the solution x, which is the
        !      correction to f for that step of the Newton iteration.
        ! KL - Number of non-zero elements below the diagional in A
        ! KU - "                           above "
        ! LDAB - as defined below
        ! IPIV - storage needed by the routine
        !
        !see dgbsv.f in the LAPACK distribution for more details
        !---------------------------------------------------------------
        integer N,NRHS,LDAB,INFO,KL,KU,LDB
        parameter (KL=1,KU=2,NRHS=1,LDAB=2*KL+KU+1)
        real*8 AB(2*KL+KU+1,Nr_max),B(Nr_max)
        integer IPIV(Nr_max)

        m=i4arg(1,5)
        r_max=r8arg(2,1.0d0)
        tolerance=r8arg(3,1.0d-8)
        max_iter=i4arg(4,1000)
        correction_weight=r8arg(5,1.0d0)
        Nr=i4arg(6,129)

        write(*,*) 'Input parameters:'
        write(*,*) 'm = ',m
        write(*,*) 'r_max = ',r_max
        write(*,*) 'tolerance = ',tolerance
        write(*,*) 'max_iter = ',max_iter
        write(*,*) 'correction_weight = ',correction_weight
        write(*,*) 'Nr = ',Nr

        LDB=Nr
        N=Nr
        dr=r_max/(Nr-1)
        dr2=2*dr
        drsqrd=dr*dr

        if (Nr.gt.Nr_max) then
           write(*,*) 'Error ... Nr_max too small'
           write(*,*) 'Nr_max=',Nr_max,' Nr=',Nr
           write(*,*) 'increase and recompile'
           stop
        end if

        !--------------------------------------------------------------
        ! Solve the equation via Newton iteration to within
        ! specified tolerance.  
        !--------------------------------------------------------------

        res_nrm=tolerance+1
        l_iter=0

        if (trace) write(*,*) 'Solving Elliptics ...'

        !--------------------------------------------------------------
        ! initial guess : f = 1, and initialize r
        !--------------------------------------------------------------
        do i=1,Nr
           f(i)=1
           r(i)=(i-1)*dr
        end do

        do while(res_nrm.gt.tolerance.and.l_iter.lt.max_iter)
           l_iter=l_iter+1
           res_nrm=0
           if (trace) write(*,*) ' l_iter ',l_iter

           !-----------------------------------------------------------
           ! Clear the storage
           !-----------------------------------------------------------
           do i=1,Nr
              B(i)=0
              do j=i-KL,i+KU
                 AB(KL+KU+1+i-min(max(j,1),Nr),min(max(j,1),Nr))=0
              end do
           end do
           
           !-----------------------------------------------------------
           ! r=0 B.C. f'=0
           !-----------------------------------------------------------
           i=1
           B(i)=-(-3*f(1)+4*f(2)-f(3))
           res_nrm=res_nrm+abs(B(i))
           j=1
           AB(KL+KU+1+i-j,j)=-3
           j=2
           AB(KL+KU+1+i-j,j)=4
           j=3
           AB(KL+KU+1+i-j,j)=-1

           !-----------------------------------------------------------
           ! interior points : f'' + f' * 2/r - f^m = 0
           !-----------------------------------------------------------
           do i=2,Nr-1
              B(i)=-( (f(i+1)-2*f(i)+f(i-1))/drsqrd
     &               +(f(i+1)-f(i-1))*2/r(i)/dr2
     &               - f(i)**m )

              res_nrm=res_nrm + abs(B(i))

              j=i-1
              AB(KL+KU+1+i-j,j)=1/drsqrd-2/r(i)/dr2
              j=i
              AB(KL+KU+1+i-j,j)=-2/drsqrd-m*f(i)**(m-1)
              j=i+1
              AB(KL+KU+1+i-j,j)=1/drsqrd+2/r(i)/dr2
           end do

           !-----------------------------------------------------------
           ! outer boundary: f=1
           !-----------------------------------------------------------
           i=Nr
           B(i)=-(f(i)-1)
           res_nrm=res_nrm + abs(B(i))
           j=Nr
           AB(KL+KU+1+i-j,j)=1

           call DGBSV(N,KL,KU,NRHS,AB,LDAB,IPIV,B,LDB,INFO)

           if (INFO.ne.0) then
              write(*,*) 'Error: elliptics.inc solve for alpha'
              write(*,*) '       DGBSV returned info=',INFO
              stop
           end if

           do i=1,Nr
              f(i)=f(i)+B(i)*correction_weight
           end do

           res_nrm=res_nrm/Nr

           if (trace) write(*,*) '    res_nrm=',res_nrm
           r_iter=l_iter
           if (vtrace) call vsxynt('f_iter',r_iter,r,f,nr)
        end do

        if (res_nrm.gt.tolerance) then
           write(*,*) '***********************************************'
           write(*,*) 'Warning ... failed to solve the equation to '
           write(*,*) tolerance,' res_nrm=',res_nrm
        end if

        call vsxynt('f',0.0d0,r,f,nr)

        end