\chapter{The Matrix Stack}

The matrix stack\index{MStack } lets you return VMatrix
references from other functions. It also keeps track of
intermediate calcuations stored in the heap, and deletes them
when necessary. The current version of YAMP uses a deep-copy
technique to push matrices faster. Assembly language programmers
should already be familiar with the method of using the stack
for storing registers or data when jumping to other sections of
code.  The method used here is little bit simpler than the
assembly stack methods. It also requires the least amount of
data structure overhead that I could come up with to solve the
problem of trapping allocations on the heap.  Besides, to quote
Knuth\cite{Kn:fa}, "A stack is the simplest kind of linked
list."

\section{How Does the Stack Work?}
The stack accomplishes three objectives:
\begin{itemize}
   \item to keep track of intermediate matrix references during
         nested function calls
   \item to keep track of the level of nesting in the
         nested function calls.
   \item explicit control of garbage collection for temporary
         matrices.
\end{itemize}
Consider the problems associated with the following code
\begin{verbatim}
   VMatrix &f1( void )
     {
       VMatrix b("local mat",5,5);
       return b;
     }
   VMatrix &f2( void )
     {
        VMatrix *b = new VMatrix("local mat",5,5);
        return *b;
     }
   main()
     {
        VMatrix c;
       
        c= f1() + f2();
     }
\end{verbatim}
This code won't work. The call to f1() attempts to return a
local matrix. However, C++ calls the destructor of b before
returning, so the returned reference points to an
abandonded location in the program stack, and hence to
garbage. The call to f2() does not quite solve the problem of
f1(). It returns a valid reference to a matrix, but the
reference is lost once the equal sign is processed in the
main function. It is lost in the sense that no variable
contains the reference, so it cannot be removed from the
heap. If you call f2() too often, the heap will fill up and
the program will crash. Both of these problems are solved
by the stack. The solution is to store a copy of the matrix
in the stack, which is a global data structure. The correct
program above then is 
\begin{verbatim}
   #include "virt.h"
   MStack *Dispatch = new MStack;

   VMatrix &f1( void )
     {
      VMatrix b("local mat",5,5);
      Dispatch->Push(b); 
      return Dispatch->ReturnMat();
     }
   VMatrix &f2( void )
     {
        VMatrix *b = new VMatrix("local mat",5,5);
        Dispatch->Push(*b);
        delete b;
        return Dispatch->ReturnMat();
     }
   main()
     {
        VMatrix c;
        c= f1() + f2();
        vclose();
     }
\end{verbatim}
This program will work as long as you do not need to
control the function nesting level. The next example shows
how you can use the function nesting level to use matrix
assignment in all functions.
\begin{verbatim}
   #include "virt.h"
   MStack *Dispatch = new MStack;

   VMatrix &f1( void )
     {
      VMatrix b = Fill(5,5,1);
      b.Nameit("Local");
      Dispatch->Push(b); 
      return Dispatch->ReturnMat();
     }
   VMatrix &f2( void )
     {
        Dispatch->Inclevel()
        VMatrix *b = new VMatrix;
        *b = Fill(5,5,2);
        b->Nameit("pointer");
        Dispatch->Push(*b);
        delete b;
        return Dispatch->Decreturn();
     }
   VMatrix &f3( void )
     {
         Dispatch->Inclevel();
         VMatrix b  = f1() + f2();
         Dispatch->Push(b);
         return Dispatch->DecReturn();
      }
   main()
     {
        VMatrix c = f1() + f3() + f2();
        vclose();
     }
\end{verbatim}
Note that a matrix assignment is used in f3(). The stack is
cleared of intermediate matrices when using matrix
assignment by calling CleanStack(). The matrix nesting level
index is used to determine which matrices are to be removed
from the stack. When f3() is called, f2() has already been
called in main(), so the stack has the result of f2() in
it. This copy will not be erased when f3() calls the equals
operator because the nesting level is incremented at the
function entry. Only the results on the stack created while
f3() is active are popped off the stack. Thus the first copy
of f2() in main() are preserved. 

The difference between Dispatch-$>$ReturnMat() and
Dispatch-$>$DecReturn() is that DecReturn decrements the matrix
nesting level and ReturnMat() does not. You use ReturnMat() if
you did not need matrix assignment in a function, and hence did
not call Inclevel() at the top. You use DecReturn() if you
needed to use matrix assignment hence called Inclevel() at the
top. 

I chose the word Dispatch for the stack root to suggest
that "messages" may be sent to the stack through the stack
dispatcher.  The stack should be declared in the global data
area as a pointer. This is expected throughout the program.  

\subsection{Garbage Collection using the Stack}
Some authors argue that an explicit stack is not necessary, and
is indeed the responsiblity of the compiler (See Eckel
\cite{EK:UC}). The problem with this approach is that using the
program stack to keep track of the temporary matrices allows the
temporary matrices to stay on the stack until they go out of
scope. They hang around and cause stack overflows in large
subroutines. This happens even if the stack only contains a
reference to a matrix instead of a copy. This problem was discussed
by Davies in his matrix program Newmat (currently available in
the Borland C++ library on Compuserve).

The advantage to having an explicit stack is that you can
control the deletion of the temporary matrices. The
\verb+Cleanstack()+ function can be used for garbage collection
of temporary matrices. It is called in matrix assignment, and
the copy constructor. It is also a public function, so it may
also be called explicitly. It's my opinion that the overhead of
manipulating an explicit stack outweighs the storage problems
and stack overflow problems of letting the program stack do it
for you. Besides, I couldn't figure out how to let the program do
it for me.

\section{The Stack}
The class MStack is derived from the VMatrix class. The
declaration is given by
\begin{verbatim}

class MStack: public VMatrix {

  private:

   int stackloc, level;
   MStack *next;
   VMatrix& Pop( void );

  public:

    MStack( void );
    void Push( VMatrix &ROp );

    void Inclevel();
    void Declevel();
    VMatrix &DecReturn( void ){
         Declevel();
         next->level = level;
         return * ((VMatrix *) this->next);
    }
    VMatrix &ReturnMat( void ){
         next->level = level;
         return * ((VMatrix *) this->next);
    }
    void Cleanstack( void );
    void PrintStack( void );
};

extern MStack *Dispatch;
\end{verbatim}
Note the last line declares Dispatch as a pointer to a
MStack object. It is also declared extern. The external
declaration is used in the header file VIRT.H so that all
modules will know that Dispatch exits, and will be created
in the module containing the main program. 

The MStack is a subclass of the VMatrix class, so it
inherits the members and data of the VMatrix class, and has
access to the protected part of the VMatrix class. The
relationship is not strictly hierachical because the matrix
class must call Dispatch to achieve some of its functions.
I call it a symbiotic relationship, for whatever that's
worth. The point to remember is that Dispatch must be
instantiated as a global data object for its parent
class to work.

\subsection{The Private Parts of MStack}
The private data of MStack are stackloc, level and next.
stackloc\index{MStack!stackloc} keeps track of how many
matrices have been pushed onto the stack. It is used as a
debugging tool, and as a loop counter in PrintStack().

The variable, level\index{MStack!level}, is the function nesting
level. When CleanStack() is called, it only pops matrices off
the stack where level is greater than or equal to
Dispatch-$>$level.  The rule of thumb for altering level is that
a function uses matrix assignment, i.e. VMatrix::operator=().
This allows each function to have its own local stack, but all
of the local stacks happen to be in the same linked list.  level
can be modified by the public member functions Declevel() and
Inclevel(), which should be used in pairs in functions.
DecReturn() also decrements level.

The pointer, next\index{MStack!next}, stores the address of
the next matrix in the stack. Thus, the stack is a singly
linked list. The last matrix in the stack stores NULL in
the value of next.  Also, the list can be traversed only by
member functions since next is a private element. The
reference stored in Dispatch-$>$next is the address returned
by DecReturn() and ReturnMat().

The private part of MStack also has one member function,
Pop()\index{MStack!Pop()}.  Pop() is only called by
CleanStack(), so it is not needed as a public function.
Pop() works like most other pop routines. It reassigns
Dispatch-$>$next, and deletes the matrix pointed to by the
original Dispatch-$>$next. Pop() also decrements stackloc.
Pop() stops the program if stackloc becomes zero or if
Dispatch-$>$next is NULL.

\subsection{The Public Parts of MStack}
All of the stack elements are allocated on the heap, so
they will not exaust the program stack space. They may
however exaust the heap if there are too many stack
matrices allocated. The constructor,
MStack()\index{MStack!constructor}, is called first by the creation of
Dispatch. The first call to the constructor sets level = 1.
A static int records the fact that Dispatch has
been created. All other calls to the constructor are
through the Push() function.\footnote{ I cannot see a reason for
directly allocating another MStack other than Dispatch,
since the matrix functions push their results onto
Dispatch. Maybe someone with a parallel processor
would want several independent stacks, but the rest of the
program would have to be rewritten. The vmm system would
need to be scrapped too. It's clear YAMP is a PC program}
Otherwise, the constructor initializes stackloc=0, level=0
and next = NULL. Push() fixes this. There is no destructor
for MStack objects. They are deleted by pop.

Push()\index{MStack!Push()} creates a MStack object on the heap,
then inserts it into the stack. However the matrix stored in
MStack is a reference to the base vector in the pushed matrix.
It also increments Dispatch-$>$stackloc and stores it in the new
MStack object. It also stores Dispatch-$>$level in the new
object.  Dispatch-$>$Push() should be called before calling
DecReturn() or ReturnMat(). 

Earlier versions of Push() had speed problems because Push()
copied its argument into the new stack object. The current
version uses a deep copy technique \index{MStack!deep copy} for
storing stack matrices. The deep copy method boosted the speed
of my programs by about 15 to 20 percent. Push() creates the
reference by a call to the VMatrix function NewReference().
NewReference() \index{VMatrix!NewReference()} increments the
variable \verb+nrefs+ in the base vector. In the in-ram version,
\verb+nrefs+ is incremented in the vdoub structure. In the
virtual version, \verb+nrefs+ is incremented in the hdr
structure. In either case, the base vector is not freed until 
\verb+nrefs+ is zero. However, the rest of the virtual matrix is deleted. 

The stack matrices share the base vector with the pushed matrix,
and hence you have indirect access to the stack matrix. This
access is different from that in Eckle \cite{EK:UC} where I
found the deep-copy technique.  Although I can't think of why
you would want this access, here is how it can be accomplished.

\begin{verbatim}
        .
        .
        .
        VMatrix B = Tran(X)*X;
        Dispatch->Push(B);
        B.M(1,1) = 2.0;
        return Dispatch->DecReturn();
  }
\end{verbatim}

Note since the stack top matrix shares the base vector with B,
the assignment writes 2.0 in both B(1,1) and the (1,1) element
of the stack top. Thus the matrix pushed onto the stack will be
different from the matrix that is returned by the function call.
The only place where this could be a problem is on the stack top
and none of the other matrices in the stack. This is because the
stack top is the matrix that is returned in function calls.  I
advise against using this trick. Again I say, {\em
Dispatch-$>$Push() should only be called immediately before calling
DecReturn() or ReturnMat()}.

Declevel()\index{MStack!Declevel()} decrements the value of
level. It is called by DecReturn(). It also stops
the program if the nesting level has been decremented to
zero. This error indicates that there is a mismatched pair
of calls to Inclevel() and DecReturn(), or Inclevel() and
Declevel(). 

Inclevel()\index{MStack!Inclevel()} increments the function
nesting level. Usually Inclevel() is called whenever the
function will use a matrix assignment, VMatrix::operator=().
DecReturn() calls Declevel(), so Inclevel() and DecReturn()
should be used in pairs.  Functions that use matrix assignment
or matrix copy constructors do not always return matrices. In
this case you should use a pair of Inclevel-Declevel calls to
control the stack.  

The code for DecReturn() and ReturnMat() is
given below by
\begin{verbatim}
    VMatrix &DecReturn( void ){
         Declevel();
         next->level = level;
         return * ((VMatrix *) this->next);
     }
    VMatrix &ReturnMat( void ){
         next->level = level;
         return * ((VMatrix *) this->next);
     }
\end{verbatim}
The two functions perform similar functions except that
DecReturn() calls Declevel(). They both are called from
Dispatch, so *this is actually Dispatch. The functions set level
in next to Dispatch-$>$level. You use ReturnMat() if you did not
need matrix assignment in a function, and hence did not call
Inclevel() at the top. You use DecReturn() if you needed to use
matrix assignment hence called Inclevel() at the top. The return
statment recasts Dispatch as a VMatrix pointer, and returns the
matrix reference stored in next.

Cleanstack()\index{MStack!Cleanstack()} is short:
\begin{verbatim}
  while( Dispatch->next->level >= Dispatch->level
         && Dispatch->next->next) Dispatch->Pop();
\end{verbatim}
It calls Pop() while the next stack matrix has a level at least
as large as Dispatch-$>$level. It also stops the popping when
there is only one matrix left on the stack which is useful when
there are problems in the program. This clears the local stack.
Recall that Cleanstack() is called by VMatrix::operator=().
The equals operator replaces its left operand by the right
before calling Cleanstack(). The left operand is not in the
stack, so it is not deleted by Pop(). 

The copy constructor\index{VMatrix!Copy constructor}
\index{MStack!VMatrix Copy Constructor}  calls
Cleanstack() so it also eliminates temporary matrices.
Therefore, you should call Inclevel() and Declevel() if you intend
to use the copy constructors.
 
You may wish to call Cleanstack() occasionally if you are in the
habit of using statements such as
\verb+(Tran(X)*X).DispalayMat();+. \footnote{By the way, just what is
this, a statement or a function call?} These
statement--function--calls leave matrices on the
stack which are removed automatically during a matrix
assignment, a copy constructor, or by an explicit call to
Cleanstack().

PrintStack()\index{MStack!PrintStack()} traverses the
stack. It reports the stackloc, level, r, c and name for
each matrix in the stack.