\section{The Strongly Implicit (SIP) algorithm\label{sip}}
Given a rectangular grid with natural ordering of points $(i,j)$ as shown below \cite{8},\\
\diagram{stencil.eps}{width=8cm}{The finite difference stencil}{stencil}
a finite difference approximation of sides $I+2$ and $J+2$ (where $i =  2, \dots . I+1$ and $j = 2,\dots .J+1$) generates one linear equation for each grid point $(i,j)$, giving $K = (I+1)(J+1)$ linear equations in total. 
The basic SIP algorithm assumes a five-point difference approximation over a rectangular mesh for the partial differential equations. This results in an $A$ matrix of pentadiagonal structure (tridiagonal matrix with two super-diagonals at a set distance $J+1$ from the main diagonal), with a bandwidth of $2(J+1)$ \cite{14, 1, 13}. For large $K$ (total mesh size) the standard $LU$ decomposition fills in a very large number of the zeros between the super and tridiagonals. From the matrix form  $A u  =  q$
the general equation at grid point $(i,j)$ is: 
\begin{equation}
 B_{i,j}u_{i-1,j} + D_{i,j}u_{i,j-1} + E_{i,j}u_{i,j} + F_{i,j}u_{i,j+1}+H_{i,j} u_{i+ 1,j}=q_{i,j} 
\end{equation}

The crux of the method is to modify the $A$ matrix with a matrix $N$ such that in factorising the $(A + N)$ matrix, the $L$ and $U$ resultant matrices each have only three non-zero elements in each row, thus avoiding matrix fill and achieving a limited computational build-up that approaches $\Theta(n )$. 


The forms of the  $L$ and $U$ matrices are shown below:-\\ 
\diagram{lu.eps}{width=12cm}{The form of the $L$ and $U$ matrices.}{lu}

 
The nonzero elements of lower triangular matrix L are in diagonals $d_{i,j}$, $c_{i,j}$ and $b_{i,j}$ corresponding to the $B_{i,j}$, $D_{i,j}$  and $E_{i,j}$  diagonals of $A$. Similarly, the upper triangular matrix $U$ has nonzero elements in diagonals corresponding to $E_{i,j}$, $F_{i,j}$  and $H_{i,j}$  diagonals of $A$ , where the main diagonal corresponding to $E_{i,j}$  is unity. The product $LU$ (the form of which is shown in figure~\ref{lutotal}) gives a matrix with seven nonzero diagonals, five corresponding to the diagonals $(B, D, E, F, H)$ of $A$  and two ($C$ and $G$ ) which fall just inside the outer diagonals $B$ and $H$, and which make up $N$. The algorithm is designed so that the terms $Nu$ from the augmented matrix are small.\\ 
\diagram{lutotal.eps}{width=8cm}{The form of the $LU$ product matrix.}{lutotal}

 The complete set of equations relating $LU$ to $(A + N)$ is:-
\begin{eqnarray}
b_{i,j}&=&B_{i,j} \label{B}\\
c_{i,j}&=&D_{i,j}\label{D} \\
d_{i,j}+b_{i,j}f_{i,j-1}+c_{i,j}e_{i-1,j}&=&E_{i,j} \label{E}\\
d_{i,j}e_{i-1,j}&=&F_{i,j}\label{F}\\
d_{i,j}f_{i,j-1}&=&H_{i,j}\label{H}\\
b_{i,j}e_{i,j-1}&=&C_{i,j}\label{C}\\
c_{i,j}f_{i-1,j}&=&G_{i,j}\label{G}
\end{eqnarray} 
For grid point $(i,j)$, there are five unknowns $(b_{i,j},c_{i,j},d_{i,j},e_{i,j},f_{i,j})$ but seven equations and therefore the system of equations is over-determined with a family of possible solutions. The simplest definition of  $(A + N)$ which can be factored into $LU$ comes from ignoring equation\ref{C} and equation\ref{G} and solving the five variables $b_{i,j}$   to $f_{i,j}$ using equations \ref{B}$\dots$\ref{H} by elimination, starting with $b_{1,1}$ (all zero subscripted variables are zero). 

Stone\cite{14}, found that the resultant $(A + N)$ matrix did not give rapid convergence on his heat conduction test problems and re-defined the $N$ matrix so that the $C_{i,j}$ and $G_{i,j}$ components are made smaller by subtracting a closely equivalent expression from them, using simple Taylor series expansions and an acceleration parameter $\alpha \ni 0<\alpha <1$ (which may be a function of the chosen grid). However, the sequential version of the SIP algorithm in the semi-conductor device simulator uses the simpler equations \ref{B}$\dots$\ref{H}   to calculate the $LU$ coefficients, with apparently satisfactory convergence properties for this particular problem. 

A SIP iteration can be defined from $Au = q$ by adding Nu to both sides of the equation and $(Au-Au)$ to the rhs giving:
\begin{equation}
 (A + N)u^{(n + 1)} = (A + N)u^{(n)}  +  (q - A u^{(n)}) \label{AN}
\end{equation}
 
Substitute $LU$ for $(A + N)$ in equation~\ref{AN} and let $d^{(n)} = u^{(n+1)}  - u^{(n)}$  and  $R^{(n)}  =  (q - A u^{(n)})$, then a complete iteration cycle consists of solving:  $LUd^{(n)}  =  R^{(n)}$ followed by $u^{(n+1)}   =   d^{(n)}+u^{(n)}$. 
 
The SIP algorithm consists of four steps:-
\begin{enumerate}
\item An initial factorisation of the A-matrix into lower and upper triangular matrices $L$ and $U$
\item  A Forward substitution step  $v^{(n)}  =  L^{-1}R^{(n)} $
\item  A Backward substitution step      $d^{(n)}  =  U^{-1}v^{(n)}$ 
\item  A maximum relative error check and new solution value update $u^{(n+1)}  =  d^{(n)} + u^{(n)}$, plus update  of residual $R^{(n)}$ from $(q - A u^{(n)})$ 
\end{enumerate} 
Steps 2 , 3 and 4 are repeated until the maximum relative difference between the old and new solution values $d^{(n)}$ are within the defined tolerance or the maximum number of iterations has been reached. \par
 ~\par
