\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{SL12F -- 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}
SL12F 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}, combined with coefficient approximation.

\section{ Specification}
\bt 
\hs {\tt SUBROUTINE SL12F(ELAM,EPS,A,B,K,N,USUB1,USUB2,TOL,XMESH,MESH,} \\
\hhs\hqs\hes {\tt \&}\hs\hhs {\tt WORK,IWORK,LWK,ILWK,NXTRAP,IFAIL) } \\
\hs {\tt INTEGER ILWK, LWK(ILWK),IWORK,IFAIL, K, MESH, N, NXTRAP }\\
\hs {\tt {\bf real } ELAM, EPS, A, B, TOL, XMESH(0:MESH), WORK(IWORK)} \\
\hs {\tt EXTERNAL USUB1,USUB2 }
\et
\section{ Description \label{sec:Descrip}}
SL12F 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}$ symmetric and $S_{2,2}$ symmetric positive definite, 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. 

The eigenvalue indexing used by SL12F is that which ensures that for Sturm-Liouville 
problems, the eigenvalues will be an infinite sequence
\[ \lambda_{0} \leq \lambda_{1} \leq \lambda_{2} \leq \cdots \]
With this indexing and the 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 SL12F 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 is performed by first approximating the
coefficients $S_{i,j}$ in (\ref{eq:1}) by piecewise constants over a suitable mesh. 
This enables an eigenvalue approximation to be found. The mesh is then subdivided and 
Wynn extrapolation applied until the error in the eigenvalue approximation appears to
be less than the user's tolerance TOL.

\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 SL12F 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 = 3, 
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[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. $S_{2,2}$ must be positive definite;
if it is not, then SL12F will fail.
\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
SL12F. 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 SL12F 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 SL12F 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
SL12F. 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[TOL] -- {\bf real} (input). 
\begin{description}   
\item[ On entry:] TOL should be strictly positive. SL12F will perform at most 10 Wynn
extrapolations in an attempt to ensure that the error in the eigenvalue approximation
returned is less than TOL.
\end{description}
\item[XMESH] -- {\bf real} array of DIMENSION (0:MESH) (workspace). 
\item[MESH] -- INTEGER (input/output)
\begin{description}
\item[On entry:] MESH must specify the first dimension of XMESH as declared in the
calling (sub)program. Since XMESH is used to store the mesh used by SL12F, and
since the number of mesh intervals required depends on the user's accuracy requirements,
it is impossible to specify a value of MESH which will always be sufficient. However 
for regular problems and modest tolerances (no smaller than $10^{-6}$), MESH=2000
will usually be more than adequate.
\item[On exit] MESH specifies the number of mesh intervals actually used.
\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*(8+60*N)
\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[NXTRAP] -- INTEGER (output).
\begin{description}
\item[On exit:] NXTRAP specifies the number of extrapolations which were required.
This gives some measure of the inherent difficulty of the problem concerned.
\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.] 

An error has occurred while converting the boundary conditions from the user's
variables to SL12F's variables. Check the routine USUB2, making sure that the
boundary conditions satisfy the rank condition and the conjointness condition.

\item[IFAIL = 3.]

Too many evaluations of the miss-distanve function have been made without an
eigenvalue being located. Try increasing EPS. It is also possible that no 
eigenvalue exists with the index specified.

\item[IFAIL = 4.]

The meshing routine is making little progress as the stepsize has become too
small. Try increasing the tolerance TOL.

\item[IFAIL = 5.]

The array XMESH is not large enough. Increase MESH or increase TOL.

\item[IFAIL = 6.]

The matrix $V+iU$ (see \cite{kn:Marletta}) has become rank deficient. Check 
the routines USUB1 and USUB2 for coding errors and consider reducing TOL.

\item[IFAIL = 7.]

The matrix $V+iU$ (see \cite{kn:Marletta}) has become ill-conditioned. Check 
the routines USUB1 and USUB2 for coding errors and consider reducing TOL.

\item[IFAIL = 8.] 

The eigenvalues of a $\Theta$-matrix (see \cite{kn:Marletta}) cannot be computed
accurately, possibly due to ill-conditioning of $V+iU$. Proceed as for IFAIL=7.

\item[IFAIL = 9.]

The eigenvalues of an H-matrix (see \cite{kn:Marletta}) cannot be computed
accurately, possibly due to the presence of large elements. Proceed as for IFAIL=7.

\item[IFAIL=10]

An H-matrix has become non-Hermitian, possibly due to an accumulation of 
round-off error. Proceed as for IFAIL=7.

\item[IFAIL=11] 

The coefficient $S_{2,2}$ cannot be diagonalised with sufficient accuracy. Check 
the routine USUB1 and make sure that $S_{2,2}$ is symmetric.

\item[IFAIL=12]

The coefficient $S_{2,2}$ is not positive definite. Check the routine USUB1.

\end{description}

\section{Accuracy}
For most problems the eigenvalue error will be less than the user's tolerance
TOL. Nevertheless for difficult problems it may be advisable to run SL12F
at two different tolerances to obtain an indication of the reliability of
the error estimates.

\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, Coefficient Approximation.

\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 SL12F 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     SL12F EXAMPLE PROGRAM TEXT.
C     MARK n RELEASE. MARCO MARLETTA COPYRIGHT 1992.
C     .. PARAMETERS ..
      INTEGER IWORK,NMAX,MAXMSH
      PARAMETER (NMAX=10,IWORK=NMAX*(8+60*NMAX),MAXMSH=5000)
C     .. LOCAL SCALARS ..
      INTEGER IFAIL,K,MESH,N,NOUT,NXTRAP
      DOUBLE PRECISION ELAM,EPS,A,B,TOL
C     .. LOCAL ARRAYS ..
      INTEGER LWK(NMAX)
      DOUBLE PRECISION WORK(IWORK),XMESH(0:MAXMSH)
C     .. EXTERNAL SUBROUTINES ..
      EXTERNAL USUB1,USUB2
C     ..
      
      DATA N /4/
      DATA NOUT /6/

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

      CALL SL12F(ELAM,EPS,A,B,K,N,USUB1,USUB2,TOL,XMESH,MESH,WORK,
     &           IWORK,LWK,NMAX,NXTRAP,IFAIL)

      IF (IFAIL.NE.0) GOTO 100

      WRITE(NOUT,9000) K,ELAM
9000  FORMAT(' LAMBDA(',I2,')= ',G18.8)
      WRITE (NOUT,9030) MESH,NXTRAP
9030  FORMAT(' MESHPOINTS: ',I5,' EXTRAPOLATIONS: ',I3)
      STOP

100   WRITE(NOUT,9010) IFAIL
9010  FORMAT(' Abnormal exit from SL12F, 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.920732
      MESHPOINTS:   171 EXTRAPOLATIONS: 2

\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}

.