\chapter{Virtual Matrices}


In the virtual memory version, the virtual matrix
clsass\index{VMatrix } is derived from the virtual vector. Thus,
before the VMatrix is constructed, the constructor for a vdoub
is called. The members of vdoub are private to the VMatrix
class, so they are available to the member functions of VMatrix.
The matrix class stores its elements in a virtual vector in row
major order. That is, the rows of the matrix are stored
sequentially in the vector. Access and storage of the matrix
elements is achieved through vdoub::operator[].

A full set of functions is created for the matrix class.
Assignment of matrices and elements is developed. Functions
are also created for
\begin{enumerate}
 \item manipulating the matrix name, 
 \item for displaying matrix information, 
 \item for arithmetic operations,
 \item for data io, 
 \item for matrix functions, 
 \item for recursive calculations, 
 \item and nested function calls.  
\end{enumerate}
Much of this is accomplished using a stack of virtual matrices
which will be discussed in more detail later. For now, it is
sufficient to say that the matrix stack is a derived class of
VMatrix, but many of the public functions of the stack class are
used in VMatrix. Thus the two classes form something of a
symbiosis rather than a strict hierarchy.  The stack base is
called Dispatch\index{MStack!Dispatch}, and all other stack
matrices are unaccessible to other functions.

\section{The Declaration of the VMatrix Class}

The virtual matrix VMatrix class is created by the following declaration
\begin{verbatim}
class VMatrix : private vdoub
 {

 protected:
   unsigned signature;
   long mindex( int i, int j){
        return   (long) (((long)(i-1))*((long)c)+((long)(j-1))) ; }
   long curvecind;       // current index for vector
   strtype name;
   void SetupVectors( int rr=1, int cc=1 );
   void PurgeVectors( void );
   void Replace( VMatrix &ROp );
   void NewReference(VMatrix &ROp );

 public:

   double Nrerror( const char *errormsg );
   void Garbage( const char *errormsg);
   void Garbage( void ){ Garbage(" ");}
   int r,c;

   double m( int i, int j ) {          // get a value
       if( i<1 || j<1 || i>r || j>c ){
       cerr << "index out of range\n";
       exit(1);
       }
       curvecind = mindex(i,j);
       return v( curvecind );
   }
   VMatrix& M( int i, int j ) { // load current buffer
        curvecind = mindex(i,j);
        vdoub::operator[](curvecind);
        return *this;
   }

   VMatrix( void );            // constructors and destructors
   VMatrix( const char *str, int i, int j);
   VMatrix( VMatrix &ROp );
   ~VMatrix( void );

   double operator = ( double d ){     // write in a value
     vdoub::operator[](curvecind) = d;
     return d;
   }
   void operator= (VMatrix &ROp);

   void Nameit( const char *str ) { name = str; }
   void Nameit( VMatrix &mat ){ name = mat.name; }
   void Nameit( strtype newname){ name = newname; }
   strtype Getname( VMatrix &mat){ return mat.name; }
   void Showname( void ){ name.Showname(); }
   char *StringAddress( void ) { return name.StringAddress(); }

   void LoadMat(void);
   void DisplayMat(void);
   void InfoMat( void );
   void Writeb (char *fid, VMatrix &mat);

   friend VMatrix& operator+ (VMatrix &LOp, VMatrix &ROp);
   friend VMatrix& operator+ (double scalar, VMatrix &ROp);
   friend VMatrix& operator+ (VMatrix &ROp, double scalar);
   friend VMatrix& operator- (VMatrix &LOp, VMatrix &ROp);
   friend VMatrix& operator- (double scalar, VMatrix &ROp);
   friend VMatrix& operator- (VMatrix &ROp, double scalar);
   friend VMatrix& operator- (VMatrix &ROp);
   friend VMatrix& operator* (VMatrix &LOp, VMatrix &ROp);
   friend VMatrix& operator* (double scalar, VMatrix &ROp);
   friend VMatrix& operator* (VMatrix &ROp, double scalar);
   friend VMatrix& operator% (VMatrix &LOp, VMatrix &ROp );
   friend VMatrix& operator/ (VMatrix &LOp, VMatrix &ROp);
   friend VMatrix& operator/ (VMatrix &ROp, double scalar);
   friend VMatrix& operator/ (double scalar, VMatrix &ROp);
};
\end{verbatim}

The in-ram version has a very similar structure, but is not
derived from any other class. The system supplies dynamic memory
so there is no need to encapsulate the memory manager in a
separate class. In the in-ram version, matrices are constructed
from vectors allocated on the heap. The in--ram \verb+vdoub+
counts the number of references to the vector. The vector is not
freed until the number of references is
zero.\index{VMatrix!IN-RAM vdoub}


\section{Protected Parts of VMatrix}
The protected part of VMatrix must be accessible to the
matrix stack. The signature functions as in the vmm and
vdoub. It provides an indication of how corrupt or active
the matrix has become.

The function mindex( int i, int j)\index{VMatrix!mindex()}
calculates a vdoub vector index.\footnote{ It appears a bit
heavy handed with all of the casts to long ints, but there
were some problems of improper type promotion encountered
during development.} Stripped of the casts, the function is
simply \verb# (i-1)*c + (j-1) #, where i and j are matrix
indices, and c is the number of columns of the matrix.
Matrices are indexed from 1 to r, and 1 to c, so there are
r*c elements in the corresponding vdoub. However the vdoub
indices range from 0 to r*c - 1, so the i-1 and j-1 force
the minimum matrix indices to be 1. This calculation forces
the matrix to be in row major order. The
curvecind\index{VMatrix!curvecind} is a long integer
containing the vdoub index of the current matrix element.

The name\index{VMatrix!name} is a string type which can be
manipulated in the matrix functions. The functions are
found in the public part, but are Nameit(), Getname(),
Showname(), and StringAddress().

The next three functions control the allocation of the
vdoub vectors. Replace()\index{VMatrix!Replace()} calls
PurgeVectors()\index{VMatrix!PurgeVectors()} and
SetupVectors()\index{VMatrix!SetupVectors()} if the VMatrix
being replaced has different dimensions than the replacing
VMatrix. Replace() is called by the
operator=()\index{VMatrix!assignment}.  PurgeVectors()
calls vfree() to free the hdr in the protected part of
vdoub. SetupVectors() calls vmalloc() to create a new hdr
for vdoub.\footnote{ This ability to operate on the hdr of
vdoub clarified the meaning of inheritance for me since a
VMatrix is really just a virtual vector with some extra
functionality. It also shows why the hdr must be protected
rather than private in vdoub.}

For the in-ram version, SetupVectors() allocates r*c+1 elements
on the heap using calloc(). The first element is set to
SIGNATURE, then the base is shifted by 1. Thus, m[-1] is a check
value. If the value is not equal to SIGNATURE, then the matrix
has been deleted or clobbered.

The function \verb+NewReference()+
\index{VMatrix!NewReference()} is used by the matrix stack to
set up a reference to a base vector instead of a copy. It is
part of the deep copy technique used by push.

\section{The Public Part of VMatrix}

\subsection{Error Detection}

Don't let OOP fool you. It's pointer {\em intensive}. The
approach taken here is to stop the program at the slightest
problem. This conservative approach results from several
years experience with rogue strings, dangling pointers to
freed space, and jumping to execution of the heap because a
stray bit of data landed in my interrupt vector table. These
are some of the less benign things that can happen when
using pointers.

Two error detection routines are provided for a graceful
exit from the program. These routines print a message of
which routine was active when the program failed.
Nrerror()\index{VMatrix!Nrerror} prints an error message
and exits to the system with error code 1. Note that it is
declared to return a double. This is so Nrerror() can be
placed in the ternary operator '?:' where a double is to be
returned. For instance, suppose x= -1.0, then the statement 
\begin{verbatim}
   double y = ( x>0 ) ? log(x) : Nrerror("LOG: negative value");
\end{verbatim}
will compile correctly since the compiler thinks Nrerror()
will return a double. The compiler does not know that the
program will exit before the return. Nrerror() provides a
similar function in Press\cite{Pr:nrc} which is an
excellent source for code to translate to VMatrix
functions.\footnote{If you force a zero division in this
routine, VAX CC will cause a print of the function call
sequence to the error output stream}

There are two functions called
Garbage()\index{VMatrix!Garbage()}. They test the integrity
of the matrix arguments in function calls. It prints the
routine where a garbage matrix was found, then calls
Nrerror(). A good example of this function is found in the
assignment function for VMatrix:
\begin{verbatim}
void VMatrix::operator= (VMatrix &ROp)
{
   Garbage("operator =");
   ROp.Garbage("operator =");
   .
   .
   .
}
\end{verbatim}
First the member function of the VMatrix *this is called.
Next the ROp member function is called. The signature is
checked for validity. If it is wrong, the program stops.
You should make a habit of checking for incoming matrices
that could cause the program to bomb by using the Garbage()
functions.

\subsection{Indexing and Assigning Elements}

The indexing range is pretty simple. Matrix elements are
indexed from 1 to c, and from 1 to r. There are some
tradeoffs here. I could not figure out how to overload a
double bracket operator, [][], so I resorted to a
functional approach to indexing. Two index functions were
created, one for each side of an equals operator. The
little m(i,j) returns a double. The big M(i,j) returns a
reference to a VMatrix. Both functions call mindex() and
set curvecind.  The difference between the two operators is
that m(i,j) calls the vdoub member function v() for
returning a matrix element. The function M(i,j) calls the
vdoub bracket operators. You can think of
M(i,j)\index{VMatrix!M()} and m(i,j)\index{VMatrix!m()} as
surrogates for vdoub::operator[] and the vdoub function
v().

To assign an element to a matrix, use statements like
\begin{verbatim}
                MAT1.M(i,j) = MAT2.m(i,j); 
                MAT1.M(i,j) = adouble;
                Matpointer1->M(i,j) = Matpointer2->m(i,j);
                Matpointer1->M(i,j) = another_double;
\end{verbatim}
where MAT1 and MAT2 are virtual matrices, and Matpointers
are pointers to virtual matrices. The basic idea at work
here is that a virtual matrix is derived from the virtual
vector. The vector functions are available to the matrix
class, and so assignment of matrix elements works just as
in virtual vectors. The only difference is that the index
has to be calculated before the vector functions can be
called. The code for assigning an element is
\begin{verbatim}
   double operator = ( double d ){     // write in a value
     vdoub::operator[](curvecind) = d;
     return d;
   }
\end{verbatim}
The scope resolution operator is required to tell the
compiler to use the equals operator from the parent class
rather than the VMatrix class. 

\subsection{Assignment of Matrices}

Matrix assignment uses statements like "A=B;" where A and B
are matrices. The operator is the function given by the code
\begin{verbatim}
void VMatrix::operator= (VMatrix &ROp)
{
    Garbage("operator =");
    ROp.Garbage("operator =");
    Replace( ROp );
    Dispatch->Cleanstack();
}
\end{verbatim}
The operands are checked for validity. The left operand has
its header replaced. Then, any intermediate stack matrices
are removed from the stack.

Cleanstack() will be discussed in the chapter on the matrix
stack. Basically it serves to dispose of the intermediate
results of a repeated calculation. Suppose you use the
statement \verb# c=a+b+d #. The sum b+d is calculated, then
added to a, and the result is assigned to c. The
intermediate calculation of b+d is still on the heap in the
stack.  Cleanstack()\index{MStack!Cleanstack()} removes b+d and
a+(b+d) from the stack. 


\subsection{Constructors and Destructors}

There are three constructors for VMatrices. The default
constructor takes a void argument. Note that since a
VMatrix is a sub-class of a vdoub, the vdoub constructor is
called first. The extra VMatrix elements are supplied next.
The default constructor\index{VMatrix!constructor} sets the
name to 't', and allocates a 1x1 matrix.

The second form of the constructor is given in the code below:
\begin{verbatim}
VMatrix::VMatrix( const char *str, int rr, int cc):
    vdoub( ((long )((long) rr)*((long) cc)) )
  {
    r=rr;
    c=cc;
    curvecind = 0;
    Nameit( str );
    signature = SIGNATURE;
  }
\end{verbatim}
The second line of this constructor calls the vdoub
constructor with a default vector length of rr*cc. The
casts avoid type promotion problems. The constructor sets
the name to the supplied string. 

The third constructor is a copy constructor. This allows 
statements like \verb#VMatrix c = b;# to be used.
This creates the new matrix c, and then copies b into it.
The assignment may also be a function call, or repeated
operations. Matrix pointers can also be created by
functions such as \verb# VMatrix *Mat = new VMatrix;#.
Note that the copy constructor calls cleanstack(), so
the stack is cleared after the constructor is called. 
Thus, you must use \verb#Dispatch->Inclevel()# and
\verb#Dispatch->Declevel()# pairs in functions that call the
copy constructor.

The destructor\index{VMatrix!destructor} first calls the
destructor for the vdoub, returns, and sets the signature
to zero.

\subsection{The Naming Functions}

The naming functions store the name, retrieve the name,
display the name, or return the address of the name. The
function Nameit()\index{VMatrix!Nameit()} accepts string
constants, VMatrix references, or string objects. The
proper functions are called to assign the name to the
matrix. 

The function Getname()\index{VMatrix!getname()} returns a
string type object containing the name of the matrix. This
can be used in concatenation and assignment with other
strings. 

The function StringAddress() returns the address of the
matrix name. This function is used in the Writeb and Writea
functions where an address is used. This function may not
be of much use otherwise. You should not manipulate the
name directly, even though you know its base address. Let
the string class member functions do it for you.

\subsection{Other Service Functions}

Some other service functions are available. The first is
Loadmat()\index{VMatrix!Loadmat()} which allows the user to
enter a matrix interactively. To use this function to enter
the elements of a matrix, use code such as \verb#
A.Loadmat()# or \verb# A->Loadmat()#. The program will stop
and let you enter the matrix elements at the terminal. 

DisplayMat()\index{VMatrix!DisplayMat()} and
InfoMat()\index{VMatrix!Infomat()} display the matrix and
the dimensions of the active matrices. You can add more
checks as you want. Use the functions
Setwid()\index{functions!Setwid()} and Setdec()\index{functions!Setdec()} to
control the displayed width of the numbers and the number
of decimals.

Writeb(char *fid, VMatrix \&mat)\index{functions!Writeb()} doesn't
really belong as member function, but needs access to some
of the private members of the matrix class. It writes a
matrix to disk in a binary form. It writes the data out in
blocks, rather than elementwise. The fid is a file name to
write to.

\section{Arithmetic Friend Functions}

Addition, subtraction, unary minus, elementwise
multiplication, elementwise division, and matrix
multiplication is supported through friend
functions\index{VMatrix!$+,-,*,\%,/$}. The functions allow
scalar arguments as left or right operands. The
only unusual feature of these operators is the elementwise
multiplication is called using the '\%' (ampersand)
operator. They all return VMatrix references. See
Shildt\cite{sh:c++} for an explanation of how and why
friend functions are needed to overload these functions.

In general, these functions follow the same precedence as
the standard functions on doubles. They also produce
intermediate values on the stack. The code for addition is
given by
\begin{verbatim}
VMatrix& operator+ (VMatrix &LOp, VMatrix &ROp)
  {
    LOp.Garbage("operator +");
    ROp.Garbage("operator +");
    if (LOp.r==ROp.r && LOp.c==ROp.c) {
      VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
      for (int i=1; i<=LOp.r; ++i) {
          for (int j=1; j<=ROp.c; ++j) {
              temp->M(i,j) = LOp.m(i,j) + ROp.m(i,j);
          }
      }
      temp->name = temp->name+LOp.name+"+"+ROp.name+")";
      Dispatch->Push(*temp);
      delete temp;
    }
  else ROp.Nrerror("Mismatched Matrix sizes in addition function\n");
  return Dispatch->ReturnMat();
  }
\end{verbatim}
The function first checks the validity of its operands. It
also checks that the matrices are conformable. If they are,
a new matrix is put on the heap. The matrices are summed.
Note the use of the little m(i,j) on the right side of the
sum, and the M(i,j) on the left. The name is created. The
matrix is pushed onto the stack, then deleted. 

The reference returned is the matrix pointed to by
\verb# Dispatch->next#. The function ReturnMat() passes the matrix
reference, and does not change the function nesting level.
Basically, the nesting level keeps track of the function depth
where addition was called. Also, Cleanstack() uses the nesting
level to determine which matrices to strip from the stack.
System functions should use ReturnMat() rather than DecReturn().
DecReturn() decrements the nesting level, then returns the
correct matrix reference. More about this in the next chapter.
The way to tell when to use DecReturn() is if you enter the
function and call Inclevel(). Call DecReturn() if you called
Inclevel(). Using the stack allows the functions to be
recursive.