%% LATEX TWICE TO GET THE REFERENCES RIGHT
%%\documentstyle[12pt]{article}
\documentstyle[12pt,epsf,psfig]{article}

\newcommand{\postscript}[2]
 {\setlength{\epsfxsize}{#2\hsize}
  \centerline{\epsfbox{#1}}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newtheorem{theorem}{Theorem}
\newtheorem{itlemma}{Lemma}[section] %number by section (set in \em by default)
%%%same numbering counter as itlemmas:
\newtheorem{itproposition}[itlemma]{Proposition}
\newtheorem{itcorollary}[itlemma]{Corollary}
\newtheorem{itremark}[itlemma]{Remark}
\newtheorem{itnote}[itlemma]{Note}
\newtheorem{itobservation}[itlemma]{Observation}
\newtheorem{itdefinition}[itlemma]{Definition}
\newtheorem{itexercise}[itlemma]{Exercise}
\newtheorem{itexample}[itlemma]{Example}

\newenvironment{lemma}{\begin{itlemma}\rm}{\end{itlemma}} %no-italics
\newenvironment{remark}{\begin{itremark}\rm}{\end{itremark}} %no-italics
\newenvironment{observation}{\begin{itobservation}\rm}{\end{itobservation}} %no-italics
\newenvironment{corollary}{\begin{itcorollary}\rm}{\end{itcorollary}}
\newenvironment{proposition}{\begin{itproposition}\rm}{\end{itproposition}}
\newenvironment{definition}{\begin{itdefinition}\rm}{\end{itdefinition}}
\newenvironment{exercise}{\begin{itexercise}\rm}{\end{itexercise}}
\newenvironment{example}{\begin{itexample}\rm}{\end{itexample}}

\tolerance=1600 %allow some tolerance in extending after line.
\parskip=3pt plus 1pt minus 1pt
\overfullrule=0pt %no dark lines if overfull

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This are for choices.  They take 4 or 6 args
%% (carriage returns allowed in between}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\twoif}[4]{
\left\{ \begin{array}{ll}#1&#2\\#3&#4\end{array}\right.
}
\newcommand{\threeif}[6]{
\left\{ \begin{array}{ll}#1&#2\\#3&#4\\#5&#6\end{array}\right.
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\partmlab#1#2#3#4#5#6#7#8{\bordermatrix{&#1&#2\cr#3&#5&#6\cr
					\noalign{\medskip}#4&#7&#8\cr}}
%% For 2by2 partitioned matrices; arguments are: two labels of top, side,
%% 4 entries (row by row).  If entries are themselves matrix or pmatrix,
%% better for readability of source to define first: \def\matrix1{...} etc.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\bigframebox}[1]{\framebox{\parbox{\smallerwidth}{#1}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Warning: use \noindent before!
% AND IN TEXT, after begin doc, do:
% \newlength{\smallerwidth}  %for bigboxes
% \setlength{\smallerwidth}{\textwidth}
% \addtolength{\smallerwidth}{-8pt}        %MUST BE DONE AFTER begin doc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def \beq {\begin{eqnarray}}
\def \eeq {\end{eqnarray}}
\def \beqn {\begin{eqnarray*}}
\def \eeqn {\end{eqnarray*}}

\newcommand{\rf}[1]{[\ref{#1}]}            %bibliographical references

\newcommand{\text}[1]{\hbox{\rm \ #1\ \/}}
\newcommand{\be}[1]{\begin{equation}\label{#1}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bl}[1]{\begin{lemma}\label{#1}}
\newcommand{\ble}[1]{\begin{lemmaex}\label{#1}}
\newcommand{\br}[1]{\begin{remark}\label{#1}}
\newcommand{\bn}[1]{\begin{note}\label{#1}}
\newcommand{\bo}[1]{\begin{observation}\label{#1}}
\newcommand{\bt}[1]{\begin{theorem}\label{#1}}
\newcommand{\bd}[1]{\begin{definition}\label{#1}}
\newcommand{\bp}[1]{\begin{proposition}\label{#1}}
\newcommand{\bc}[1]{\begin{corollary}\label{#1}}
\newcommand{\bfact}[1]{\begin{fact}\label{#1}}
\newcommand{\ber}[1]{\begin{exercise}\label{#1}}
\newcommand{\bex}[1]{\begin{example}\label{#1}}
\newcommand{\bem}[1]{\begin{example}\label{#1}}    %Yes, 2 different ones...
\newcommand{\ec}{\mybox\end{corollary}}
\newcommand{\efact}{\mybox\end{fact}}
\newcommand{\eer}{\mybox\end{exercise}}
\newcommand{\eex}{\mybox\end{example}}
\newcommand{\eem}{\mybox\end{example}}
\newcommand{\el}{\mybox\end{lemma}}
\newcommand{\ele}{\mybox\end{lemmaex}}
\newcommand{\er}{\mybox\end{remark}}
\newcommand{\en}{\mybox\end{note}}
\newcommand{\eo}{\mybox\end{observation}}
\newcommand{\et}{\qed\end{theorem}}
\newcommand{\ed}{\mybox\end{definition}}
\newcommand{\ep}{\mybox\end{proposition}}
\newcommand{\epr}{\end{proof}}
\newcommand{\bpr}{\begin{proof}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  These are just for statements, no "qed" box is put in.
%%  (Also for when I have to finish something after a displayed equation
%%        or itemized.)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\ecs}{\end{corollary}}
\newcommand{\eers}{\end{exercise}}
\newcommand{\eexs}{\end{example}}
\newcommand{\eems}{\end{example}}
\newcommand{\els}{\end{lemma}}
\newcommand{\eles}{\end{lemmaex}}
\newcommand{\ers}{\end{remark}}
\newcommand{\ens}{\end{note}}
\newcommand{\eos}{\end{observation}}
\newcommand{\ets}{\end{theorem}}
\newcommand{\eds}{\end{definition}}
\newcommand{\eps}{\end{proposition}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\qed}{\hfill \halmos} %put \qed at right margin
\newcommand{\mybox}{\hfill $\Box$} %put \qed at right margin (white square)

%%%\newcommand{\pmat}[2]{\left( \begin{array}{#1} #2 \end{array} \right) }
%%%use as e.g. \pmat{cc}{ a&b\\ a &c }  Note 'cc' for 2 Centered cols
%%%can also use \pmatrix{ a&b\cr a &c \cr} for plainTeX standard construct
%%%\setcounter{page}{27} should set the CURRENT page to 27

\def \R {\bb R} % reals
\def \C {\bb C} % complexes
\def \Z {\bb Z} % integers
\def \X {\bb X} % for this paper only!
\def \inter {{\rm int \,}}
\def \rank {{\rm rank \,}}
%%%\def \inf {{\rm inf \,}}  %better in Latex -- subscripts OK!
%%%\def \min {{\rm min \,}}
\def \im {{\rm im \,}}
\def \ker {{\rm ker \,}}
\def \norm #1{\Vert #1 \Vert}
\def \ip #1\by #2{\langle #1 , #2 \rangle}
\def \Sone {{\bf S}^1} % circle
\def \wigglearrow {\sim \mkern -6mu >}

\def \today {\ifcase \month \or January\or February\or  
  March\or April\or May\or June\or July\or August\or 
  September\or October\or November\or December\fi\   
  \number \day, \number \year}                       
\newcommand{\comment}[1]{}
\newcommand{\halmos}{\rule{1ex}{1.4ex}}
%%%alternative to \qed (see below)
\newenvironment{proof}{\noindent {\em Proof}.\ }{\hspace*{\fill}$\halmos$\medskip}

%%%\def\qed{\kern 1em \vrule height .5em width .5em depth 0em}
%%%(I now use new qed: in right margin)

\def\vbar{\mathchoice{\vrule height6.3ptdepth-.5ptwidth.8pt\kern-.8pt}
   {\vrule height6.3ptdepth-.5ptwidth.8pt\kern-.8pt}
   {\vrule height4.1ptdepth-.35ptwidth.6pt\kern-.6pt}
   {\vrule height3.1ptdepth-.25ptwidth.5pt\kern-.5pt}}
\def\fudge{\mathchoice{}{}{\mkern.5mu}{\mkern.8mu}}
\def\bbc#1#2{{\rm \mkern#2mu\vbar\mkern-#2mu#1}}
\def\bbb#1{{\rm I\mkern-3.5mu #1}}
\def\bba#1#2{{\rm #1\mkern-#2mu\fudge #1}}
\def\bb#1{{\count4=`#1 \advance\count4by-64 \ifcase\count4\or\bba A{11.5}\or
   \bbb B\or\bbc C{5}\or\bbb D\or\bbb E\or\bbb F \or\bbc G{5}\or\bbb H\or
   \bbb I\or\bbc J{3}\or\bbb K\or\bbb L \or\bbb M\or\bbb N\or\bbc O{5} \or
   \bbb P\or\bbc Q{5}\or\bbb R\or\bbc S{4.2}\or\bba T{10.5}\or\bbc U{5}\or
%%%   \bbb P\or\bbc Q{5}\or\bbb R\or\bba S{8}\or\bba T{10.5}\or\bbc U{5}\or
   \bba V{12}\or\bba W{16.5}\or\bba X{11}\or\bba Y{11.7}\or\bba Z{7.5}\fi}}
%
%
\renewcommand{\thefootnote}{\fnsymbol{footnote}}
%%% read as: theFOOTNOTE (i.e. the command  executed when doing a footnote) is
%%% now to be set as a FOOTNOTESYMBOL (type of number) value of counter
%%% FOOTNOTE (Instead of Roman).  See LaTeX manual pp.91-92 as well as 175.
\def\bi{\begin{itemize}}
\def\ei{\end{itemize}}
\def\ben{\begin{enumerate}}
\def\een{\end{enumerate}}
\def\st{\, | \,}

\newcommand{\myboxineq}{\eqno\Box} %put \qed at margin (white square) in eqn
\newcommand{\qedineq}{\eqno\halmos} %put \qed at margin in eqn

\topmargin -0.5in %default is 0.0, i.e. 1 inch from top.  Now it is 0.5
\textheight 9in
\oddsidemargin 0.25in %more reasonable margins!  (This is 1.25; can use <0 too)
\evensidemargin 0.25in %(same as above)
\textwidth 6.25in %(could use also -.5,-.5,7.5 for instance)
\tolerance=1600 %allow some tolerance in extending after line.
\parskip=3pt plus 1pt minus 1pt
\overfullrule=0pt %no dark lines if overfull

%%%% FROM HERE, macros special for paper
\def\mm{m}
\def\pp{p}
\def\proj{\pi _1} % projection on first group of coordinates
\def\x{u}  % vector argument
\def\SEC{(SEC)}
\def\SECo{(SEC$^0$)}
\def\INV{(INV)}
\def\Fp{{\cal F}_\pp}
\def\Fpm{{\cal F}_\pp^\mm}  %there is a 3,3 and a 2,1 in text

\def\Rm{\R^m}
\def\Rn{\R^n}
\newcommand{\ball}[2]{{\bb B}(#1,#2)}
\newcommand{\cball}[2]{\overline {\bb B} (#1,#2)}
\def\Fnm{{\cal F}_n^\mm}
\def\clos{{\rm clos\,}}
\def\inter{{\rm int\,}}

\def\CCF{compatible class of functions}
\begin{document}

\pagestyle{empty}
\begin{flushright}
{\sl CPS615-fall.94.hmk3}
\end{flushright}

\vspace{1cm}
\vspace{1cm}

\hrule
\hrule
\hrule
\hrule
\hrule
\hrule
\hrule
\hrule
\hrule
\hrule
\vspace{.05cm}
\hrule
\vspace{.25cm}
{\bf 
Fortran-90 programs for the solution of linear systems
}
\vspace{.25cm}
\hrule
\vspace{.05cm}
\hrule
\hrule
\hrule
\hrule
\hrule
\vspace{1cm}

\begin{center}
Saleh Elmohamed\\
Northeast Parallel Architectures Center,\\
Syracuse University\\
Syracuse, NY 13244-4100\\
saleh@npac.syr.edu\\
\end{center}
\newpage
\pagestyle{empty}
\tableofcontents

\newpage

\pagestyle{plain}
\pagestyle{myheadings}
\markright{\sl CPS615-fall.94.hmk3}
\section{Solution to problem [1]}
\setcounter{page}{1}

\subsection{The programs}
The actual code which consists of two parts:
\begin{itemize}
\item The main routine which is the driver of the subroutine
\item The Gaussian elimination subroutine which uses {\bf partial
pivoting} and {\bf back substitution} to a linear system of equations
in the form of matrices.
\end{itemize}
\begin{verbatim}
!
      program xgauss_elim
!     
! Driver for routine gauss_elim
!
! M_a(1...n,1...n) is the input matrix. M_b(1...n,1...m) 
! is the input containing the m right-hand side vectors. On output, 
! matrix M_a is replaced by its inverse matrix, and matrix M_b is replaced 
! by the corresponding set of solution vectors. 
! Matrix M_ai carries matrix M_a, and matrix M_x carries matrix M_b locally.
! Matrix M_u carries the result of matrix M_a * inverse-matrix M_a.
! Matrix M_t carries the result ?.....
! Let m also stands for the number of equations in the linear system.
!
      interface 
         subroutine gauss_elim(M_ai,n,NP,M_x,m,MP)
           integer :: MP, NP, n, m
           real, dimension (1:NP,1:NP) :: M_ai
           real, dimension (1:NP,1:MP) :: M_x
         end subroutine
      end interface
      integer, parameter :: MP = 20, NP = 20
      integer  :: j, k, l, m, n
      real, dimension (NP, NP) :: M_ai
      real, dimension (NP, MP) :: M_x
!
      real, dimension (NP, NP) :: M_a = 0.0
      real, dimension (NP, MP) :: M_b = 0.0
      real, dimension (NP, NP) :: M_u = 0.0      
      real, dimension (NP, MP) :: M_t = 0.0
      character*3 :: indicator
      open(7,file='matrix.dat.gauss',status='old')
10    read(7,'(a)') indicator
      if (indicator == 'END') then
         goto 99
      end if
      read(7,*)
      read(7,*) n,m
      read(7,*)
      read(7,*) ((M_a(k,l), l=1,n), k=1,n)
      read(7,*)
      read(7,*) ((M_b(k,l), k=1,n), l=1,m)
!
! Save matrices for later testing of results ...
!
      M_ai = M_a
      M_x = M_b
!
! Invert the matrix ...
! ----------------------
      call gauss_elim(M_ai,n,NP,M_x,m,MP)
! ----------------------
      print *, "Inverse of Matrix M_a : "
      do k=1,n,1
        write(*,'(1h ,(6f12.6))') (M_ai(k,l), l=1,n)
      end do
!
! Test Results and check Inverse.
! number of operations are O(n^3) ... If m = n then it in O(m^3) ...
!
      print *, "M_a multiplied by Inverse of M_a ... (unit matrix) ..."
      do  k=1,n,1
        do  l=1,n,1
           M_u(k,l)=0.0
           do  j=1,n,1
              M_u(k,l)=M_u(k,l)+M_a(k,j)*M_ai(j,l)
           end do
         end do
         write(*,'(1h ,(6f12.6))') (M_u(k,l), l=1,n)
      end do
!
! Check Vector Solutions ...
!
! Operation: below, we have 2n * n * m, or for m = n, 
! it is = O(m^3) for n >> m, it is O(n^2) ...
!
      write(*,*) 'x-vector is:'
      do k=1,m,1
       write(*,'(1x,6f12.6)') (M_x(l,k), l=1,n)
      end do
      print *, "Check the following vectors for equality:"
      write(*,'(t12,a10,t23,a18)') ' Original ','Matrix*Solution'
      do l=1,m,1
         write(*,'(1x,a,i2,a)') 'Vector ',l,':'
         do k=1,n,1
            M_t(k,l)=0.0
            do j=1,n,1
               M_t(k,l)=M_t(k,l)+M_a(k,j)*M_x(j,l)
            end do
            write(*,'(8x,2f12.6)') M_b(k,l),M_t(k,l)
         end do 
      end do
      print *, "***********************************"
      goto 10
99    close(7)
      end program xgauss_elim
!
! So, we have O(m^3) from gauss_elim and O(m^3) + (O(n^2) V O(m^3))
! yeilding an overall of O(m^3) operations where n is the number
! of the equations ...
!
|******************************************************************
      subroutine gauss_elim(M_a,n,np,M_b,m,mp)
!
! This routine solves a set of linear equations using Gaussian elimination
! method. M_a(1...n,1...n) is the input matrix. M_b(1...n,1...m)
! is the input containing the m right-hand side vectors. On output, 
! matrix M_a is replaced by its inverse matrix, and M_b is replaced by the 
! corresponding set of solution vectors.
!
      integer, parameter :: NMAX = 50
      integer  :: m, mp, n, np, i, icol, irow, j, k, l, ll
      real, dimension (1:np,1:np) :: M_a, M_b
      integer, dimension (1:NMAX) :: indexc = 0      !These three arrays are
      integer, dimension (1:NMAX) :: indexr = 0      !used for bookkeeping 
      integer, dimension (1:NMAX) :: ipivot = 0      !purposes for the pivot.
      real  :: large = 0.0
      real  :: temp = 0.0
      real  :: pivotinv = 0.0
!
! Main loop over the columns to be reduced.
!
       do i=1,n,1
         large = 0.0
!
! The outer loop of the search for a pivot element.
!
        do j=1,n,1
          if (ipivot(j) /= 1) then
             do k=1,n,1
                if (ipivot(k) == 0) then
                   if (abs(M_a(j,k)) >= large) then
                      large=abs(M_a(j,k))
                      irow=j
                      icol=k
                    end if
                    else if (ipivot(k) > 1) then
                          print *, "error: singular matrix in gauss_elim..."
                          stop
                 end if
             end do 
           end if
        end do
        ipivot(icol)=ipivot(icol)+1
!
! We now have the pivot element, so we interchange rows, if needed, to
! put the pivot element on the diagonal. The columns are not physically 
! interchanged, only relabeled:
! indexc(i), the column of the ith pivot element, is the ith column that
! is reduced, while indexr(i) is the row in which that pivot element was
! originally located. If indexr(i) != indexc(i) there is an implied column 
! interchange. With this form of bookkeeping, the solution b's will end 
! up in the correct order, and the inverse matrix will be scrambled by 
! columns.
!
        if (irow /= icol) then
          do l=1,n,1
             temp=M_a(irow,l)
             M_a(irow,l)=M_a(icol,l)
             M_a(icol,l)=temp
          end do
          do l=1,m,1
             temp=M_b(irow,l)
             M_b(irow,l)=M_b(icol,l)
             M_b(icol,l)=temp
          end do
        end if
!
! Below, we have 1 + n + m. If m == n, then 1+2m ===> O(m), else
! If m << n then the number is in O(n) ...
!
        indexr(i)=irow
        indexc(i)=icol
        if (M_a(icol,icol) == 0) then 
          print *, "singular matrix in gauss_elim ..."
          stop
        end if
        pivotinv=1.0/M_a(icol,icol)
        M_a(icol,icol)=1.0
        M_a(icol,1:n)=M_a(icol,1:n)*pivotinv
        M_b(icol,1:m)=M_b(icol,1:m)*pivotinv
!
! Below, we have the following:
! (2n + 2m) n = 2n^2 + 2mn, for m = n ===> 2m^2 + 2m^2 = 4 m^2 ==>
! O(m^2), but for m << n ===> O(n^2)
!
         do ll=1,n,1
           if (ll /= icol) then
             temp=M_a(ll,icol)
             M_a(ll,icol)=0.0
             M_a(ll,1:n)=M_a(ll,1:n)-M_a(icol,1:n)*temp
             M_b(ll,1:m)=M_b(ll,1:m)-M_b(icol,1:m)*temp
           end if
         end do
        end do
!
! Computing a solution vector from an augmented
! given matrix that has gone through lower decomposition,
! compute a solution vector from it.
!
        do l=n,1,-1
          if (indexr(l) /= indexc(l)) then
            do k=1,n,1
               temp=M_a(k,indexr(l))
               M_a(k,indexr(l))=M_a(k,indexc(l))
               M_a(k,indexc(l))=temp
            end do
          end if
        end do
        end subroutine gauss_elim
!
! So, if m == n, the overall number of operations is
! (1 + 2m + 4 m^2)*m which is in  O(m^3). 
!
\end{verbatim}

\subsection{The input system}

It consists of the following linear equations:
\be{eq.1}
x_1 + x_2 + x_3 = 2
\ee
\be{eq.2}
- x_1 + 2 x_2 - 3 x_3 = 32
\ee
\be{eq.3}
3 x_1 - 4 x_3 = 17
\ee

$A x = b$ where the coefficient matrix $A = [a_{ik}]$ is the $n x n$
matrix:

\vspace{0.1cm}
$\left(\begin{array}{ccc}1&1&1\\ -1&2&-3\\ 3&0&-4\end{array}\right)$

\vspace{0.1cm}
and $b$ is the RHS vector which can be written as also a matrix:

\vspace{0.1cm}
$\left(\begin{array}{ccc}2&32&17\end{array}\right)$

\subsection{The output values}

\begin{verbatim}

Inverse of Matrix M_a : 
     0.296296   -0.148148    0.185185
     0.481481    0.259259   -0.074074
     0.222222   -0.111111   -0.111111

M_a multiplied by Inverse of M_a ... (unit matrix) ...
     1.000000    0.000000    0.000000
     0.000000    1.000000    0.000000
     0.000000    0.000000    1.000000

x-vector (SOLUTION) is:
    -1.000000    8.000000   -5.000000

Check the following vectors for equality:
            Original     Matrix*Solution
 Vector  1:
            2.000000    2.000000
           32.000000   32.000000
           17.000000   17.000002
***********************************

\end{verbatim}

\section{Solution to problem [2]}
\subsection{The program}
The code consists of two parts:
\begin{itemize}
\item The main routine which is the driver of the subroutine.
\item The LU decomposition routine using the {\it Crout Algorithm}.
\end{itemize}
\begin{verbatim}
!*******************************************************************
      subroutine LUdecomp_sub(M_a,n,np,indexer,d)
!
! Given the matrix M_a(1...n,1...n) , this subroutine computes the LU
! decomp of a row-wise permutation of M_a. Then, the original M_a with
! this decomposition. So, first, M_a and n are input, then the output
! goes to M_a. indexer(1...n) is an integer vector (output) listing the
! row permutation effected by the partial pivoting. d is output as
! +/-1 depending on whether the number of row interchanges was even or
! odd (even +1, odd -1).
!
      integer, parameter :: NMAX = 500
      real, parameter :: very_small= 1.0e-20
      integer :: n, np, i, imaximum, j, k
      integer, dimension(n) :: indexer
      real, dimension(np,np) :: M_a
      real :: s = 0.0, amaximum = 0.0, temp = 0.0, sum = 0.0, d
      real, dimension(NMAX) :: V_v = 0.0

!
! floating point operations are performed inside the following 
! n loop:
!
      d = 1.0
      do i=1,n,1                 ! loop over rows ...
        amaximum=0.0
        do j=1,n,1
          if (abs(M_a(i,j)) > amaximum) then
           amaximum=abs(M_a(i,j))
          end if
        end do
        if (amaximum == 0.0) then
          print *, "error: singular matrix in LUdecomp ..."
          stop
        end if
        V_v(i)=1.0/amaximum      ! keep track of the scaling value ...
      end do
!
! Implementation of Crout's Algorithm:
! begin
!
! (1):
!
      do j=1,n,1
         do i=1,j-1,1
           sum=M_a(i,j)
           do k=1,i-1,1
             sum=sum-M_a(i,k)*M_a(k,j)
           end do
           M_a(i,j)=sum
         end do
!
! Operations performed above is in O(n^2) or 2(n-1)(n-1) = 
! 2n^2 - 4n + 2 ===> for n = m, that is O(m^2).
!
! (2):
!
        amaximum=0.0             ! next, search for the largest pivot ...
        do i=j,n,1
          sum=M_a(i,j)
          do k=1,j-1,1
            sum=sum-M_a(i,k)*M_a(k,j)
          end do
          M_a(i,j)=sum
          temp=V_v(i)*abs(sum)
          if (temp >= amaximum) then    ! Is the pivot number is better than
            imaximum=i                  ! the best we obtained so far? if so
            amaximum=temp               ! then do the update ...
          endif
        end do
!
! Operations done = (2n + 1) n which is on O(n^2) ...
!
! (3):
!
        if (j /= imaximum) then   ! If an interchange of rows is necessary
          do k=1,n,1              ! then it is done ...
            temp=M_a(imaximum,k)
            M_a(imaximum,k)=M_a(j,k)
            M_a(j,k)=temp
          end do
          d=-d                    ! also the parity of d is changed and
          V_v(imaximum)=V_v(j)    ! the scaling factor is interchanged ...
        endif
        indexer(j)=imaximum
!
! (4):
!
        if(M_a(j,j) == 0.0) then   ! If the pivot == 0.0 then the matrix
         M_a(j,j)=very_small       ! is singular.
        endif
        if(j /= n) then
          temp=1.0/M_a(j,j)
          M_a(j+1:n,j)=M_a(j+1:n,j)*temp
        endif                      ! Iterate over the next column in
      end do                       ! the reduction till the nth one ...
!
!  fp operations done above: 1 + n  ===> in O(n) ...
!  Now, the overall number of FP operations is
!  ((2n^2 - 4n +2) + n(2n+1)) n = n(2n^2 + 2n^2 -4n +n +2)
!  n(4n^2 -3n +2) = 4n^3 -3n^2 + 2n  ===> is in O(m^3) where
!  m is the number of equations in the linear system.
!
!  End Crout method.
!
      end subroutine LUdecomp_sub
\end{verbatim}

\subsection{Another LU decomposition}
Here, we uses  the LU decomposition method:

\begin{verbatim}
      program LU
!
!     M_a(1...n,1...n) is the input matrix. M_b(1...n,1...n) 
!     is the input containing the m right-hand side vectors. On output, 
!     matrix M_a is replaced by its inverse matrix, and matrix M_b is 
!     replaced by the corresponding set of solution vectors. 
!
      interface 
         subroutine LUdecomp_sub(M_c,n,NP,indexer,p)
           integer :: NP, n
           real :: p
           real, dimension (1:NP,1:NP) :: M_c
           integer, dimension (1:NP) :: indexer
         end subroutine
         subroutine LUbacksub_sub(M_c,n,NP,indexer,M_x)
           integer :: MP, NP, n, m
           real, dimension (1:NP,1:NP) :: M_c
           real, dimension (1:NP) :: M_x
         end subroutine
      end interface
      parameter (NP = 20)
      real p
      integer  j,k,l,m,n
      character*3 indicator
      real, dimension (NP,NP) :: M_a = 0.0
      real, dimension (NP,NP) :: M_b = 0.0
      real, dimension (NP,NP) :: M_c = 0.0
      real, dimension (NP) :: M_x = 0.0
      integer, dimension (NP) :: indexer = 0
      open(7,file='matrix.dat',status='old')
      read(7,*)
10    read(7,*)
      read(7,*) n,m
      read(7,*)
      read(7,*) M_a(1:n,1:n)
      read(7,*)
      read(7,*) ((M_b(k,l), k=1,n), l=1,m)
!
!     Save matrix a for later testing ...
!
      M_c = M_a
!
!     first: do LU decomposition ...
!--------------------------------------
      call LUdecomp_sub(M_c,n,NP,indexer,p)
!--------------------------------------
!     Solve equations for each right-hand vector ...
!
      do k=1,m,1
        M_x(1:n)=M_b(1:n,k)
!
!     second: do LU backsubstitution ...
!--------------------------------------
        call LUbacksub_sub(M_c,n,NP,indexer,M_x)
!--------------------------------------
!     test results with original matrix
!
        print *, "Right-hand side vector:"
        write(*,'(1x,6f12.6)') M_b(1:n,k)
        print *, "Result of matrix applied to solution vector"
        do l=1,n,1
          M_b(l,k)=0.0
          do j=1,n,1
            M_b(l,k)=M_b(l,k)+M_a(l,j)*M_x(j)
          end do
        end do
!
!  2n * n * m , for m << n, this is O(n*n), but for m = n
!  will end up in O(n^3) + O(n^2) (from LUbacksub) ===> O(n^3)
!
        write(*,'(1x,6f12.6)') M_b(1:n,k)
        print *, "***********************************"
      end do
      read(7,'(a3)') indicator
      if (indicator /= 'END') then
         goto 10
      end if
      close(7)
      end program LUbacksub

!
!*******************************************************************
      subroutine LUbacksub_sub(M_a,n,np,indexer,M_b)
! Given a set of n linear equations; Ax = b . Here M_a(1..n,1..n) is
! input, not as the matrix A but rather as its LU decomposition,
! determined by the LUdecomp_sub routine. indexer(1..n) is input as the
! permutation vector returned by LUdecomp_sub. M_b(1..n) is input as
! the RHS vector b, and returns with solution vector x. No
! modification is done to matrix M_a, n and the index (indexer)...
!
      integer, dimension (1:n) :: indexer
      integer n, np, i, ii, j, ll
      real, dimension (np, np) :: M_a
      real, dimension (n) :: M_b
      real sum = 0.0
      ii = 0
!
! 2n * n  ===> O(n^2), since m = n, O(m^2) ...
!
      do i=1,n,1
        ll=indexer(i)
        sum=M_b(ll)
        M_b(ll)=M_b(i)
        if (ii /= 0) then
          do j=ii,i-1,-1
            sum=sum-M_a(i,j)*M_b(j)
          end do
        else if (sum /= 0.0) then
          ii=i
        endif
        M_b(i)=sum
      end do
!
! (2n + 1) n  ===> O(n^2) ===> O(m^2) ...
!
      do i=n,1,-1                ! Do the backsubstitution ...
        sum=M_b(i)
        do j=i+1,n,1
          sum=sum-M_a(i,j)*M_b(j)
        end do
        M_b(i)=sum/M_a(i,i)      ! Part of the solution vector is
      end do                     ! is computed ...
      end subroutine LUbacksub_sub
\end{verbatim}

\subsection{The input system}

It consists of the following linear equations:
\be{eq.1}
3x_1 + 5x_2 + 2x_3 = 8
\ee
\be{eq.2}
8x_2 + 2x_3 = -7
\ee
\be{eq.3}
6x_1 + 2x_2 + 8x_3 = 26
\ee

$A x = b$ where the coefficient matrix $A = [a_{ik}]$ is the $n x n$
matrix:

\vspace{0.1cm}
$\left(\begin{array}{ccc}3&5&2\\ 0&8&2\\ 6&2&8\end{array}\right)$

\vspace{0.1cm}
and $b$ is the RHS vector which can be written as also a matrix:

\vspace{0.1cm}
$\left(\begin{array}{ccc}8&-7&26\end{array}\right)$

\subsection{The output values}
\begin{verbatim}
Right-hand side vector:
     8.000000   -7.000000   26.000000

x-vector (SOLUTION) is:
    -6.000000    1.791667    4.333333

Result of matrix applied to solution vector
     8.000000   -7.000000   26.250000
\end{verbatim}

\section{Number of floating-point operations}

In {\it gauss elim}, we have $1+n+m$. If $m == n$, then $1+2m$
 $\longrightarrow O(m)$, else If $m \ll n$ then the number is in
O(n). 
$(2n + 2m) n = 2n^{2} + 2mn$, for $m = n \longrightarrow  2m^{2} +
2m^{2} = 4 m^{2} \longrightarrow  O(m^2)$, but for $m \ll n
\longrightarrow O(n^2)$ ...

So, if $m == n$, the overall number of operations is $(1 + 2m + 4 m^2)m$
which is in  $O(m^3)$.

The total number (including those of the main program): 
$O(m^3)$ from {\it gauss elim} and $O(m^3) + (O(n^2) \bigvee O(m^3))$
yeilding an overall of $O(m^3)$ operations where m is the number
of the equations. 

In the LUdecomp driver we have $(2n * n * n)$ operations and since 
$m = n$, that is equivalent to $O(m^3)$. In the LUdecomp subroutine,
the overall number of FP operations is $((2n^{2} - 4n +2) + n(2n+1)) n
= n(2n^{2} + 2n^{2} -4n +n +2) = n(4n^{2} -3n +2) = 4n^{3} -3n^{2} + 2n
\longrightarrow O(m^3)$. So, there are in the order of $m^3$
operations in this system.

\section{References}
\begin{enumerate}
\item Gene Golub and Charles van Loan, {\it Matrix Computations}, Johns
Hopkins Univ. Press, 1989.

\item W. Press, S. Teukolsky, W. Vettrling, and B. Flannery, {\it
Numerical Recipes - The Art of Scientific Computing}, 2nd edition,
Cambridge Univ. Press, 1992.
\end{enumerate}

\end{document}

