\documentstyle{article}
\setlength{\headheight}{0.0in}
\setlength{\textheight}{8.75in}
\setlength{\textwidth}{6.0in}
\setlength{\oddsidemargin}{0.15in}
\newcommand{\hs}{\hspace{1.0in}}
\newcommand{\hhs}{\hspace{0.5in}}
\newcommand{\hqs}{\hspace{0.25in}}
\newcommand{\hes}{\hspace{0.125in}}
\newcommand{\bt}{\begin{tabbing}}
\newcommand{\et}{\end{tabbing}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\ds}{\displaystyle}
\newcommand{\ver}{\verbatim}
\newtheorem{lemma}{Lemma}
\begin{document}
\title{SL11F -- MARCO FORTRAN Library Routine Document}
\author{}
\date{}
\maketitle
{\scriptsize NOTE: before using this routine, please read the 
appropriate implementation document to check the interpretation of {\bf bold 
italicised } terms and other implementation-dependent details. The routine 
name may be precision dependent.}
\section{ Purpose}
SL11F finds a specified eigenvalue of a differential equation  eigenvalue 
problem posed in the form of an eigen-parameter dependent Hamiltonian 
system with suitable separated boundary conditions. The method used is 
a shooting method based on matrix oscillation theory, using the algorithm
described in \cite{kn:Marletta}.

\section{ Specification}
\bt 
\hs {\tt SUBROUTINE SL11F(ELAM,EPS,A,B,K,N,TOL,USUB1,USUB2,KNOBS,} \\
\hhs\hqs\hes {\tt \&}\hs\hhs {\tt WORK,IWORK,LWK,ILWK,IFAIL) } \\
\hs {\tt INTEGER ILWK, LWK(ILWK),IWORK,IFAIL, K, KNOBS, N }\\
\hs {\tt {\bf real } ELAM, EPS, A, B, TOL, WORK(IWORK)} \\
\hs {\tt EXTERNAL USUB1,USUB2 }
\et
\section{ Description \label{sec:Descrip}}
SL11F finds a specified eigenvalue $\lambda_{k} $ of a Hamiltonian eigenvalue problem of
the form
\be
\left( \begin{array}{c} u' \\ v' \end{array} \right)  = 
\left(\begin{array}{ll} S_{1,1}(x,\lambda) & S_{1,2}(x,\lambda) \\
                        S_{1,2}^{T}(x,\lambda) & S_{2,2}(x,\lambda) \end{array} \right)
\left( \begin{array}{cc} u \\ v \end{array} \right), \;\;\; x\in(a,b) \label{eq:1} \ee
\be A_{1}^{T}u(a) + A_{2}^{T}v(a)  =  0, \label{eq:2} \ee
\be B_{1}^{T}u(b) + B_{2}^{T}v(b)  =  0, \label{eq:3} \ee
where the $S_{i,j}$ are real $n\times n$ matrix functions of $x$ and $\lambda$, with 
$S_{1,1}$ and $S_{2,2}$ symmetric, and $A_{1}$, $A_{2}$, $B_{1}$ and $B_{2}$ are 
$n\times n$ matrices satisfying the conjointness conditions
\be A_{1}^{T}A_{2}-A_{2}^{T}A_{1} = 0, \;\;\; B_{1}^{T}B_{2}-B_{2}^{T}B_{1} = 0, \ee
and the rank conditions
\be \mbox{rank}(A_{1}|A_{2}) = n, \;\;\; \mbox{rank}(B_{1}|B_{2}) = n. \label{eq:rank} \ee
Here $(A_{1}|A_{2})$ denotes the $n\times 2n$ matrix whose first $n$ columns are the
columns of $A_{1}$ and whose $n+1$st to $2n$th columns are the columns of $A_{2}$.

There is no a-priori reason for the eigenvalue problem which we have just posed to have
any solutions. For a number of well-known cases (see, e.g., Atkinson \cite{kn:Atkinson}) 
the existence of eigenvalues is established. These cases include Hamiltonian 
systems derived from $2m$-th order linear formally self-adjoint differential 
operators.

The eigenvalue indexing used by SL11F is that which ensures that for Sturm-Liouville 
problems and for problems derived from $2m$-th order linear self-adjoint operators, 
the eigenvalues will be an infinite sequence
\[ \lambda_{0} \leq \lambda_{1} \leq \lambda_{2} \leq \cdots \]
For other problems, possibly with nonlinear dependence of (\ref{eq:1}) 
on the eigenparameter $\lambda$, it cannot be guaranteed that an eigenvalue will exist 
with any particular index. It is possible that the first few eigenvalues will be missing, 
or that there will be only finitely many eigenvalues; it is even possible for the eigenvalues 
to form a sequence which extends both to $-\infty$ and to $+\infty$; such a sequence requires 
an indexing $(\lambda_{k})_{k=-\infty}^{+\infty}$.

The shooting method used by SL11F integrates a differential equation to find a function
from which the position of the eigenvalues may be deduced by a rootfinding process. 
The integration of the differential equation requires a tolerance TOL, which is supplied
by the user. The accuracy of the eigenvalue approximation obtained depends on TOL, 
although the error in the eigenvalue approximation may well exceed TOL. It is therefore
recommended that SL11F be called with two different values of TOL in order to assess the
accuracy of the eigenvalue approximations obtained.

\section{References}
See the bibliography at the end of this document.

\section{ Parameters}
\begin{description}   
\item[ELAM] -- {\bf real} scalar (input/output).
\begin{description}   
\item[On entry:] ELAM must 
specify an initial approximation to the eigenvalue $ \lambda_{k} $ which is 
sought. This need not be at all accurate since the routine will always find an 
approximation to $\lambda_{k} $ even if this is not the eigenvalue closest to 
the user's initial approximation. 
\item[On exit:] ELAM will contain the approximation to $\lambda_{k}$.
\end{description}
\item[EPS] -- {\bf real} scalar (input).
\begin{description}   
\item[On entry:] EPS must be set to a {\bf strictly positive} order-of-maginitude 
estimate of the error in the initial estimate ELAM. Over-estimates are definitely 
preferable to under-estimates.
\end{description}
\item[A] -- {\bf real} (input). 
\begin{description}   
\item[On entry:]  A must specify the left-hand endpoint $a$ of the interval $(a,b)$
 on which the problem is posed.
\end{description}
\item[B] -- {\bf real} (input). 
\begin{description}   
\item[On entry:] B must specify the right-hand endpoint $b$ of the interval $(a,b)$
 on which the problem is posed. B $>$ A is required.
\end{description}
\item[K] -- INTEGER (input). 
\begin{description}   
\item[On entry:] K must specify the index $k$ of the eigenvalue sought. To check for a 
sequence of eigenvalues decreasing to $-\infty$, call SL11F with a large negative value of
K (e.g. -100) and a large value of EPS (e.g. EPS=10); if failure occurs with IFAIL = 2, 
then such a sequence is unlikely.
\end{description}
\item[N] -- INTEGER (input).
\begin{description}   
\item[On entry:] N must specify the dimension of the matrices $S_{i,j}$ in (\ref{eq:1}).
\end{description}
\item[TOL] -- {\bf real} (input). 
\begin{description}   
\item[ On entry:] TOL should be strictly positive. SL11F will use TOL as an error 
 control parameter during the shooting process. Reducing TOL should reduce the 
 error in the eigenvalue approximation.
\end{description}
\item[USUB1] -- SUBROUTINE, supplied by the user. USUB1 must supply the 
values of the coefficient matrices $S_{i,j}$ in (\ref{eq:1}).
The specification of USUB1 is as follows. 
\bt
\hs {\tt SUBROUTINE USUB1(X,ELAM,S11,S12,S21,S22,N)} \\
\hs {\bf real} {\tt X, ELAM} \\
\hs {\tt INTEGER N} \\
\hs {\bf real} {\tt S11(N,N), S12(N,N), S21(N,N), S22(N,N) } 
\et
\begin{description}   
\item[X] -- {\bf real} (input).
\begin{description}   
\item[On entry:] X contains the value of a point in $(a,b)$. X must not be 
changed by USUB1.
\end{description}
\item[ELAM] -- {\bf real} (input).
\begin{description}   
\item[On entry:] ELAM contains the current value of $\lambda$. ELAM must 
not be changed by USUB1.
\end{description}
\item[S11] -- {\bf real} ARRAY of DIMENSION (N,N) (output).
\begin{description}   
\item[On exit:] S11 must contain the value of $S_{1,1}$ at the point X for the
current value of $\lambda$ specified by ELAM.
\end{description}
\item[S12] -- {\bf real} ARRAY of DIMENSION (N,N) (output).
\begin{description}   
\item[On exit:] S12 must contain the value of $S_{1,2}$ at the point X for the
current value of $\lambda$ specified by ELAM.
\end{description}
\item[S21] -- {\bf real} ARRAY of DIMENSION (N,N) (output).
\begin{description}   
\item[On exit:] S21 must contain the value of $S_{1,2}^{T}$ at the point X for the
current value of $\lambda$ specified by ELAM.
\end{description}
\item[S22] -- {\bf real} ARRAY of DIMENSION (N,N) (output).
\begin{description}   
\item[On exit:] S22 must contain the value of $S_{2,2}$ at the point X for the
current value of $\lambda$ specified by ELAM.
\end{description}
\item[N] -- INTEGER (input).
\begin{description}   
\item[On entry:] N specifies the value of the variable N in the user's call to
SL11F. N must not be changed by USUB1.
\end{description}
\end{description}
\item[USUB2] -- SUBROUTINE, supplied by the user. This subroutine supplies the 
boundary conditions for the problem as specified by equations 
(\ref{eq:2},\ref{eq:3}). The specification of USUB2 is as follows.
\bt 
\hs {\tt SUBROUTINE USUB2(IEND,ISING,X,U,V,NDIM,N,ELAM) } \\
\hs {\tt INTEGER IEND, NDIM, N } \\
\hs {\tt LOGICAL ISING } \\
\hs {\bf real }{\tt X, U(NDIM,N), V(NDIM,N), ELAM }
\et
\begin{description}   
\item[IEND] -- INTEGER (input). 
\begin{description}   
\item[On entry:] IEND specifies the end of the interval 
$ (a,b) $ at which the boundary conditions are to be imposed. If IEND 
= 0 then boundary conditions are required at the end $ x = a$; if IEND = 1 
then boundary conditions are required at $ x = b$. IEND must not be 
changed by the routine UUSB1.
\end{description}
\item[ISING] -- LOGICAL (output).
\begin{description}   
\item[On exit:] ISING must have the value .FALSE. This parameter is
included to allow for a future upgrade of SL11F to handle singular
problems.
\end{description}
\item[X] -- {\bf real} scalar (input).
\begin{description}   
\item[On entry:] X specifies the actual endpoint at which the 
boundary condition is required (i.e. the value A or B). Like
ISING, this parameter is included to allow for a future upgrade 
of SL11F to handle singular problems.
\end{description}
\item[U] -- {\bf real} ARRAY of DIMENSION (NDIM,N) (output).
\begin{description}   
\item[On exit:] U must specify the matrix $A_{2}$ if IEND=0
or $B_{2}$ if IEND=1.
\end{description}
\item[V] -- {\bf real} ARRAY of DIMENSION (NDIM,N) (output).
\begin{description}   
\item[On exit:] V must specify the matrix $-A_{1}$ if IEND=0
or $-B_{1}$ if IEND=1.
\end{description}
\item[NDIM] -- INTEGER (input).
\begin{description}   
\item[On entry:]  NDIM specifies the first dimensions of the 
arrays U and V. The value of NDIM must not be changed by USUB2.
\end{description}
\item[N] -- INTEGER (input).
\begin{description}   
\item[On entry:] N specifies the value of the variable $N$ in the user's call to
SL11F. N must not be changed by USUB2.
\end{description}
\item[ELAM] -- {\bf real} scalar (input).
\begin{description}   
\item[On entry:] ELAM specifies the current value of $\lambda$, allowing
for $\lambda$-dependent boundary conditions. The value of ELAM must not 
be changed by USUB2.
\end{description}
\end{description}
\item[KNOBS] -- INTEGER (input). 
\begin{description}   
\item[On entry:] KNOBS allows the user to increase the maximum number of 
miss-distance function evaluations used to find the eigenvalue. The 
maximum number of miss-distance function evaluations is max(60,KNOBS). 
In general, it is recommended that KNOBS be assigned the value 0. 
\end{description}                                         
\item[WORK] -- {\bf real} array of DIMENSION IWORK (workspace). 
\item[IWORK] -- INTEGER (input). 
\begin{description}   
\item[On entry:] IWORK must specify the dimension of WORK as declared in the
calling (sub)program; IWORK $\geq $ N*(16+73*N)+65.
\end{description}
\item[LWK] -- INTEGER array of DIMENSION ILWK (workspace). 
\item[ILWK] -- INTEGER (input). 
\begin{description}   
\item[ON entry:] ILWK must specify the dimension of LWK as declared in the 
calling (sub)program; ILWK $\geq $ N.
\end{description}
\item[IFAIL] -- INTEGER (input/output).
\begin{description}   
\item[On entry:] IFAIL must be set to -1, 0 or 1. For users not familiar with 
this parameter (described in Chapter P01 of the NAG Library documentation) the 
recommended value is 0. 
\item[On exit:] IFAIL = 0 unless the routine detects an error (see Section \ref{section:6}).
\end{description}
\end{description}

\section{ Error Indicators and Warnings \label{section:6}}
Errors detected by the routine:
\begin{description}   
\item[IFAIL = 1.] 
A parameter error on input. This error may be trigerred by all of the following:
N $<$ 2, TOL $\leq$ 0, EPS $\leq$ 0 and B $\leq$ A, IWORK or LWK less than
the minimum values specified above.  
If IFAIL was -1 or 0 on entry then an error message will give more details.

\item[IFAIL = 2.] 

Too many attempts have been made to find the required eigenvalue. Try increasing
the parameter KNOBS (to a value greater than 60) or increase TOL.

\item[IFAIL = 3.]

An error has occurred in F02BJF when converting the boundary conditions supplied
by the user into the variables used by SL11F. This may occur if the matrices
$(A_{1}|A_{2})$ or $(B_{1}|B_{2})$ are almost rank deficient, or if they have 
some extremely large entries. Check the routine USUB2 for coding errors.

\item[IFAIL = 4.]

The boundary conditions most recently supplied do not appear to satisfy the
rank condition (\ref{eq:rank}). Check the routine USUB2 for coding errors.

\item[IFAIL = 5.]

An $H$-matrix (see \cite{kn:Marletta}) has occurred whose exponential cannot 
be accurately computed. Check the routine USUB1 for coding errors and consider 
reducing the value of TOL.

\item[IFAIL = 6.]

A failure has occurred on a call to F02AJF; the eigenvalues of the central 
miss-distance unitary matrix $\Theta_{R}^{*}\Theta_{L}$ (see \cite{kn:Marletta})
cannot be computed accurately. Check the routine USUB1 for coding errors and
consider reducing the value of TOL.

\item[IFAIL = 7.]

An $H$-matrix (see \cite{kn:Marletta}) has occurred which cannot be 
accurately diagonalised. Check the routine USUB1 for coding errors and consider
reducing the value of TOL.

\item[IFAIL = 8.] 

Too many unsuccessful attempts have been made to find a suitable initial 
stepsize to integrate the shooting equations. Check the routine USUB1 for
coding errors and consider changing the value of TOL (increasing TOL is
more likely to yield success).

\item[IFAIL = 9.]

The integration of the shooting equations has ground to a halt because
a suitable stepsize cannot be found for the next step. Check the routine
USUB1 for coding errors and consider changing the value of TOL (increasing
TOL is more likely to yield success).
\end{description}

\section{Accuracy}
The error in the eigenvalue approximation obtained is found empirically to be
a modest multiple of TOL for most problems. However the user is always advised 
to run SL11F at two different tolerances to obtain a more reliable indication
of the accuracy which is being achieved.

\section{Run times}
The run times are proportional to $N^{3}$ and increase with decreasing TOL.

\section{ Keywords}
Eigenvalue Problems, Hamiltonian Equations, Sturm-Liouville Problems,
Matrix Oscillation Theory, Shooting Methods.


\section{ Example}
We consider the problem
\[ -y'' + Q(x)y = \lambda y, \;\;\; x\in (0.1,1], \;\;\; y(0.1) = y(1) = 0, \]
where $Q(x)$ is the following $4\times 4$ matrix:
\be Q_{i,j}(x) = \frac{1}{\max(i,j)}\cos x + \frac{\delta_{i,j}}{x^{i}}. \label{eq:exeq} \ee
This may be written in the form of a Hamiltonian system with $S_{1,1}=\lambda I - Q$,
$S_{1,2}=S_{2,1}=0$, and $S_{2,2}=I$, where $I$ denotes the identity matrix.
The symbol $\delta_{i,j}$ in (\ref{eq:exeq}) is the usual Kronecker $\delta$:
\[ \delta_{i,j} = \left\{ \begin{array}{ll} 0 & \mbox{for $i\neq j$} \\
                                            1 & \mbox{for $i=j$} 
                          \end{array}{ll} \right. \]
We ask SL11F to compute an approximation to $\lambda_{3}$. 

\subsection{ Program Text}
\noindent {\scriptsize WARNING: This {\bf double precision} program may 
require amendment for certain implementations. The results produced may not be 
the same. If in doubt, please seek further advice. }
\begin{verbatim}

C     SL11F EXAMPLE PROGRAM TEXT.
C     MARK n RELEASE. MARCO MARLETTA COPYRIGHT 1992.
C     .. PARAMETERS ..
      INTEGER IWORK,NMAX
      PARAMETER (NMAX=10,IWORK=NMAX*(16+73*NMAX)+65)
C     .. LOCAL SCALARS ..
      INTEGER IFAIL,K,KNOBS,N,NOUT
      DOUBLE PRECISION ELAM,EPS,A,B,TOL
C     .. LOCAL ARRAYS ..
      INTEGER LWK(NMAX)
      DOUBLE PRECISION WORK(IWORK)
C     .. EXTERNAL SUBROUTINES ..
      EXTERNAL USUB1,USUB2
C     ..
      
      DATA N /4/
      DATA NOUT /6/

      ELAM = 0.0D0
      EPS = 1.0D0
      TOL = 1.D-7
      A = 0.1D0
      B = 1.0D0
      K = 3
      IFAIL = 1

      CALL SL11F(ELAM,EPS,A,B,K,N,TOL,USUB1,USUB2,KNOBS,WORK,
     &           IWORK,LWK,NMAX,IFAIL)

      IF (IFAIL.NE.0) GOTO 100

      WRITE(NOUT,9000) K,ELAM
9000  FORMAT(' LAMBDA(',I2,')= ',G18.8)
      STOP

100   WRITE(NOUT,9010) IFAIL
9010  FORMAT(' Abnormal exit from SL11F, IFAIL = ',I2)

      WRITE(NOUT,9020) ELAM
9020  FORMAT(' ELAM = ',G18.8)

      STOP

      END

      SUBROUTINE USUB1(X,ELAM,S11,S12,S21,S22,N)
      INTEGER I,J,N
      DOUBLE PRECISION X,ELAM,S11(N,N),S12(N,N),S21(N,N),S22(N,N)

      DO 10 J=1,N
         DO 20 I=1,N
            S12(I,J) = 0.0D0
            S21(I,J) = 0.0D0
            S22(I,J) = 0.0D0
            S11(I,J) = -COS(X)/DBLE(MAX(I,J))       
20       CONTINUE
         S22(J,J) = 1.0D0
         S11(J,J) = S11(J,J) + ELAM - X**(-J)
10    CONTINUE

      RETURN
      END

      SUBROUTINE USUB2(IEND,ISING,X,U,V,NDIM,N,ELAM)
      INTEGER I,J,IEND,NDIM,N
      LOGICAL ISING
      DOUBLE PRECISION ELAM,X,U(NDIM,N),V(NDIM,N)

      ISING = .FALSE.
      DO 10 J=1,N
         DO 20 I=1,N
            U(I,J) = 0.0D0
            V(I,J) = 0.0D0
20       CONTINUE    
         V(J,J) = 1.0D0
10    CONTINUE

C REMARK: WE COULD ACTUALLY ASSIGN V TO BE ANY NON-SINGULAR MATRIX.

      RETURN
      END

\end{verbatim}

\subsection{Results}
\begin{verbatim}

      LAMBDA( 3) = 26.920694

\end{verbatim}

\begin{thebibliography}{99}
\bibitem{kn:Atkinson} F.V. Atkinson, {\em Discrete and Continuous Boundary
Value Problems}, Academic Press, New York (1964).
\bibitem{kn:Marletta} M. Marletta, {\em Numerical Solution of Eigenvalue
Problems for Hamiltonian Systems}, University of Leicester Mathematics
Department Technical Report 1992/17.
\end{thebibliography}
\end{document}

.