%
% part of numlib docs. In time this won't be a standalone pdf anymore, but part of a larger part.
% for now I keep it directly compliable. Uses fpc.sty from fpc-docs pkg.
%

\documentclass{report}

\usepackage{fpc}
\lstset{%
  basicstyle=\small,
  language=delphi,
  commentstyle=\itshape,
  keywordstyle=\bfseries,
  showstringspaces=false,
  frame=
}

\makeindex

\newcommand{\FunctionDescription}{\item[Description]\rmfamily}
\newcommand{\Dataorganisation}{\item[Data Struct]\rmfamily}
\newcommand{\DeclarationandParams}{\item[Declaration]\rmfamily}
\newcommand{\References}{\item[References]\rmfamily}
\newcommand{\Method}{\item[Method]\rmfamily}
\newcommand{\Precision}{\item[Precision]\rmfamily}
\newcommand{\Remarks}{\item[Remarks]\rmfamily}
\newcommand{\Example}{\item[Example]\rmfamily}
\newcommand{\ProgramData}{\item[Example Data]\rmfamily}
\newcommand{\ProgramResults}{\item[Example Result]\rmfamily}

\makeatletter
\@ifpackageloaded{tex4ht}{%
\newcommand{\NUMLIBexample}[1]{
\par \file{\textbf{Listing:} \exampledir/#1.pas}%
\HCode{<HR/>}%
\listinginput[9999]{5000}{\exampledir/#1.pas}%
\HCode{<HR/>}%
}%
}{% else ifpackageloaded
\newcommand{\NUMLIBexample}[1]{%
\par \file{\textbf{Listing:} \exampledir/#1.pas}%
\lstinputlisting{\exampledir/#1.pas}%
}% End of FPCExample
}% End of ifpackageloaded.
\makeatother
%

\begin{document}
\FPCexampledir{../examples}
\section{general comments}

\textbf{Original comments:} \\ 
The used floating point type \textbf{real} depends on the used version,
see the general introduction for more information. You'll need to USE
units typ an inv to use these routines. 

\textbf{MvdV notes:} \\
Integers used for parameters are of type "ArbInt" to avoid problems with
systems that define integer differently depending on mode. 

Floating point values are of type "ArbFloat" to allow writing code
that is independent of the exact real type. (Contrary to directly
specifying single, real, double or extended in library and
examples). Typ.pas and the central includefile have some conditional
code to switch between floating point types.

These changes were already prepared somewhat when I got the lib, but weren't
consequently applied. I did that while porting to FPC.

\section{Unit inv}


\begin{procedure}{invgen}

\FunctionDescription

Procedure to calculate the inverse of a matrix.

\Dataorganisation

The procedure assumes that the calling program has declared a two
dimensional matrix containing the matrix elements in a square
partial block.

\DeclarationandParams

\lstinline|procedure invgen(n, rwidth: ArbInt; var ai: ArbFloat;var term: ArbInt);|

\begin{description}
 \item[n: integer] \mbox{ } \\
    The parameter {\bf n} describes the size of the matrix
 \item[rwidth: integer] \mbox{} \\
    The parameter {\bf rwidth} describes the declared rowlength of the two dimensional
    matrix.
 \item[var ai: real] \mbox{} \\
    The parameter {\bf ai} must contain to the two dimensional array element
    that is top-left in the matrix.
    After the function has ended successfully (\textbf{term=1}) then
    the input matrix has been changed into its inverse, otherwise the contents 
    of the input matrix are undefined.
 \item[var term: integer]  \mbox{} \\
    After the procedure has run, this variable contains the reason for 
    the termination of the procedure:\\
      {\bf term}=1, successful termination, the inverse has been calculated
      {\bf term}=2, inverse matrix could not be calculated because the matrix
		    is (very close to) being singular.
      {\bf term}=3, wrong input n$<$1
\end{description}
\Remarks
  This procedure does not do array range checks. When called with invalid
  parameters, invalid/nonsense responses or even crashes may be the result.

\Example

Calculate the inverse of 
\[
 A=
 \left(
 \begin{array}{rrrr}
   4 & 2 & 4 & 1  \\
  30 & 20 & 45 & 12 \\
  20 & 15 & 36 & 10 \\
  35 & 28 & 70 & 20
 \end{array}
 \right)
 .
\]

Below is the source of the invgenex demo that demonstrates invgenex,
some routines of iom were used to construct matrices.

\NUMLIBexample{invgenex}

\ProgramData

The input datafile looks like:
\begin{verbatim}
 4  2  4  1
30 20 45 12
20 15 36 10
35 28 70 20
\end{verbatim}

\ProgramResults

\begin{verbatim}
program results invgenex

A =
 4.0000000000000000E+0000   2.0000000000000000E+0000
 3.0000000000000000E+0001   2.0000000000000000E+0001
 2.0000000000000000E+0001   1.5000000000000000E+0001
 3.5000000000000000E+0001   2.8000000000000000E+0001

 4.0000000000000000E+0000   1.0000000000000000E+0000
 4.5000000000000000E+0001   1.2000000000000000E+0001
 3.6000000000000000E+0001   1.0000000000000000E+0001
 7.0000000000000000E+0001   2.0000000000000000E+0001

term= 1

inverse of A =
 4.0000000000000000E+0000  -2.0000000000000000E+0000
-3.0000000000000000E+0001   2.0000000000000000E+0001
 2.0000000000000000E+0001  -1.5000000000000000E+0001
-3.4999999999999999E+0001   2.7999999999999999E+0001

 3.9999999999999999E+0000  -1.0000000000000000E+0000
-4.4999999999999999E+0001   1.2000000000000000E+0001
 3.5999999999999999E+0001  -1.0000000000000000E+0001
-6.9999999999999999E+0001   2.0000000000000000E+0001
\end{verbatim}

\Precision

The procedure doesn't supply information about the precision of the
operation after termination.

\Method

The calculation of the inverse is based on LU decomposition with partial
pivoting.

\References

 Wilkinson, J.H. and Reinsch, C.\\
 Handbook for Automatic Computation.\\
 Volume II, Linear Algebra.\\
 Springer Verlag, Contribution I/7, 1971.

\end{procedure}

\begin{procedure}{invgpd}

\FunctionDescription

Procedure to calculate the inverse of a symmetrical, positive
definitive matrix

\Dataorganisation The procedure assumes that the calling program has
declared a two dimensional matrix containing the matrix elements in
a square partial block.

\DeclarationandParams

\lstinline|procedure invgpd(n, rwidth: ArbInt; var ai: ArbFloat; var term: ArbInt);|

\begin{description}
 \item[n: integer] \mbox{ } \\
    The parameter {\bf n} describes the size of the matrix
 \item[rwidth: integer] \mbox{} \\
    The parameter {\bf rwidth} describes the declared row length of the two dimensional
    matrix.
 \item[var ai: real] \mbox{} \\
    The parameter {\bf ai} must contain to the two dimensional array element
    that is top-left in the matrix.
    After the function has ended successfully (\textbf{term=1}) then
    the input matrix has been changed into its inverse, otherwise the contents 
    of the input matrix are undefined.
 \item[var term: integer]  \mbox{} \\
    After the procedure has run, this variable contains the reason for 
    the termination of the procedure:\\
      {\bf term}=1, successful termination, the inverse has been calculated
      {\bf term}=2, inverse matrix could not be calculated because the matrix
		    is (very close to) being singular.
      {\bf term}=3, wrong input n$<$1
\end{description}
\Remarks

\begin{itemize}
\item Only the bottom left part of the matrix $A$ needs to be passed.
\item \textbf{Warning} This procedure does not do array range checks. When called with invalid
parameters, invalid/nonsense responses or even crashes may be the result.
\end{itemize}

\Example

Calculate the inverse of 
\[
 A=
 \left(
 \begin{array}{rrrr}
   5 & 7 & 6 & 5  \\
   7 & 10 & 8 & 7 \\
   6 & 8 & 10 & 9 \\
   5 & 7 & 9 & 10
 \end{array}
 \right)
 .
\]

\NUMLIBexample{invgpdex}

\ProgramData

\begin{verbatim}
5
7 10
6 8 10
5 7 9 10
\end{verbatim}

\ProgramResults
\begin{verbatim}

program results invgpdex

A =
 5.0000000000000000E+0000   7.0000000000000000E+0000
 7.0000000000000000E+0000   1.0000000000000000E+0001
 6.0000000000000000E+0000   8.0000000000000000E+0000
 5.0000000000000000E+0000   7.0000000000000000E+0000

 6.0000000000000000E+0000   5.0000000000000000E+0000
 8.0000000000000000E+0000   7.0000000000000000E+0000
 1.0000000000000000E+0001   9.0000000000000000E+0000
 9.0000000000000000E+0000   1.0000000000000000E+0001

term= 1

inverse of A =
 6.8000000000000000E+0001  -4.1000000000000000E+0001
-4.1000000000000000E+0001   2.5000000000000000E+0001
-1.7000000000000000E+0001   1.0000000000000000E+0001
 1.0000000000000000E+0001  -6.0000000000000000E+0000

-1.7000000000000000E+0001   1.0000000000000000E+0001
 1.0000000000000000E+0001  -6.0000000000000000E+0000
 5.0000000000000000E+0000  -3.0000000000000000E+0000
-3.0000000000000000E+0000   2.0000000000000000E+0000
\end{verbatim}

\Precision

The procedure doesn't supply information about the precision of the
operation after termination.

\Method

The calculation of the inverse is based on Cholesky decomposition
for a symmetrical positive definitive matrix.

\References

 Wilkinson, J.H. and Reinsch, C.\\ Handbook for Automatic Computation.\\
 Volume II, Linear Algebra.\\ Springer Verlag, Contribution I/7, 1971.

\end{procedure}

\begin{procedure}{invgsy}

\FunctionDescription

Procedure to calculate the inverse of a symmetrical matrix.

\Dataorganisation

The procedure assumes that the calling program has declared a two
dimensional matrix containing the matrix elements in (the bottom
left part of) a square partial block.

\DeclarationandParams

\lstinline|procedure invgsy(n, rwidth: ArbInt; var ai: ArbFloat;var term:ArbInt);|

\begin{description}
 \item[n: integer] \mbox{ } \\
    The parameter {\bf n} describes the size of the matrix
 \item[rwidth: integer] \mbox{} \\
    The parameter {\bf rwidth} describes the declared row length of the two dimensional
    matrix.
 \item[var ai: real] \mbox{} \\
    The parameter {\bf ai} must contain to the two dimensional array element
    that is top-left in the matrix.
    After the function has ended successfully (\textbf{term=1}) then
    the input matrix has been changed into its inverse, otherwise the contents
    of the input matrix are undefined.
 \item[var term: integer]  \mbox{} \\
    After the procedure has run, this variable contains the reason for
    the termination of the procedure:\\
      {\bf term}=1, successful termination, the inverse has been calculated
      {\bf term}=2, inverse matrix could not be calculated because the matrix
                    is (very close to) being singular.
      {\bf term}=3, wrong input n$<$1

\end{description}

\Remarks
\begin{itemize}
\item Only the bottom left part of the matrix $A$ needs to be passed.
\item \textbf{Warning} This procedure does not do array range checks. When called with invalid
parameters, invalid/nonsense responses or even crashes may be the result.
\end{itemize}

\Example

  Calculating the inverse of
\[
 A=
 \left(
 \begin{array}{rrrr}
   5 & 7 & 6 & 5  \\
   7 & 10 & 8 & 7 \\
   6 & 8 & 10 & 9 \\
   5 & 7 & 9 & 10
 \end{array}
 \right)
 .
\]

\NUMLIBexample{invgsyex}


\ProgramData

\begin{verbatim}
5
7 10
6 8 10
5 7 9 10
\end{verbatim}

\ProgramResults

\begin{verbatim}
program results invgsyex

A =
 5.0000000000000000E+0000   7.0000000000000000E+0000
 7.0000000000000000E+0000   1.0000000000000000E+0001
 6.0000000000000000E+0000   8.0000000000000000E+0000
 5.0000000000000000E+0000   7.0000000000000000E+0000

 6.0000000000000000E+0000   5.0000000000000000E+0000
 8.0000000000000000E+0000   7.0000000000000000E+0000
 1.0000000000000000E+0001   9.0000000000000000E+0000
 9.0000000000000000E+0000   1.0000000000000000E+0001

term= 1

inverse of A =
 6.8000000000000001E+0001  -4.1000000000000001E+0001
-4.1000000000000001E+0001   2.5000000000000000E+0001
-1.7000000000000000E+0001   1.0000000000000000E+0001
 1.0000000000000000E+0001  -6.0000000000000001E+0000

-1.7000000000000000E+0001   1.0000000000000000E+0001
 1.0000000000000000E+0001  -6.0000000000000001E+0000
 5.0000000000000001E+0000  -3.0000000000000000E+0000
-3.0000000000000000E+0000   2.0000000000000000E+0000
\end{verbatim}

\Precision

The procedure doesn't supply information about the precision of the
operation after termination.

\Method

The calculation of the inverse is based on reduction of a symmetrical
matrix to a tridiagonal form.

\References

Aasen, J. O. \\
On the reduction of a symmetric matrix to tridiagonal form. \\
BIT, 11, (1971), pag. 223-242.

\end{procedure}

\end{document}