%% 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-fall96}
\end{flushright}

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

\hrule
\hrule
\hrule
\hrule
\hrule
\hrule
\hrule
\hrule
\hrule
\hrule
\vspace{.05cm}
\hrule
\vspace{.25cm}
{\bf 
HPF/Fortran-90 programs for the solution of a linear systems and
Laplace heat equation.
}
\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-fall96}
\section{Solution to problem [1]}
\setcounter{page}{1}

In the previous assignment, a Fortran-90 was written to solve the
Gaussian elimination problem. For this assignment we convert that code
HPF by the DISTRIBUTE compiler directives and running it using PGI
HPF compiler on the SP2. Comparsion of timing between the use of a
BLOCK distribution vs. CYCLIC one. These directives are clearly shown
and commented in the programs.

\subsection{The programs}
The actual code which consists of the following:
\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.
\item Gauss elimination routine divides the task and assign it to
three other subs. 
\end{itemize}
\begin{verbatim}
!
! ********* HPF: Using the BLOCK directives *********

!!! The below subroutine gauss_elimination solves the input linear equations
!!! via gaussian elimination with partial pivoting.
!!! Input: A set of coefficients of this linear system forming a
!!! matrix, and the number of equations of the system. This number is
!!! also the matrix size.
!!!
!!! The main program (gauss) ....

       program gauss

       interface
         subroutine gauss_elim(c,n)
           integer n
           real, dimension(1:n, 1:n+1) :: c
         end subroutine
       end interface

!!! Here, we will assume that 20 is the max number of equations that
!!! we will deal with ...
       
       parameter (n = 20)
       integer i
       real, dimension(n, n+1)::augmented
       real*8 time1, time2, elapsed_time
       character*50 filename

!HPF$ DISTRIBUTE augmented (BLOCK, *)

!!! request name of input data file ...
       print '("/ Data file> "$)'
       read *, filename
       open (10, file=filename, status = "OLD")

       print *, ' Number of equations in the system: ', n
!!! read rows of the augmented matrix ...
       do i = 1, n
         read(10,*) augmented(i,1:n+1)
       enddo
       print *, '==========================='
       print *, ' The original matrix size: ', n
       print *, '==========================='
       print *, ' Original augmented matrix is : '
       print *, '==========================='
       do i = 1, n
         print *, augmented(i, 1:n+1)
       end do

       time1 = timef()
       call gauss_elim(augmented(1:n, 1:n+1), n)
       time2 = timef()
       elapsed_time = time2 - time1

       print *, ' Solution Vector : ' 
       print *, '==========================='
       print *, augmented( 1:n, 1:n+1)
       print *, ' Elapsed time = ', elapsed_time

       end

!**********************************************************
!**********************************************************
       subroutine gauss_elim(matrix,n)
         integer n
         real, dimension (1:n, 1:n+1)  :: matrix

!HPF$  INHERIT MATRIX

!!! local var ...
         integer i

!!! Also find_and_check_the_pivot, check_for_singularity, and back_sub
!!! are local subroutines to gauss_elim ...

!!! Do for as many equations in the linear system ... 

         do i = 1, n
           call find_and_check_the_pivot(matrix, i, n)
           call check_for_singularity(matrix, i, n)
         enddo
         
         call back_sub(matrix, n)

       end 

!**********************************************************
! See Numerical Recipes book.
!**********************************************************
       subroutine find_and_check_the_pivot(a, i, nsize)
         integer i, nsize
         real, dimension (1:nsize, 1:nsize+1) :: a

!!! local variables ...
         real, dimension (1:nsize+1) :: a_row
         integer pindex
         integer, dimension (1) :: pindex_locations

!!! the pivot locations is ...
         pindex_locations = maxloc(abs(a(i:nsize,i)))
         pindex = pindex_locations(1) + i - 1

!!! adjust the pivot row ...
         if(pindex /= i) then
            a_row = a(i,1:nsize+1)
            a(i,1:nsize+1) = a(pindex,1:nsize+1)
            a(pindex,1:nsize+1) = a_row
         endif
       end subroutine

!*******************************************************
!*******************************************************
! For more info on matrices and matrix singularity, (see Numerical
! Recipes book).

       subroutine check_for_singularity(a, i, nsize)
         integer i, nsize
         real, dimension (1:nsize, 1:nsize+1) :: a

!!! local variables ...
         integer k

         if(abs(a(i,i)) < 1.0e-6) then
            write(*,*) 'error: matrix is singular ... bye bye.'
            stop
         endif
         a(i,i:nsize+1) = a(i,i:nsize+1)/a(i,i)
         do k = i+1, nsize
            a(k,i:nsize+1) = a(k,i:nsize+1) - a(k,i) * a(i,i:nsize+1)
         enddo
       end subroutine

!********************************************************
!********************************************************
! This routine will do simple back substitution ... (see Numerical
! Recipes book.) 
!
       subroutine back_sub(a, nsize)
         integer nsize
         real, dimension (1:nsize, 1:nsize+1) :: a

!!! local variables ...
         integer k
         real temp

         do k = nsize-1, 1,-1
            temp = sum(a(k,k+1:nsize) * a(k+1:nsize,nsize+1))
            a(k,nsize+1) = a(k,nsize+1) - temp
            a(k,k+1:nsize) = 0.0
         enddo
       end subroutine

\end{verbatim}

\subsection{The input system}

It consists of the following set of 20 equations:
\begin{verbatim}
3. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 6.
0. 2. 2. -.3 .1 .5 .5 .5 .5 0. 0. 0. 0. -.4 -.4 -.4 -.4 -.4 -.4 -.4 12.
0. 0. 5. 1. 4. 1. 3. 1. 3. 1. 3. 1. 3. 1. 3. 1. 3. 1. 3. 1. 7.
0. 1. 0. 9. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 6.
0. 0. .2 .3 -7 .5 .5 .5 .5 0. 0. 0. 0. -.4 -.4 -.4 -.4 -.4 -.4 -.4 46.
0. 0. 0. 1. .2 8. 3. 1. 3. 1. 3. 1. 3. 1. 3. 1. 3. 1. 3. 1. -75.
0. 1. 1. 1. 1. 0. -9. 1. 1. .2 1. 1. 1. 1. .7 1. 1. 1. -.5 1. 19.
0. 2. 2. -.3 0. 0. 0. 19. .5 0. 4. -9. 0. .4 3. 5. .13 -.4 2. -.4 -15.
0. 0. 5. 1. 4. 1. 3. 1. -42. 1. 3. 1. 17. 1. 3. .11 3. 1. 9. 0. 5.
1. 1. 1. 1. 1. 1. 1. 1. 1. 51. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. -13.
0. 2. 7. 0. .1 0. .5 .5 .5 0. -11. 8. -3. .1 -.4 0. -.4 13. -.2 -.4 22.
0. 0. 0. 1. -.3 1. 2. 0. .3 0. 0. 65. 3. 19. .7 0. 19. 1. 3. 0. 10.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. -29. 1. 5. 1. 2. 1. 1. 1. -14.
0. 0. 2. 9. .1 2. .5 .17 .5 0. 0. 0. 0. -47. -.2 6. 0. 18. 0. -.4 33.
0. 0. 0. 1. 0. 1. 0. 1. 0. 0. 3. 1. 0. 1. 59. 1. 0. 1. 3. 0. -7.
0. 1. 1. 0. 1. 1. 0. 1. 1. 1. 0. 1. 0. 1. 1. 17. 0. 1. 1. 0. 18.
0. 0. 2. -.1 .1 0. 0. .8 .2 0. 0. .9 0. -.4 0. -.4 33. 0. 0. 0. 12.
0. 0. 0. 0. 0. 0. 3. 1. 0. 1. 0. 1. 3. 0. 3. 1. 0. 17. 7. 0. 14.
1. 1. 1. 1. 1. 1. .3 1. 1. 2. 1. 1. 0. 1. 0. 1. 0. 1. 8. 1. -29.
0. 2. 1. 0 .1 .5 .1 .8 .5 0. 0. 0. 0. -.7 0. -.4 -.4 -.1 -.4 97. 11.
\end{verbatim}

\subsection{The output values}

On one of the SP2 nodes (merlin1, merlin2, $\ldots$), after updating
the environment (SETENV), the following will take care of the
compilation:
{\bf pghpf -O gauss-BLOCK.hpf gauss-BLOCK.pgi} then {gauss-BLOCK.pgi}
will give us the output which consists of the original matrix, the
solution vector and the elapsed time:

\subsubsection{Using BLOCK distribution}

\begin{verbatim}

saleh:merlin8[16] gauss_BLOCK.pgi
/ Data file> inputdata.v1
  Number of equations in the system:           20
 ===========================
  The original matrix size:           20
 ===========================
  Original augmented matrix is : 
 ===========================
    3.000000       1.000000       1.000000       1.000000       1.000000    
    1.000000       1.000000       1.000000       1.000000       1.000000    
    1.000000       1.000000       1.000000       1.000000       1.000000    
    1.000000       1.000000       1.000000       1.000000       1.000000    
    6.000000    
       0.0000000   2.000000       2.000000     -0.3000000      0.1000000    
   0.5000000      0.5000000      0.5000000      0.5000000          0.0000000
       0.0000000      0.0000000      0.0000000 -0.4000000     -0.4000000    
  -0.4000000     -0.4000000     -0.4000000     -0.4000000     -0.4000000    
    12.00000    
                    ... etc. ...
  Solution Vector : 

    2.552054       2.830099       7.536022       1.899023      -7.336004    
   -10.31382     -0.3766738      -3.037920     -0.4032587     -0.2779179    
    6.267160      0.1674865      0.4818895      0.7439715     -0.1189687    
    1.751745      3.6666505E-02   2.642349      -4.211722     6.3707352E-02

  Elapsed time =    55.70101737976074     
\end{verbatim}

\subsubsection{Using the CYCLIC distribution}

\begin{verbatim}
saleh:merlin8[18] gauss_CYCLIC.pgi
/ Data file> inputdata.v1
  Number of equations in the system:           20
 ===========================
  The original matrix size:           20
 ===========================
  Original augmented matrix is : 

    3.000000       1.000000       1.000000       1.000000       1.000000    
    1.000000       1.000000       1.000000       1.000000       1.000000    
    1.000000       1.000000       1.000000       1.000000       1.000000    
    1.000000       1.000000       1.000000       1.000000       1.000000    
    6.000000    
    0.0000000   2.000000       2.000000     -0.3000000      0.1000000    
   0.5000000      0.5000000      0.5000000      0.5000000      0.0000000
       0.0000000      0.0000000      0.0000000 -0.4000000     -0.4000000    
  -0.4000000     -0.4000000     -0.4000000     -0.4000000     -0.4000000    

                    .... etc. ....

  Solution Vector : 

    2.552054       2.830099       7.536022       1.899023      -7.336004    
   -10.31382     -0.3766738      -3.037920     -0.4032587     -0.2779179    
    6.267160      0.1674865      0.4818895      0.7439715     -0.1189687    
    1.751745      3.6666505E-02   2.642349      -4.211722     6.3707352E-02

  Elapsed time =    64.38326835632324     
\end{verbatim}

As shown above, the elapsed time of using the CYCLIC distribution is
higher than the time of using BLOCK distribution. They differ by
8.68225097656250.

\subsection{The second version}
Another version of the above code (using SPREAD constructs) but with
only one subroutine doing the entire task will yield the following
output for the same input/size matrix as above:

\subsubsection{Using the BOLCK distribution}
\begin{verbatim}

  Solution Vector : 

    2.552054       2.830098       7.536021       1.899023      -7.336003    
   -10.31382     -0.3766738      -3.037920     -0.4032587     -0.2779179    
    6.267161      0.1674865      0.4818895      0.7439714     -0.1189687    
    1.751745      3.6666509E-02   2.642349      -4.211722      6.3707352E-02

  Elapsed time =    35.05432605743408     

\end{verbatim}

\subsubsection{Using the CYCLIC distribution}
\begin{verbatim}


  Solution Vector : 

    2.552054       2.830098       7.536021       1.899023      -7.336003    
   -10.31382     -0.3766738      -3.037920     -0.4032587     -0.2779179    
    6.267161      0.1674865      0.4818895      0.7439714     -0.1189687    
    1.751745      3.6666509E-02   2.642349      -4.211722      6.3707352E-02

  Elapsed time =    38.29395771026611     

\end{verbatim}

Again, as shown, the elapsed time using CYCLIC is higher than that of
using the BLOCK distribution. The difference is 3.23963165283203 which
is less than half of the previous time difference.


\newpage
\section{Solution to problem [2]}
\subsection{Summary of the problem}
Given a metal (x,y) of dimensions $0 \le x \le 2$ and $0 \le y \le 1$,
the theoretical steady-state heat distribution is given by the
following solution:
\be{heat-dist}
U(x,y) = -4.7 + 0.2 x - 0.3 y + 0.09 (x^2 - y^2) + \ln{(\sqrt{(x +
0.2)^2 + (y + 0.2)^2})}
\ee

The ``Dirichlet boundary'' conditions for the metal may be calculated
using the theoretical solution $U(x,y)$, along the following
intervals: {\it top}: $0 \le x \le 2$, $y = 1$, {\it bottom}: $0 \le x
\le 2$, $y = 0$, {\it right}: $x = 2$ , $0 \le y \le 1$, and {\it
left}: $x = 0$ , $0 \le y \le 1$.
\subsection{What need to be done?}
\begin{itemize}
\item [A] Solve Laplace heat equation using a mesh size of 0.10 and an
initial value of 0.0.
\item [B] Solve Laplace heat equation using a mesh size of 0.05 and
initial value of 0.0.
\item [C] Compare the number of iterations necessary to solve the Laplace
heat equation using the two mesh sizes above.
\item [D] Solve the Laplace heat equation using a mesh size of 0.10 and
initial internal value composed of the average of the ``Dirichelt
boundary values''.
\item [E] Compare the number of iterations necessary to solve the Laplace
heat equation using the initial values specified in the above [A] and
[D] items.
\item [F] A tolerance of 0.000003 should be used in the Jacobi iteration
algorithm.
\end{itemize}

\subsection{Solutions to Problem 2}

What is done for this part is to partition the mesh into intervals of
size $0.1$ starting from the point $(0,0)$ going along the x-axis to
point $(2,0)$. Then from $(0,0)$ going up along the y-axis to point
$(0,1)$. Given this mesh partition, the solution of U is determined. 

\subsection{The program}

Here is a listing for the program code:

\begin{verbatim}

       program laplace_heat_equation

       interface
          subroutine compute_points_xANDy(x_list,y_list,int_x,int_y)
             integer int_x,int_y
             double precision, dimension(1:int_x) :: x_list
             double precision, dimension(1:int_y) :: y_list
          end subroutine

          subroutine avg_via_Dirichlet(U,int_x,int_y)
             integer int_x, int_y
             double precision, dimension(1:int_x,1:int_y)  :: U
          end subroutine

          subroutine one_iteration_ala_jacobi(U,Uold,int_x,int_y)
             integer int_x, int_y
             double precision, dimension(1:int_x,1:int_y) :: U, Uold
          end subroutine

          subroutine assign_boundary_value(U,int_x,int_y,x_list,y_list)
             double precision, dimension(1:int_x,1:int_y) :: U
             integer int_x, int_y
             double precision, dimension(1:int_x) :: x_list
             double precision, dimension(1:int_y) :: y_list
          end subroutine

          subroutine analytic_error(error,U,int_x,int_y,x_list,y_list)
             double precision error
             integer int_x, int_y
             double precision, dimension(1:int_x, 1:int_y):: U
             double precision, dimension(1:int_x) :: x_list
             double precision, dimension(1:int_y) :: y_list
          end subroutine
       end interface

!!! The parameters used to set the mesh size, 0.1 or 0.05 ...
!!!
       parameter(int_x = 21)     !! = 41 for 0.05 mesh spacing.
       parameter(int_y = 11)     !! = 21 for 0.05 mesh spacing.
       integer counter
       double precision error
       double precision, dimension(1:int_x, 1:int_y):: U, Uold
       double precision, dimension(1:int_x) :: x_list
       double precision, dimension(1:int_y) :: y_list

!!! Compute the point (x,y) ...
       call compute_points_xANDy(x_list,y_list,int_x,int_y)

!!! Set (U, Uold) to (0.0, 1.0) ...
       U = 0.0d0                     !! "d" for double precision.     
       Uold = 1.0d0

!!! Put boundary values on the boundaries ... This is quite crucial
!!! for computing the other internal points of the grid/mesh ...
!!!
       call assign_boundary_value(U,int_x,int_y,x_list,y_list)

!!!********* the following call is used when solving 2.2a *****
!!!... otherwise it can be commented out ....
!!!
       call avg_via_Dirichlet(U,int_x,int_y)

       error = 1.0d0
       counter = 0
!!!
!!! As it is specified in the question, a tolerance of 0.000003 (or
!!! 3.0d-06) should be used in the Jacobi iteration algorithm.
!!!
       do while( error > 3.0d-06 )
          counter = counter + 1
          call one_iteration_ala_jacobi(U,Uold,int_x,int_y)
          error = maxval(abs(U - Uold))
!!! write(*,*) 'error norm between U & Uold after jacobi step ',
!!! counter,' is ',error 
       enddo

       write(*,*) 'Program ... Done compiling!!!'

!!! The use of "d" and dble for double precision computations ...
!!!
       write(*,*) 'The mesh resolution is ...'
       write(*,*) 'x-interval: ',2.0d0/dble(int_x - 1)
       write(*,*) 'y-interval: ',1.0d0/dble(int_y - 1)
!!
       write(*,*) 'Total number of Jacobi iteration steps: ', counter
       write(*,*) '   '
       write(*,*) 'The error norm (U, Uold) for last iteration step:', error
!!
       call analytic_error(error,U,int_x,int_y,x_list,y_list)
!!
       write(*,*) '   '
       write(*,*) 'The error norm (U, analytical_solution) is:', error
       end

!*************** avg_of_internal_points_via_Dirichlet ******
!
! The internal points of the grid is calculated based on the boundary
! points. Specifically, these points are the average of those grid
! boundary points.
!!!
       subroutine avg_via_Dirichlet(U,int_x,int_y)
          integer int_x,int_y
          double precision, dimension(1:int_x,1:int_y)  :: U

!!! local var ...
!!!
       double precision average
       integer k, j, istarting_point

       average = 0.0d0
       istarting_point = 1
       do k = istarting_point, int_x
          average = average + U(k,1) + U(k,int_y)
       enddo
       do k = istarting_point, int_y
          average = average + U(1,k) + U(int_x,k)
       enddo
       average = average / dble(2 * (int_x + int_y))
       do k = istarting_point + 1, int_x - 1
        do j = istarting_point + 1, int_y - 1
          U(k,j) = average
        enddo
       enddo

       end

!*********************** analytic_error ********************
! Returns the norm between the input grid U and the already known
! analytic solution.
!!!
       subroutine analytic_error(error,U,int_x,int_y,x_list,y_list)
          double precision error
          integer int_x,int_y
          double precision, dimension(1:int_x, 1:int_y):: U
          double precision, dimension(1:int_x) :: x_list
          double precision, dimension(1:int_y) :: y_list

!!! local variables ...
       double precision, dimension(1:int_x,1:int_y) :: Utemp
       integer k, j

       do k=1, int_x
         do j = 1, int_y
           call compute_val_of_U( Utemp(k,j),x_list(k),y_list(j))
         enddo
       enddo

!!! determine the error between the given U and the just derived one.
       error = maxval(abs(U - Utemp))

       end

!*************** one-iteration-ala-jacobi ******************
!
! This sub basically performs one iteration on the input grid U. The
! new values are returned in U. Upon termination, Uold holds the input U.
!!!
       subroutine one_iteration_ala_jacobi(U,Uold,int_x,int_y)
          integer int_x,int_y
          double precision, dimension(1:int_x,1:int_y) :: U,Uold

!!! local vars ...
       integer i, j, istarting_point
       double precision avg_point

       Uold = U
       istarting_point = 1
       avg_point = 4.0d0

       do i = istarting_point + 1, int_x - 1
         do j = istarting_point + 1, int_y - 1
          U(i,j)=(Uold(i-1,j)+Uold(i+1,j)+Uold(i,j-1)+Uold(i,j+1))/avg_point
         enddo
       enddo
       
       end

!***************** assign_boundary_value **********************
!
! assigns the boundary values on the boundary of the mesh
!!!
       subroutine assign_boundary_value(U,int_x,int_y,x_list,y_list)
          interface
             subroutine compute_val_of_U( uval,x,y)
                double precision uval,x,y
             end subroutine
          end interface

          double precision, dimension(1:int_x,1:int_y) :: U
          integer int_x,int_y
          double precision, dimension(1:int_x) :: x_list
          double precision, dimension(1:int_y) :: y_list
!!! local vars ...
       integer k
       double precision x, y

       do k=1,int_x
          y = 0.0d0
          call compute_val_of_U( U(k,1) ,x_list(k), y)
          y = 1.0d0
          call compute_val_of_U( U(k,int_y), x_list(k), y)
       enddo
       do k=1,int_y
          x = 0.0d0
          call compute_val_of_U( U(1,k), x, y_list(k))
          x = 2.0d0
          call compute_val_of_U( U(int_x,k), x, y_list(k))
       enddo

       end
!******************* compute_val_of_U **********************
!
! This routine calculates the Uval given the x and y points ...
!!!
       subroutine compute_val_of_U( Uval,x,y)
          double precision Uval,x,y
!!! local vars ...
       double precision part1

       part1 = -4.7 + 0.2 * x - 0.3 * y + 0.09 * (x * x - y * y)
       Uval = part1 + log(sqrt( ((x+0.2)**2.0) + ((y+0.2)**2.0) ) )
       write(*,*) 'U =',Uval,'x =',x,'y =', y

       end
!****************** compute_points_xANDy *******************
!
! Computes the set of x and y points along the x-axis and y-axis,
! respectively. 
!!!
       subroutine compute_points_xANDy(x_list,y_list,int_x,int_y)
          integer int_x,int_y
          double precision, dimension(1:int_x) :: x_list
          double precision, dimension(1:int_y) :: y_list

!!! local vars ...
       integer k

       do k = 1,int_x
          x_list(k) = 2.0 * (dble(k-1)/dble(int_x - 1))
       enddo
       do k = 1,int_y
          y_list(k) = 1.0 * (dble(k-1)/dble(int_y - 1))
       enddo
       
       end
\end{verbatim}

\subsection{The Program's Output for problem (2.1a)}

The program was run on of the SP2 nodes (merlin8) using the local F90
compiler. Here is a script of the actual run:

\begin{verbatim}
saleh:merlin8[43] f90 -o merlin8.laplace.0.1-no-avg equation.laplace.f90
** laplace_equation   === End of Compilation 1 ===
** avg_via_dirichlet   === End of Compilation 2 ===
** analytic_error   === End of Compilation 3 ===
** one_iteration_ala_jacobi   === End of Compilation 4 ===
** assign_boundary_value   === End of Compilation 5 ===
** compute_val_of_u   === End of Compilation 6 ===
** compute_points_xandy   === End of Compilation 7 ===
1501-510  Compilation successful for file equation.laplace.f90.

saleh:merlin8[44] merlin8.laplace.0.1-no-avg

 U = -5.96286411651810333 x = 0.000000000000000000E+00 y = 0.000000000000000000E+00
 U = -4.89397877805519421 x = 0.000000000000000000E+00 y = 1.00000000000000000
 U = -5.69921021173219611 x = 0.100000000000000006 y = 0.000000000000000000E+00
 U = -4.85646595380458912 x = 0.100000000000000006 y = 1.00000000000000000
 U = -5.46111875580239214 x = 0.200000000000000011 y = 0.000000000000000000E+00
 U = -4.81139800642014581 x = 0.200000000000000011 y = 1.00000000000000000
 U = -5.25083697885634670 x = 0.299999999999999989 y = 0.000000000000000000E+00
 U = -4.75953555608105194 x = 0.299999999999999989 y = 1.00000000000000000
 U = -5.06374516747745229 x = 0.400000000000000022 y = 0.000000000000000000E+00
 U = -4.70170648756675558 x = 0.400000000000000022 y = 1.00000000000000000

... etc...

 U = -3.32827202488223772 x = 2.00000000000000000 y = 0.699999999999999956
 U = -4.97798946107532814 x = 0.000000000000000000E+00 y = 0.800000000000000044
 U = -3.35523440077204782 x = 2.00000000000000000 y = 0.800000000000000044
 U = -4.93132804413421244 x = 0.000000000000000000E+00 y = 0.900000000000000022
 U = -3.38287066497822320 x = 2.00000000000000000 y = 0.900000000000000022
 U = -4.89397877805519421 x = 0.000000000000000000E+00 y = 1.00000000000000000
 U = -3.41131481264321001 x = 2.00000000000000000 y = 1.00000000000000000

 error norm between U & Uold after jacobi step  1  is  2.86255510625650889
 error norm between U & Uold after jacobi step  2  is  1.05420358806227865
 error norm between U & Uold after jacobi step  3  is  0.691190449620399994
 error norm between U & Uold after jacobi step  4  is  0.466032699785735804
 error norm between U & Uold after jacobi step  5  is  0.426660661717004119
 error norm between U & Uold after jacobi step  6  is  0.323497777440174827
 error norm between U & Uold after jacobi step  7  is  0.277283551707262710
 error norm between U & Uold after jacobi step  8  is  0.241818005008382908
... etc ...

 error norm between U & Uold after jacobi step  354  is  0.359948967165735212E-05
 error norm between U & Uold after jacobi step  355  is  0.344652620398733234E-05
 error norm between U & Uold after jacobi step  356  is  0.338237891650550182E-05
 error norm between U & Uold after jacobi step  357  is  0.323863336770813248E-05
 error norm between U & Uold after jacobi step  358  is  0.317836364960299989E-05
 error norm between U & Uold after jacobi step  359  is  0.304328082556537538E-05
 error norm between U & Uold after jacobi step  360  is  0.298665399078856808E-05

 The mesh resolution is ...
 x-interval:  0.100000000000000006
 y-interval:  0.100000000000000006
 Total number of Jacobi iteration steps:  360

 The error norm (U, Uold) for last iteration step: 0.298665399078856808E-05

 U = -5.96286411651810333 x = 0.000000000000000000E+00 y = 0.000000000000000000E+00
 U = -5.75101021329383855 x = 0.000000000000000000E+00 y = 0.100000000000000006
 U = -5.56831875906872753 x = 0.000000000000000000E+00 y = 0.200000000000000011
 U = -5.41703698397042466 x = 0.000000000000000000E+00 y = 0.299999999999999989
 U = -5.29254517458232598 x = 0.000000000000000000E+00 y = 0.400000000000000022
 U = -5.18993894727688421 x = 0.000000000000000000E+00 y = 0.500000000000000000
 U = -5.10523105372845176 x = 0.000000000000000000E+00 y = 0.599999999999999978
 U = -5.03535928025427992 x = 0.000000000000000000E+00 y = 0.699999999999999956
 U = -4.97798946107532814 x = 0.000000000000000000E+00 y = 0.800000000000000044
 U = -4.93132804413421244 x = 0.000000000000000000E+00 y = 0.900000000000000022
 U = -4.89397877805519421 x = 0.000000000000000000E+00 y = 1.00000000000000000
 U = -5.69921021173219611 x = 0.100000000000000006 y = 0.000000000000000000E+00

... etc ...

 U = -3.44857525438098467 x = 1.89999999999999991 y = 0.800000000000000044
 U = -3.47483397066822253 x = 1.89999999999999991 y = 0.900000000000000022
 U = -3.50187897388649860 x = 1.89999999999999991 y = 1.00000000000000000
 U = -3.14742717760134161 x = 2.00000000000000000 y = 0.000000000000000000E+00
 U = -3.17323029468883977 x = 2.00000000000000000 y = 0.100000000000000006
 U = -3.19888083376002363 x = 2.00000000000000000 y = 0.200000000000000011
 U = -3.22446087603561571 x = 2.00000000000000000 y = 0.299999999999999989
 U = -3.25007047994170239 x = 2.00000000000000000 y = 0.400000000000000022
 U = -3.27582417514370050 x = 2.00000000000000000 y = 0.500000000000000000
 U = -3.30184724532809870 x = 2.00000000000000000 y = 0.599999999999999978
 U = -3.32827202488223772 x = 2.00000000000000000 y = 0.699999999999999956
 U = -3.35523440077204782 x = 2.00000000000000000 y = 0.800000000000000044
 U = -3.38287066497822320 x = 2.00000000000000000 y = 0.900000000000000022
 U = -3.41131481264321001 x = 2.00000000000000000 y = 1.00000000000000000
    
 The error norm (U, analytical_solution) is: 0.163660439956903048E-02
\end{verbatim}

\subsection{The Program's Output for problem (2.1b)}

\begin{verbatim}

saleh:merlin8[65] f90 -o merlin8.laplace.0.05.no_avg equation.laplace.f90
** laplace_equation   === End of Compilation 1 ===
** avg_via_dirichlet   === End of Compilation 2 ===
** analytic_error   === End of Compilation 3 ===
** one_iteration_ala_jacobi   === End of Compilation 4 ===
** assign_boundary_value   === End of Compilation 5 ===
** compute_val_of_u   === End of Compilation 6 ===
** compute_points_xandy   === End of Compilation 7 ===
1501-510  Compilation successful for file equation.laplace.f90.

saleh:merlin8[66] merlin8.laplace.0.05.no_avg

U = -5.96286411651810333 x = 0.000000000000000000E+00 y = 0.000000000000000000E+00
 U = -4.89397877805519421 x = 0.000000000000000000E+00 y = 1.00000000000000000
 U = -5.82872103622507520 x = 0.500000000000000028E-01 y = 0.000000000000000000E+00
 U = -4.87620962801468227 x = 0.500000000000000028E-01 y = 1.00000000000000000
 U = -5.69921021173219611 x = 0.100000000000000006 y = 0.000000000000000000E+00
 U = -4.85646595380458912 x = 0.100000000000000006 y = 1.00000000000000000
 U = -5.57651343725686832 x = 0.149999999999999994 y = 0.000000000000000000E+00
 U = -4.83483126996424328 x = 0.149999999999999994 y = 1.00000000000000000
 U = -5.46111875580239214 x = 0.200000000000000011 y = 0.000000000000000000E+00
 U = -4.81139800642014581 x = 0.200000000000000011 y = 1.00000000000000000

 ... etc ...

 U = -4.93132804413421244 x = 0.000000000000000000E+00 y = 0.900000000000000022
 U = -3.38287066497822320 x = 2.00000000000000000 y = 0.900000000000000022
 U = -4.91156419708708913 x = 0.000000000000000000E+00 y = 0.949999999999999956
 U = -3.39698353224754168 x = 2.00000000000000000 y = 0.949999999999999956
 U = -4.89397877805519421 x = 0.000000000000000000E+00 y = 1.00000000000000000
 U = -3.41131481264321001 x = 2.00000000000000000 y = 1.00000000000000000
 error norm between U & Uold after jacobi step  1  is  2.92072301830327241
 error norm between U & Uold after jacobi step  2  is  1.08365346955412889
 error norm between U & Uold after jacobi step  3  is  0.716092779691341663
 error norm between U & Uold after jacobi step  4  is  0.487070609159459433
 error norm between U & Uold after jacobi step  5  is  0.444272288479668376
 error norm between U & Uold after jacobi step  6  is  0.337855956469814700
 error norm between U & Uold after jacobi step  7  is  0.289850475748066216
 error norm between U & Uold after jacobi step  8  is  0.252233827068212069
 error norm between U & Uold after jacobi step  9  is  0.228889384328354373

 ... etc ...

 error norm between U & Uold after jacobi step  1261  is  0.314099208775786565E-05
 error norm between U & Uold after jacobi step  1262  is  0.310746442711007376E-05
 error norm between U & Uold after jacobi step  1263  is  0.309282473498484478E-05
 error norm between U & Uold after jacobi step  1264  is  0.305980887294055037E-05
 error norm between U & Uold after jacobi step  1265  is  0.304539603135367543E-05
 error norm between U & Uold after jacobi step  1266  is  0.301288417947631615E-05
 error norm between U & Uold after jacobi step  1267  is  0.299869465170132798E-05

 The mesh resolution is ...
 x-interval:  0.500000000000000028E-01
 y-interval:  0.500000000000000028E-01
 Total number of Jacobi iteration steps:  1267

 The error norm (U, Uold) for last iteration step:
0.299869465170132798E-05

 U = -5.96286411651810333 x = 0.000000000000000000E+00 y = 0.000000000000000000E+00
 U = -5.85417103698801444 x = 0.000000000000000000E+00 y = 0.500000000000000028E-01
 U = -5.75101021329383855 x = 0.000000000000000000E+00 y = 0.100000000000000006
 U = -5.65556343965297437 x = 0.000000000000000000E+00 y = 0.149999999999999994
 U = -5.56831875906872753 x = 0.000000000000000000E+00 y = 0.200000000000000011
 U = -5.48900158878293443 x = 0.000000000000000000E+00 y = 0.250000000000000000

... etc ...

 U = -3.32827202488223772 x = 2.00000000000000000 y = 0.699999999999999956
 U = -3.34167747052657171 x = 2.00000000000000000 y = 0.750000000000000000
 U = -3.35523440077204782 x = 2.00000000000000000 y = 0.800000000000000044
 U = -3.36895983996019011 x = 2.00000000000000000 y = 0.849999999999999978
 U = -3.38287066497822320 x = 2.00000000000000000 y = 0.900000000000000022
 U = -3.39698353224754168 x = 2.00000000000000000 y = 0.949999999999999956
 U = -3.41131481264321001 x = 2.00000000000000000 y = 1.00000000000000000
    
 The error norm (U, analytical_solution) is: 0.510305779594766307E-03

\end{verbatim}


\subsection{Comparsion of the above two results (2.1c)}

For the case of x-interval:  $0.1$ and  y-interval: $0.1$ :

\begin{itemize}
\item Total number of Jacobi iteration steps is 360, and the error norm
between current $U$ and $Uold$ for the last iteration step is
$0.298665399078856808E-05$. While the error norm between $U$ and the
analytical solution is $0.163660439956903048E-02$. 
\end{itemize}

For the second case of mesh resolution equivalent to x-interval:
$0.05$  and y-interval: $0.05$:

\begin{itemize}
\item Total number of Jacobi iteration steps is 1267, and the error
norm between current $U$ and $Uold$ for the last iteration step is 
$0.299869465170132798E-05$. While the error norm between $U$ and the
analytical solution is  $0.510305779594766307E-03$.
\end{itemize}

As we can see above, the increase in the number of x-interval and
y-interval points (or as the resolution becomes finer and finer), the
Jacobi iterations increase about 4 folds. While the error norm between
the value of $U$ (the theoretical) and the analytical solution
decreases. The ratio of $0.510305779594766307E-03$ to
$0.163660439956903048E-02$ is about $0.31180725$. So the error
reduction is close to $70\%$.

\subsection{The Program's Output (2.2a)}

Using a mesh resolution of size 0.1 and the internal values composed
of averaging over the Dirichlet boundary values we obtain:

\begin{verbatim}
saleh:merlin8[71] f90 -o merlin8.laplace.0.1-w-avg equation.laplace.f90
** laplace_equation   === End of Compilation 1 ===
** avg_via_dirichlet   === End of Compilation 2 ===
** analytic_error   === End of Compilation 3 ===
** one_iteration_ala_jacobi   === End of Compilation 4 ===
** assign_boundary_value   === End of Compilation 5 ===
** compute_val_of_u   === End of Compilation 6 ===
** compute_points_xandy   === End of Compilation 7 ===
1501-510  Compilation successful for file equation.laplace.f90.

saleh:merlin8[75] merlin8.laplace.0.1-w-avg

 U = -5.96286411651810333 x = 0.000000000000000000E+00 y = 0.000000000000000000E+00
 U = -4.89397877805519421 x = 0.000000000000000000E+00 y = 1.00000000000000000
 U = -5.69921021173219611 x = 0.100000000000000006 y = 0.000000000000000000E+00
 U = -4.85646595380458912 x = 0.100000000000000006 y = 1.00000000000000000
 U = -5.46111875580239214 x = 0.200000000000000011 y = 0.000000000000000000E+00
 U = -4.81139800642014581 x = 0.200000000000000011 y = 1.00000000000000000
 U = -5.25083697885634670 x = 0.299999999999999989 y = 0.000000000000000000E+00
 U = -4.75953555608105194 x = 0.299999999999999989 y = 1.00000000000000000
 U = -5.06374516747745229 x = 0.400000000000000022 y = 0.000000000000000000E+00
 U = -4.70170648756675558 x = 0.400000000000000022 y = 1.00000000000000000
 U = -4.89493893803816427 x = 0.500000000000000000 y = 0.000000000000000000E+00
 U = -4.63873981798585433 x = 0.500000000000000000 y = 1.00000000000000000
 U = -4.74043104221283507 x = 0.599999999999999978 y = 0.000000000000000000E+00

... etc ...

 U = -4.93132804413421244 x = 0.000000000000000000E+00 y = 0.900000000000000022
 U = -3.38287066497822320 x = 2.00000000000000000 y = 0.900000000000000022
 U = -4.89397877805519421 x = 0.000000000000000000E+00 y = 1.00000000000000000
 U = -3.41131481264321001 x = 2.00000000000000000 y = 1.00000000000000000
 error norm between U & Uold after jacobi step  1  is  0.723135224902835994
 error norm between U & Uold after jacobi step  2  is  0.251921132554651095
 error norm between U & Uold after jacobi step  3  is  0.156335479281981549

... etc ...

 error norm between U & Uold after jacobi step  202  is  0.314593155037456995E-05
 error norm between U & Uold after jacobi step  203  is  0.305137468759042463E-05
 error norm between U & Uold after jacobi step  204  is  0.290039736983160878E-05

 x-interval:  0.100000000000000006
 y-interval:  0.100000000000000006
 Total number of Jacobi iteration steps:  204
    
 The error norm (U, Uold) for last iteration step: 0.290039736983160878E-05

... etc ...

 U = -3.14742717760134161 x = 2.00000000000000000 y = 0.000000000000000000E+00
 U = -3.17323029468883977 x = 2.00000000000000000 y = 0.100000000000000006
 U = -3.19888083376002363 x = 2.00000000000000000 y = 0.200000000000000011
 U = -3.22446087603561571 x = 2.00000000000000000 y = 0.299999999999999989
 U = -3.25007047994170239 x = 2.00000000000000000 y = 0.400000000000000022
 U = -3.27582417514370050 x = 2.00000000000000000 y = 0.500000000000000000
 U = -3.30184724532809870 x = 2.00000000000000000 y = 0.599999999999999978
 U = -3.32827202488223772 x = 2.00000000000000000 y = 0.699999999999999956
 U = -3.35523440077204782 x = 2.00000000000000000 y = 0.800000000000000044
 U = -3.38287066497822320 x = 2.00000000000000000 y = 0.900000000000000022
 U = -3.41131481264321001 x = 2.00000000000000000 y = 1.00000000000000000
    
 The error norm (U, analytical_solution) is: 0.161867772884427552E-02
\end{verbatim}


\subsection{Comparsion (2.2b)}

\begin{enumerate}
\item Using a mesh size of $0.1$ and no averaging to compute the
internal values of the mesh, we obtain the following:

\begin{verbatim}
 x-interval:    0.100000000000000006
 y-interval:    0.100000000000000006

 Total number of Jacobi iteration steps:  
                
                360

 The error norm between U and Uold for last iteration step: 
                
                0.298665399078856808E-05
    
 The error norm between $U$ and the analytical solution is: 
                
                0.163660439956903048E-02
\end{verbatim}

\item Using a mesh of also size $0.1$ and computing the internal
values by averaging the {\it Dirichlet boundary} values, we obtain:
\begin{verbatim}

 x-interval:   0.100000000000000006
 y-interval:   0.100000000000000006

 Total number of Jacobi iteration steps:  

               204
    
 The error norm between U and Uold for last iteration step: 

               0.290039736983160878E-05

 The error norm between U and the analytical solution is: 

               0.161867772884427552E-02
\end{verbatim}
\end{enumerate}
In summary, ss shown above, with averaging we obtain fewer Jacobi
iterations and and overall smaller error.

\end{document}

