Implement fill-stack

knodlseder
1 parent 4f37bc98
Showing 16 changed files with 897 additions and 390 deletions
doc/gammalib_um.tex
src/GMatrix.cpp
src/GMatrix.hpp
src/GSparseMatrix.cpp
src/GSparseMatrix.hpp
src/GSparseSymbolic.cpp
src/GSymMatrix.cpp
src/GSymMatrix.hpp
src/GammaLib.hpp
src/Makefile.am
test/Makefile.am
test/Makefile.in
test/test_GMatrix.cpp
test/test_GSparseMatrix.cpp
test/test_GSymMatrix.cpp
test/test_GVector.cpp
-%---------------------------------------------------------------------------- 
%                            GammaLib User Manual
%
% Version .........: 0.0.0
% Last modification: 26 October 2006
%----------------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definitions for manual package
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\task}{\mbox{GammaLib}}
\newcommand{\this}{\mbox{\tt \task}}
\newcommand{\version}{\mbox{0.0.0}}
\newcommand{\calendar}{\mbox{26 October 2006}}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Document definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass{article}[12pt,a4]
\usepackage{epsfig}
\usepackage{manual}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Begin of document body
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\frontpage


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Introduction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Motivation}

Again another library?
Indeed, 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Philisophy}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Implementation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GammaLib objects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{\this\ objects}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Vectors and matrices}

%----------------------------------------------------------------------------
\subsubsection{General}

Most of the calculations that are needed in data analysis can be 
expressed as vector and matrix operations.
For this reason, \this\ provides a vector and matrix package that is 
implemeted as the {\tt GVector} and the {\tt GMatrix} classes.


%----------------------------------------------------------------------------
\subsubsection{Matrix types}

The matrix base class is {\tt GMatrix} which stores all elements of 
the matrix (including zeros) in a full array.
The full array has the dimensions {\tt rows} $\times$ {\tt cols}, 
matrix elements are stored in column major format.

{\tt GSymMatrix} symmetric matrix.

{\tt GSparseMatrix} sparse matrix.

{\tt GSymSparseMatrix} symmetric sparse matrix.

%----------------------------------------------------------------------------
\subsubsection{Matrix arithmetics}

{\tt GMatrix} and its derived objects can be added, subtracted, 
scaled and negated, very much as other data types.
The following code fragment illustrates the usage:
\begin{verbatim}
 // Matrix allocation
 GMatrix a(2,3);
 GMatrix b(3,2);
 GMatrix c(2,2);
 GMatrix d(2,2);
 
 // Element assignment
 a(0,0) = a(1,1) = a(1,2) = 1.0;
 b(0,0) = b(1,1) = b(2,1) = 1.0;
 d(0,0) = d(1,1) = 1.0;
 
 // Arithmetics
 c = a * b;           // Matrix multiplication
 d = c;               // Assignment
 c = c * 2.0;         // Multiplication of all matrix elements by 2.0
 a = a + 1.0;         // Addition of 1.0 to all matrix elements
 d = d + d - c;       // Matrix addition and subtraction
 c = -d;              // Negation of all matrix elements
 a *= 5.0;            // Multiplication of all matrix elements by 5.0
\end{verbatim}
Matrix allocation requires {\bf always} the specification of the number 
of rows and columns of the matrix.
This means that there is no such thing as a matrix object with zero 
rows or columns.
All matrix elements are initialised to 0.0 by the matrix allocation.

Matrix elements are accessed by the {\tt matrix(row,col)} function,
where {\tt row} and {\tt col} start from 0 for the first row or column 
and run up to the number of rows or columns minus 1.

Matrix arithmetics comprise binary operators {\tt C = A operator B} and
unary operators {\tt C operator A}. {\tt A} and {\tt B} may be either
{\tt GMatrix} (or derived) objects or {\tt double} values. In case of
{\tt double} values the operation is performed on each matrix element,
i.e.
\begin{verbatim}
 c = c * 2.0;
 d = d / 3.0;
\end{verbatim}
multiply each matrix element of {\tt c} by 2.0 and divide each matrix 
element of {\tt d} by 3.0.
These operations can be writted also more compactly using unary 
operators:
\begin{verbatim}
 c *= 2.0;
 d /= 3.0;
\end{verbatim}
While the addition, multiplication

%----------------------------------------------------------------------------
\subsubsection{Matrix functions}


%----------------------------------------------------------------------------
\subsubsection{Matrix compression}

The {\tt GSymMatrix} Cholesky decomposition, solver and inversion
routines may also be applied to matrices that contain rows or
columns that are filled by zeros.
In this case the functions provide the option to (logically)
compress the matrices by skipping the zero rows and columns during
the calculation.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code reference
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Code reference}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{{\tt GVector}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{{\tt GMatrix}}

\end{document} 
 \ No newline at end of file
+%---------------------------------------------------------------------------- 
%                             GammaLib User Manual
%
% Version .........: 0.0.0
% Last modification: 26 October 2006
%----------------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definitions for manual package
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\task}{\mbox{GammaLib}}
\newcommand{\this}{\mbox{\tt \task}}
\newcommand{\version}{\mbox{0.0.0}}
\newcommand{\calendar}{\mbox{26 October 2006}}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Document definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass{article}[12pt,a4]
\usepackage{epsfig}
\usepackage{manual}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Begin of document body
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\frontpage


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Introduction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Motivation}

Again another library?
Indeed, 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Philisophy}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Implementation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GammaLib objects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{\this\ objects}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Vectors and matrices}

%----------------------------------------------------------------------------
\subsubsection{General}

Most of the calculations that are needed in data analysis can be 
expressed as vector and matrix operations.
For this reason, \this\ provides a vector and matrix package that is 
implemeted as the {\tt GVector} and the {\tt GMatrix} classes.


%----------------------------------------------------------------------------
\subsubsection{Matrix storage classes}

The matrix base class is {\tt GMatrix} which stores all elements of 
the matrix (including zeros) in a full array.
The full array has the dimensions {\tt rows} $\times$ {\tt cols}, 
matrix elements are stored in column major format.

Derived matrix classes exist for more compact storage modes, reducing 
the memory requirements for handling and storing the information.

{\tt GSymMatrix} symmetric matrix.

{\tt GSparseMatrix} sparse matrix.

The following description of the matrix class properties apply to all 
matrix storage classes, if not stated otherwise. 


%----------------------------------------------------------------------------
\subsubsection{Matrix arithmetics}

{\tt GMatrix} and its derived objects can be added, subtracted, 
scaled and negated, very much as other data types.
The following code fragment illustrates the usage:
\begin{verbatim}
 // Matrix allocation
 GMatrix A(2,3);
 GMatrix B(3,2);
 GMatrix C(2,2);
 GMatrix D(2,2);
 
 // Element assignment
 A(0,0) = A(1,1) = A(1,2) = 1.0;
 B(0,0) = B(1,1) = B(2,1) = 1.0;
 D(0,0) = D(1,1) = 1.0;
 
 // Arithmetics
 C = A * B;                 // Matrix multiplication
 D = C;                     // Assignment
 C = C * 2.0;               // Multiplication of all matrix elements by 2.0
 D = D + D - C;             // Matrix addition and subtraction
 C = -D;                    // Negation of all matrix elements
 A *= 5.0;                  // Multiplication of all matrix elements by 5.0
\end{verbatim}
Matrix allocation requires {\bf always} the specification of the number 
of rows and columns of the matrix.
This means that there is no such thing as a matrix object with zero 
rows or columns.
All matrix elements are initialised to 0.0 by the matrix allocation.

Matrix elements are accessed by the {\tt matrix(row,col)} function,
where {\tt row} and {\tt col} start from 0 for the first row or column 
and run up to the number of rows or columns minus 1.

Matrix arithmetics comprise binary operators {\tt C = A operator B} and
unary operators {\tt C operator A}. {\tt A} and {\tt B} may be either
{\tt GMatrix} (or derived) objects or {\tt double} values. The 
following fundamental operations exist for all matrices:
\begin{verbatim}
 C = A + B;                 // Matrix addition
 C = A - B;                 // Matrix subtraction
 C = A * B;                 // Matrix multiplication
 C = A * v;                 // Matrix vector multiplication
 C = A * s;                 // Matrix scalar multiplication
 C = s * A;                 // Matrix scalar multiplication (left sided)
 C = A * s;                 // Matrix scalar multiplication (right sided)
 C = A / s;                 // Matrix scalar division
 C = -A;                    // Negation
 A += B;                    // Matrix inplace addition
 A -= B;                    // Matrix inplace subtraction
 A *= B;                    // Matrix inplace multiplications
 A *= s;                    // Matrix inplace scalar multiplication
 A /= s;                    // Matrix inplace scalar division
\end{verbatim}
Note that all {\bf inplace operations}, such as the 
{\tt +=},
{\tt -=},
{\tt *=}, and
{\tt /=} operators do not require the allocation of additional
memory to hold the matrix result elements, hence inplace operations
may be preferable if memory requirements are an issue.

Matrices may be compared to each other using
\begin{verbatim}
 int identical = (A == B);  // True if all matrix elements are identical
                            // and if the matrices have identical dimensions
 int different = (A != B);  // True if at least one matrix element differs
                            // or if the matrices have different dimensions
\end{verbatim}

%----------------------------------------------------------------------------
\subsubsection{Matrix functions}

Matrix attributes:
\begin{verbatim}
 unsigned rows = M.rows();   // Returns number of rows in matrix
 unsigned cols = M.cols();   // Returns number of columns in matrix
\end{verbatim}

Matrix information:
\begin{verbatim}
 double sum = M.sum();       // Sum of all elements in matrix
 double min = M.min();       // Returns minimum element of matrix
 double max = M.max();       // Returns maximum element of matrix
\end{verbatim}

%----------------------------------------------------------------------------
\subsubsection{Matrix factorisations}

A general tool of numeric matrix calculs is factorisation.

Solve linear equation {\tt Ax = b}.
Inverse a matrix (by solving successively {\tt Ax = e}, where {\tt e} 
are the unit vectors for all dimensions).

%
\paragraph{Matrix factorisations}

For symmetric and positive definite matrices the most efficient 
factorisation is the Cholesky decomposition.
The following code fragment illustrates the usage:
\begin{verbatim}
 GMatrix A(n_rows, n_cols};
 GVector x(n_rows};
 GVector b(n_rows};
 ...
 A.cholesky_decompose();     // Perform Cholesky factorisation
 x = A.cholesky_solver(b);   // Solve Ax=b for x
\end{verbatim}
Note that once the function {\tt A.cholesky\_decompose()} has been 
applied, the original matrix content has been replaced by its Cholesky 
decomposition.
Since the Cholesky decomposition can be performed inplace (i.e. 
without the allocation of additional memory to hold the result), the 
matrix replacement is most memory economic.
In case that the original matrix should be kept, one may either copy 
it before into another {\tt GMatrix} object or use the function
\begin{verbatim}
 GMatrix L = cholesky_decompose(A);
 x = L.cholesky_solver(b);
\end{verbatim}

A symmetric and positif definite matrix can be inverted using the 
Cholesky decomposition using
\begin{verbatim}
 A.cholesky_invert();        // Inverse matrix using Cholesky factorisation
\end{verbatim}
Alternatively, the function 
\begin{verbatim}
 GMatrix A_inv = cholesky_invert(A);
\end{verbatim}
may be used.

The Cholesky decomposition, solver and inversion
routines may also be applied to matrices that contain rows or
columns that are filled by zeros.
In this case the functions provide the option to (logically)
compress the matrices by skipping the zero rows and columns during
the calculation.

For compressed matrix Cholesky factorisation, only the non-zero rows 
and columns have to be symmetric and positive definite.
In particular, the full matrix may even be non-symmetric.



%----------------------------------------------------------------------------
\subsubsection{Sparse matrix}

%%%%%
\paragraph{General}

\this\ provides a sparse matrix storage type that can save 
considerable amounts of memory for calculations.

Assigment possible.

The only exception that does not work is
\begin{verbatim}
 GSparseMatrix A;
 A(0,0) = A(1,1) = A(2,2) = 1.0;
\end{verbatim}
In this case the value {\tt 1.0} is only assigned to the last 
element, i.e. {\tt A(2,2)}.

Inverting a sparse matrix produces in general a full matrix, so the 
inversion function should be used with caution.
Note that a full matrix that is stored in sparse format takes roughly
twice the memory than a normal {\tt GMatrix} object.
If nevertheless the inverse of a sparse matrix should be examined, it 
is recommended to perform the analysis column-wise:
\begin{verbatim}
 GSparseMatrix A(rows,cols);                // Allocate sparse matrix
 GVector       unit(rows);                  // Allocate vector
 ...
 A.cholesky_decompose();                    // Factorise matrix

 // Column-wise solving the matrix equation
 for (int col = 0; col < m_cols; ++col) {
   unit(col) = 1.0;                         // Set unit vector
   GVector x = cholesky_solver(unit);       // Get column x of inverse
   ...
   unit(col) = 0.0;                         // Clear unit vector for next round
 }
\end{verbatim}


%%%%%
\paragraph{Filling of a sparse matrix}

The filling of a sparse matrix is a tricky issue since the storage
of the elements depends on their distribution in the matrix.
If one would know beforehand this distribution, sparse matrix filling
would be easy and fast.
In general, however, the distribution is not known a priori, and
matrix filling may become a quite time consuming task.

If a matrix has to be filled element by element, the access through
the operator
\begin{verbatim}
 m(row,col) = value;
\end{verbatim}
may be mandatory.
In principle, if a new element is inserted into a matrix a new memory
cell has to be allocated for this element, and other elements may be 
moved.
Memory allocation is quite time consuming, and to reduce the overhead,
{\tt GSparseMatrix} can be configured to allocate memory in bunches.
By default, each time more matrix memory is needed, {\tt GSparseMatrix}
allocates 512 cells at once (or 6144 Bytes since each element requires 
a {\tt double} and a {\tt int} storage location).
If this amount of memory is not adequat one may change this value by
using
\begin{verbatim}
 m.set_mem_block(size);
\end{verbatim}
where {\tt size} is the number of matrix elements that should be
allocated at once (corresponding to a total memory of 
$12 \times {\tt size}$ Bytes).

Alternatively, a matrix may be filled column-wise using the functions
\begin{verbatim}
 m.insert_col(vector,col);                  // Insert a vector in column
 m.add_col(vector,col);                     // Add content of a vector to column
\end{verbatim}
While {\tt insert\_col} sets the values of column {\tt col} (deleting 
thus any previously existing entries), {\tt add\_col} adds the content 
of {\tt vector} to all elements of column {\tt col}.
Using these functions is considerably more rapid than filling individual 
values.

Still, if the matrix is big (i.e. severeal thousands of rows and 
columns), filling individual columns may still be slow.
To speed-up dynamical matrix filling, an internal fill-stack has been 
implemented in {\tt GSparseMatrix}.
Instead of inserting values column-by-column, the columns are stored 
in a stack and filled into the matrix once the stack is full.
This reduces the number of dynamic memory allocations to let the 
matrix grow as it is built.
By default, the internal stack is disabled.
The stack can be enabled and used as follows:
\begin{verbatim}
 m.stack_init(size, entries);               // Initialise stack
 ...
 m.add_col(vector,col);                     // Add columns
 ...
 m.stack_destroy();                         // Flush and destory stack
\end{verbatim}
The method {\tt stack\_init} initialises a stack with a number of 
{\tt size} elements and a maximum of {\tt entries} columns.
The larger the values {\tt size} and {\tt entries} are chosen, the
more efficient the stack works.
The total amount of memory of the stack can be estimated as
$12 \times {\tt size} + 8 \times {\tt entries}$ Bytes.
If a rough estimate of the total number of non-zero elements is
available it is recommended to set {\tt size} to this value.
As a rule of thumb, {\tt size} should be at least of the dimension of 
either the number of rows or the number of columns of the matrix 
(take the maximum of both).
{\tt entries} is best set to the number of columns of the matrix.
If memory limits are an issue smaller values may be set, but if the
values are too small, the speed increase may become negligible (or 
even, stack-filling may become slower than normal filling).

Stack-filling only works with the method {\tt add\_col}.
Note also that filling sub-sequently the same row leads to stack 
flushing, so if the same column is filled in a loop, i.e.
\begin{verbatim}
 for (int col = 0; col < 100; ++col) {
   column    = 0.0;                         // Reset column
   column(i) = i;                           // Set column
   m.add_col(column,col);                   // Add column
 }   
\end{verbatim}
it is better to fill the column before adding, i.e.
\begin{verbatim}
 column    = 0.0;                           // Reset column
 for (int col = 0; col < 100; ++col)
   column(i) = i;                           // Set column
\end{verbatim}
The method {\tt add\_col} calls the method {\tt stack\_push\_column}
for stack filling.
{\tt add\_col} is more general than {\tt stack\_push\_column} in that 
it decides which of stack- or direct filling is more adequate.
In particular, {\tt stack\_push\_column} may refuse pushing a column 
onto the stack if there is not enough space.
In that case, {\tt stack\_push\_column} returns a non-zero value that 
corresponds to the number of non-zero elements in the vector that 
should be added.
However, it is recommended to not use {\tt stack\_push\_column} and 
call instead {\tt add\_col}.

The method {\tt stack\_destroy} is used to flush and destroy the
stack. 
After this call the stack memory is liberated.
If the stack should be flushed without destroying it, the method
{\tt stack\_flush} may be used:
\begin{verbatim}
 m.stack_init(size, entries);               // Initialise stack
 ...
 m.add_col(vector,col);                     // Add columns
 ...
 m.stack_flush();                           // Simply flush stack
\end{verbatim}
Once flushed, the stack can be filled anew.

Note that stack flushing is not automatic!
This means, if one trys to use a matrix for calculs without flushing, 
the calculs may be work.
{\bf If a stack is used for filling, always flush the stack before using 
the matrix.}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code reference
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Code reference}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{{\tt GVector}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{{\tt GMatrix}}

\end{document} 
 \ No newline at end of file
@@ -20,6 +20,96 @@ using namespace std;
 /*==========================================================================
  =                                                                         =
+ =                            Private functions                            =
+ =                                                                         =
+ ==========================================================================*/
+void GMatrix::init_members(void)
+{
+  // Initialise GMatrix members
+  m_rows       = 0;      // Number of rows
+  m_cols       = 0;      // Number of columns
+  m_elements   = 0;      // Logically used number of elements
+  m_alloc      = 0;      // Allocated # of elements (>= m_elements)
+  m_num_rowsel = 0;      // Number of selected rows (for comp. decomp.)
+  m_num_colsel = 0;      // Number of selected columns (for comp. decomp.)
+  m_data       = NULL;   // Matrix data
+  m_colstart   = NULL;   // Column start indices (m_cols+1)
+  m_rowsel     = NULL;   // Row selection (for compressed decomposition)
+  m_colsel     = NULL;   // Column selection (for compressed decomposition)
+  
+  // Return
+  return;
+}
+/***************************************************************************/
+void GMatrix::copy_members(const GMatrix& m)
+{
+  // Copy matrix attributes
+  m_rows       = m.m_rows;
+  m_cols       = m.m_cols;
+  m_elements   = m.m_elements;
+  m_alloc      = m.m_alloc;
+  m_num_rowsel = m.m_num_rowsel;
+  m_num_colsel = m.m_num_colsel;
+
+  // Allocate memory for column start array and copy content
+  m_colstart = new int[m_cols+1];
+  if (m_colstart == NULL)
+    throw mem_alloc("GMatrix::copy_members(const Matrix&)", m_cols+1);
+  for (int i = 0; i <= m_cols; ++i)
+    m_colstart[i] = m.m_colstart[i];
+  
+  // Allocate memory for elements and copy them (only if there are elements)
+  if (m.m_data != NULL && m_alloc > 0) {
+    m_data = new double[m_alloc];
+    if (m_data == NULL)
+	  throw mem_alloc("GMatrix::copy_members(const Matrix&)", m_alloc);
+    for (int i = 0; i < m_elements; ++i)
+      m_data[i] = m.m_data[i];
+  }
+  
+  // If there is a row selection then copy it
+  if (m.m_rowsel != NULL && m_num_rowsel > 0) {
+    m_rowsel = new int[m_num_rowsel];
+    if (m_rowsel == NULL)
+	  throw mem_alloc("GMatrix::copy_members(const Matrix&)", m_num_rowsel);
+    for (int i = 0; i < m_num_rowsel; ++i)
+      m_rowsel[i] = m.m_rowsel[i];
+  }
+
+  // If there is a column selection then copy it
+  if (m.m_colsel != NULL && m_num_colsel > 0) {
+    m_colsel = new int[m_num_colsel];
+    if (m_colsel == NULL)
+	  throw mem_alloc("GMatrix::copy_members(const Matrix&)", m_num_colsel);
+    for (int i = 0; i < m_num_colsel; ++i)
+      m_colsel[i] = m.m_colsel[i];
+  }
+
+  // Return
+  return;
+}
+/***************************************************************************/
+void GMatrix::free_members(void)
+{
+  // De-allocate only if memory has indeed been allocated
+  if (m_colsel   != NULL) delete [] m_colsel;
+  if (m_rowsel   != NULL) delete [] m_rowsel;
+  if (m_data     != NULL) delete [] m_data;
+  if (m_colstart != NULL) delete [] m_colstart;
+
+  // Properly mark members as free
+  m_colstart = NULL;
+  m_data     = NULL;
+  m_rowsel   = NULL;
+  m_colsel   = NULL;
+  
+  // Return
+  return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
  =                      GMatrix constructors/destructors                   =
  =                                                                         =
  ==========================================================================*/
@@ -30,17 +120,8 @@ using namespace std;
 GMatrix::GMatrix(void)
 {
   // Initialise private members for clean destruction
-  m_rows       = 0;
-  m_cols       = 0;
-  m_elements   = 0;
-  m_alloc      = 0;
-  m_num_rowsel = 0;
-  m_num_colsel = 0;
-  m_data       = NULL;
-  m_colstart   = NULL;
-  m_rowsel     = NULL;
-  m_colsel     = NULL;
-	
+  init_members();
+
   // Return
   return;
 }
@@ -52,16 +133,7 @@ GMatrix::GMatrix(void)
 GMatrix::GMatrix(int rows, int cols)
 {
   // Initialise private members for clean destruction
-  m_rows       = 0;
-  m_cols       = 0;
-  m_elements   = 0;
-  m_alloc      = 0;
-  m_num_rowsel = 0;
-  m_num_colsel = 0;
-  m_data       = NULL;
-  m_colstart   = NULL;
-  m_rowsel     = NULL;
-  m_colsel     = NULL;
+  init_members();
   // Determine number of elements to store in matrix
   int elements = rows*cols;
@@ -106,59 +178,11 @@ GMatrix::GMatrix(int rows, int cols)
 GMatrix::GMatrix(const GMatrix& m)
 {
   // Initialise private members for clean destruction
-  m_rows       = 0;
-  m_cols       = 0;
-  m_elements   = 0;
-  m_alloc      = 0;
-  m_num_rowsel = 0;
-  m_num_colsel = 0;
-  m_data       = NULL;
-  m_colstart   = NULL;
-  m_rowsel     = NULL;
-  m_colsel     = NULL;
-  
-  // Allocate memory for column start array and copy content
-  m_colstart = new int[m.m_cols+1];
-  if (m_colstart == NULL)
-    throw mem_alloc("GMatrix copy constructor", m.m_cols+1);
-  for (int i = 0; i <= m.m_cols; ++i)
-    m_colstart[i] = m.m_colstart[i];
-  
-  // Allocate memory for elements and copy them (only if there are elements)
-  if (m.m_alloc > 0) {
-    m_data = new double[m.m_alloc];
-    if (m_data == NULL)
-	  throw mem_alloc("GMatrix copy constructor", m.m_alloc);
-    for (int i = 0; i < m.m_elements; ++i)
-      m_data[i] = m.m_data[i];
-  }
-  
-  // If there is a row selection then copy it
-  if (m.m_rowsel != NULL && m.m_num_rowsel > 0) {
-    m_rowsel = new int[m.m_num_rowsel];
-    if (m_rowsel == NULL)
-	  throw mem_alloc("GMatrix copy constructor", m.m_num_rowsel);
-    for (int i = 0; i < m.m_num_rowsel; ++i)
-      m_rowsel[i] = m.m_rowsel[i];
-	m_num_rowsel = m.m_num_rowsel;
-  }
+  init_members();
-  // If there is a column selection then copy it
-  if (m.m_colsel != NULL && m.m_num_colsel > 0) {
-    m_colsel = new int[m.m_num_colsel];
-    if (m_colsel == NULL)
-	  throw mem_alloc("GMatrix copy constructor", m.m_num_colsel);
-    for (int i = 0; i < m.m_num_colsel; ++i)
-      m_colsel[i] = m.m_colsel[i];
-	m_num_colsel = m.m_num_colsel;
-  }
+  // Copy members
+  copy_members(m);
-  // Copy matrix attributes
-  m_rows     = m.m_rows;
-  m_cols     = m.m_cols;
-  m_elements = m.m_elements;
-  m_alloc    = m.m_alloc;
-  
   // Return
   return;
 }
@@ -169,11 +193,8 @@ GMatrix::GMatrix(const GMatrix&amp; m)
  ***************************************************************************/
 GMatrix::~GMatrix()
 {
-  // De-allocate only if memory has indeed been allocated
-  if (m_rowsel   != NULL) delete [] m_rowsel;
-  if (m_colsel   != NULL) delete [] m_colsel;
-  if (m_colstart != NULL) delete [] m_colstart;
-  if (m_data     != NULL) delete [] m_data;
+  // Free members
+  free_members();
   // Return
   return;
@@ -193,66 +214,15 @@ GMatrix&amp; GMatrix::operator= (const GMatrix&amp; m)
 {
   // Execute only if object is not identical
   if (this != &m) { 
-  
-    // De-allocate memory if it has indeed been allocated
-    if (m_colstart != NULL) delete [] m_colstart;
-    if (m_rowsel   != NULL) delete [] m_rowsel;
-    if (m_colsel   != NULL) delete [] m_colsel;
-    if (m_data     != NULL) delete [] m_data;
+    // Free members
+    free_members();
+  
     // Initialise private members for clean destruction
-    m_rows       = 0;
-    m_cols       = 0;
-    m_elements   = 0;
-	m_alloc      = 0;
-    m_num_rowsel = 0;
-    m_num_colsel = 0;
-    m_data       = NULL;
-    m_colstart   = NULL;
-    m_rowsel     = NULL;
-    m_colsel     = NULL;
-  
-    // Allocate memory for column start array and copy content
-    m_colstart = new int[m.m_cols+1];
-    if (m_colstart == NULL)
-      throw mem_alloc("GMatrix::operator= (const GMatrix&)", m.m_cols+1);
-    for (int i = 0; i <= m.m_cols; ++i)
-      m_colstart[i] = m.m_colstart[i];
-  
-    // Allocate memory for elements and copy them (only if there are elements)
-    if (m.m_alloc > 0) {
-      m_data = new double[m.m_alloc];
-      if (m_data == NULL)
-	    throw mem_alloc("GMatrix::operator= (const GMatrix&)", m.m_alloc);
-      for (int i = 0; i < m.m_elements; ++i)
-        m_data[i] = m.m_data[i];
-    }
-
-    // If there is a row selection then copy it
-    if (m.m_rowsel != NULL && m.m_num_rowsel > 0) {
-      m_rowsel = new int[m.m_num_rowsel];
-      if (m_rowsel == NULL)
-  	    throw mem_alloc("GMatrix::operator= (const GMatrix&)", m.m_num_rowsel);
-      for (int i = 0; i < m.m_num_rowsel; ++i)
-        m_rowsel[i] = m.m_rowsel[i];
-	  m_num_rowsel = m.m_num_rowsel;
-    }
-
-    // If there is a column selection then copy it
-    if (m.m_colsel != NULL && m.m_num_colsel > 0) {
-      m_colsel = new int[m.m_num_colsel];
-      if (m_colsel == NULL)
-	    throw mem_alloc("GMatrix::operator= (const GMatrix&)", m.m_num_colsel);
-      for (int i = 0; i < m.m_num_colsel; ++i)
-        m_colsel[i] = m.m_colsel[i];
-	  m_num_colsel = m.m_num_colsel;
-    }
+    init_members();
-    // Copy matrix attributes
-    m_rows     = m.m_rows;
-    m_cols     = m.m_cols;
-    m_elements = m.m_elements;
-	m_alloc    = m.m_alloc;
+    // Copy members
+    copy_members(m);
   } // endif: object was not identical
@@ -137,6 +137,9 @@ protected:
   GMatrix();
   // Private member functions
+  void init_members(void);
+  void copy_members(const GMatrix& m);
+  void free_members(void);
   void select_non_zero(void);
   // Private data area
@@ -17,7 +17,7 @@
  * m_rowinx   : holds the row indices for all elements                (Ai) *
  * Column 'col' covers therefore [m_colstart[j], ..., m_colstart[j+1]-1]   *
  * Note that 'm_colstart' has m_cols+1 elements.                           *
-/***************************************************************************
+ ***************************************************************************/
 /* __ Includes ___________________________________________________________ */
 #include "GSparseMatrix.hpp"
@@ -27,31 +27,26 @@
 /* __ Namespaces _________________________________________________________ */
 using namespace std;
-/* __ Definitions ________________________________________________________ */
+/* __ Debug definitions __________________________________________________ */
 //#define G_DEBUG_SPARSE_PENDING                    // Analyse pending values
 //#define G_DEBUG_SPARSE_INSERTION                 // Analyse value insertion
 //#define G_DEBUG_SPARSE_COMPRESSION   // Analyse zero row/column compression
 //#define G_DEBUG_SPARSE_MALLOC                  // Analyse memory management
+//#define G_DEBUG_SPARSE_STACK_PUSH           // Analyse stack column pushing
+//#define G_DEBUG_SPARSE_STACK_FLUSH         // Analyse stack column flushing
+
+/* __ Definitions ________________________________________________________ */
 #define G_USE_MEMCPY
+
 /*==========================================================================
  =                                                                         =
- =                   GSparseMatrix constructors/destructors                =
+ =                            Private functions                            =
  =                                                                         =
  ==========================================================================*/
-
-/***************************************************************************
- *                        GSparseMatrix constructor                        *
- ***************************************************************************/
-GSparseMatrix::GSparseMatrix(int rows, int cols, int elements) : GMatrix()
+void GSparseMatrix::init_members(void)
 {
-  // Initialise private members for clean destruction
-  m_rows      = 0;
-  m_cols      = 0;
-  m_elements  = 0;             // Logically used # of elements
-  m_alloc     = 0;             // Allocated # of elements (>= m_elements)
-  m_data      = NULL;
-  m_colstart  = NULL;
+  // Initialise sparse matrix members
   m_rowinx    = NULL;
   m_mem_block = G_SPARSE_MATRIX_DEFAULT_MEM_BLOCK;
   m_zero      = 0.0;
@@ -60,7 +55,106 @@ GSparseMatrix::GSparseMatrix(int rows, int cols, int elements) : GMatrix()
   m_fill_col  = 0;
   m_symbolic  = NULL;
   m_numeric   = NULL;
-	
+  
+  // Initialise stack members
+  m_stack_max_entries = 0;
+  m_stack_size        = 0;
+  m_stack_entries     = 0;
+  m_stack_colinx      = NULL;
+  m_stack_start       = NULL;
+  m_stack_data        = NULL;
+  m_stack_rowinx      = NULL;
+  
+  // Return
+  return;
+}
+/***************************************************************************/
+void GSparseMatrix::copy_members(const GSparseMatrix& m)
+{
+  // Copy data members
+  m_mem_block = m.m_mem_block;
+  m_zero      = m.m_zero;
+  m_fill_val  = m.m_fill_val;
+  m_fill_row  = m.m_fill_row;
+  m_fill_col  = m.m_fill_col;
+
+  // Copy row indices
+  if (m.m_rowinx != NULL && m.m_alloc > 0) {
+    m_rowinx = new int[m.m_alloc];
+    if (m_rowinx == NULL)
+	  throw mem_alloc("GSparseMatrix::copy_members(const GSparseMatrix&)", m.m_alloc);
+    for (int i = 0; i < m.m_elements; ++i)
+      m_rowinx[i] = m.m_rowinx[i];
+  }
+
+  // Copy symbolic decomposition
+  if (m.m_symbolic != NULL) {
+    GSparseSymbolic* symbolic = new GSparseSymbolic();
+    if (symbolic == NULL)
+	  throw mem_alloc("GSparseMatrix::copy_members(const GSparseMatrix&)", 1);
+    *symbolic  = *((GSparseSymbolic*)m.m_symbolic);
+    m_symbolic = (void*)symbolic;
+  }
+
+  // Copy numeric decomposition
+  if (m.m_numeric != NULL) {
+    GSparseNumeric* numeric = new GSparseNumeric();
+    if (numeric == NULL)
+	  throw mem_alloc("GSparseMatrix::copy_members(const GSparseMatrix&)", 1);
+    *numeric  = *((GSparseNumeric*)m.m_numeric);
+    m_numeric = (void*)numeric;
+  }
+  
+  // Return
+  return;
+}
+/***************************************************************************/
+void GSparseMatrix::free_members(void)
+{
+  // De-allocate only if memory has indeed been allocated
+  if (m_numeric  != NULL) delete (GSparseNumeric*)m_numeric;
+  if (m_symbolic != NULL) delete (GSparseSymbolic*)m_symbolic;
+  if (m_rowinx   != NULL) delete [] m_rowinx;
+
+  // Properly mark members as free
+  m_rowinx   = NULL;
+  m_symbolic = NULL;
+  m_numeric  = NULL;
+
+  // Free stack members
+  if (m_stack_colinx != NULL) delete [] m_stack_colinx;
+  if (m_stack_start  != NULL) delete [] m_stack_start;
+  if (m_stack_data   != NULL) delete [] m_stack_data;
+  if (m_stack_rowinx != NULL) delete [] m_stack_rowinx;
+
+  // Properly mark members as free
+  m_stack_colinx = NULL;
+  m_stack_start  = NULL;
+  m_stack_data   = NULL;
+  m_stack_rowinx = NULL;
+
+  // Return
+  return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                   GSparseMatrix constructors/destructors                =
+ =                                                                         =
+ ==========================================================================*/
+
+/***************************************************************************
+ *                        GSparseMatrix constructor                        *
+ * ----------------------------------------------------------------------- *
+ * First invoques GMatrix void constructor which initialises all GMatrix   *
+ * members but which does not allocate memory or assign any value.         *
+ ***************************************************************************/
+GSparseMatrix::GSparseMatrix(int rows, int cols, int elements) : GMatrix()
+{
+  // Initialise private members for clean destruction
+  init_members();
+
   // Allocate column start array. This is the only array that we can
   // allocate at this time. The other arrays can only be allocated during
   // filling of the matrix
@@ -96,44 +190,11 @@ GSparseMatrix::GSparseMatrix(int rows, int cols, int elements) : GMatrix()
 GSparseMatrix::GSparseMatrix(const GSparseMatrix& m) : GMatrix(m)
 { 
   // Initialise private members for clean destruction
-  m_rowinx   = NULL;
-  m_symbolic = NULL;
-  m_numeric  = NULL;
-
-  // Copy data members
-  m_mem_block = m.m_mem_block;
-  m_zero      = m.m_zero;
-  m_fill_val  = m.m_fill_val;
-  m_fill_row  = m.m_fill_row;
-  m_fill_col  = m.m_fill_col;
-  
-  // Allocate memory for row indices and copy them (only if there are elements)
-  if (m_alloc > 0) {
-    m_rowinx = new int[m_alloc];
-    if (m_rowinx == NULL)
-	  throw mem_alloc("GSparseMatrix copy constructor", m_alloc);
-    for (int i = 0; i < m_elements; ++i)
-      m_rowinx[i] = m.m_rowinx[i];
-  }
+  init_members();
-  // If there is a symbolic decomposition then copy it
-  if (m.m_symbolic != NULL) {
-    GSparseSymbolic* symbolic = new GSparseSymbolic();
-    if (symbolic == NULL)
-	  throw mem_alloc("GSparseMatrix copy constructor", 1);
-    *symbolic  = *((GSparseSymbolic*)m.m_symbolic);
-    m_symbolic = (void*)symbolic;
-  }
+  // Copy members
+  copy_members(m);
-  // If there is a numeric decomposition then copy it
-  if (m.m_numeric != NULL) {
-    GSparseNumeric* numeric = new GSparseNumeric();
-    if (numeric == NULL)
-	  throw mem_alloc("GSparseMatrix copy constructor", 1);
-    *numeric  = *((GSparseNumeric*)m.m_numeric);
-    m_numeric = (void*)numeric;
-  }
-  
   // Return
   return;
 }
@@ -147,10 +208,8 @@ GSparseMatrix::GSparseMatrix(const GSparseMatrix&amp; m) : GMatrix(m)
  ***************************************************************************/
 GSparseMatrix::~GSparseMatrix()
 {
-  // De-allocate only if memory has indeed been allocated by derived class
-  if (m_rowinx   != NULL) delete [] m_rowinx;
-  if (m_symbolic != NULL) delete (GSparseSymbolic*)m_symbolic;
-  if (m_numeric  != NULL) delete (GSparseNumeric*)m_numeric;
+  // Free members
+  free_members();
   // Return
   return;
@@ -174,50 +233,14 @@ GSparseMatrix&amp; GSparseMatrix::operator= (const GSparseMatrix&amp; m)
     // Invoque base class assignment operator
     this->GMatrix::operator=(m);
-    // De-allocate memory if it has indeed been allocated
-    if (m_rowinx   != NULL) delete [] m_rowinx;
-    if (m_symbolic != NULL) delete (GSparseSymbolic*)m_symbolic;
-    if (m_numeric  != NULL) delete (GSparseNumeric*)m_numeric;
+    // Free members
+    free_members();
     // Initialise private members for clean destruction
-    m_rowinx   = NULL;
-	m_symbolic = NULL;
-	m_numeric  = NULL;
-
-    // Copy data members
-	m_mem_block = m.m_mem_block;
-    m_zero      = m.m_zero;
-    m_fill_val  = m.m_fill_val;
-    m_fill_row  = m.m_fill_row;
-    m_fill_col  = m.m_fill_col;
-
-    // Allocate memory for row indices and copy them (only if there are elements)
-    if (m_alloc > 0) {
-      m_rowinx = new int[m_alloc];
-      if (m_rowinx == NULL)
-	    throw mem_alloc("GSparseMatrix::operator= (const GSparseMatrix&)", 
-		                m_alloc);
-      for (int i = 0; i < m_elements; ++i)
-        m_rowinx[i] = m.m_rowinx[i];
-    }
+    init_members();
-    // If there is a symbolic decomposition then copy it
-    if (m.m_symbolic != NULL) {
-      GSparseSymbolic* symbolic = new GSparseSymbolic();
-      if (symbolic == NULL)
-	    throw mem_alloc("GSparseMatrix copy constructor", 1);
-      *symbolic  = *((GSparseSymbolic*)m.m_symbolic);
-      m_symbolic = (void*)symbolic;
-    }
-
-    // If there is a numeric decomposition then copy it
-    if (m.m_numeric != NULL) {
-      GSparseNumeric* numeric = new GSparseNumeric();
-      if (numeric == NULL)
-	    throw mem_alloc("GSparseMatrix copy constructor", 1);
-      *numeric  = *((GSparseNumeric*)m.m_numeric);
-      m_numeric = (void*)numeric;
-    }
+    // Copy members
+    copy_members(m);
   } // endif: object was not identical
@@ -742,22 +765,40 @@ void GSparseMatrix::insert_col(const GVector&amp; v, int col)
  ***************************************************************************/
 void GSparseMatrix::add_col(const GVector& v, int col)
 {
-  // Raise an exception if the column index is invalid
-  #if defined(G_RANGE_CHECK)
-  if (col >= m_cols)
-    throw out_of_range("GSparseMatrix::add_col(const GVector&, int)", 
-	                   0, col, m_rows, m_cols);
-  #endif
+  // Initialise number of non-zero elements to 0
+  int non_zero = 0;
+
+  // If we have a stack then try to push elements on stack first ...
+  if (m_stack_data != NULL) {
+    non_zero = stack_push_column(v, col);
+	if (non_zero == 0)
+	  return;
+  }
+   
+  // ... otherwise check first the arguments and determine the number of
+  // non-zero elements in vector
+  else {
-  // Raise an exception if the matrix and vector dimensions are incompatible
-  if (m_rows != v.m_num)
-    throw matrix_vector_mismatch("GSparseMatrix::add_col(const GVector&, int)", 
-	                             v.m_num, m_rows, m_cols);
+    // Raise an exception if the column index is invalid
+    #if defined(G_RANGE_CHECK)
+    if (col >= m_cols)
+      throw out_of_range("GSparseMatrix::add_col(const GVector&, int)", 
+	                     0, col, m_rows, m_cols);
+    #endif
+
+    // Raise an exception if the matrix and vector dimensions are incompatible
+    if (m_rows != v.m_num)
+      throw matrix_vector_mismatch("GSparseMatrix::add_col(const GVector&, int)", 
+	                               v.m_num, m_rows, m_cols);
+
+    // Determine number of non-zero elements in vector
+	non_zero = v.non_zeros();	
+  }
   // Extract vector for column, add elements, and re-insert vector (only if
   // vector to insert has non-zeros)
-  if (v.non_zeros() > 0) {
-  
+  if (non_zero > 0) {
+
     // Copy input vector
     GVector column = v;
@@ -771,9 +812,10 @@ void GSparseMatrix::add_col(const GVector&amp; v, int col)
 	// Insert vector into matrix
     insert_col(column, col);
+	
   }
-
+  
   // Return
   return;
 }
@@ -1119,6 +1161,405 @@ void GSparseMatrix::cholesky_invert(int compress)
 }
+/***************************************************************************
+ *                     Initialises matrix filling stack                    *
+ * ----------------------------------------------------------------------- *
+ * The matrix filling stack is used to allow for a fast column-wise        *
+ * filling of a sparse matrix. Columns are successively appended to a      *
+ * stack which is regularily flushed when it is full. This reduces memory  *
+ * copies and movements and increases filling speed.                       *
+ ***************************************************************************/
+void GSparseMatrix::stack_init(int size, int entries)
+{
+  // If there is an existing stack the delete it now
+  if (m_stack_colinx != NULL) delete [] m_stack_colinx;
+  if (m_stack_start  != NULL) delete [] m_stack_start;
+  if (m_stack_data   != NULL) delete [] m_stack_data;
+  if (m_stack_rowinx != NULL) delete [] m_stack_rowinx;
+  
+  // Initialise stack members
+  m_stack_max_entries = (entries > 0) ? entries : G_SPARSE_MATRIX_DEFAULT_STACK_ENTRIES;
+  m_stack_size        = (size > 0) ? size : G_SPARSE_MATRIX_DEFAULT_STACK_SIZE;
+  m_stack_entries     = 0;
+  m_stack_colinx      = NULL;
+  m_stack_start       = NULL;
+  m_stack_data        = NULL;
+  m_stack_rowinx      = NULL;
+  
+  // Allocate stack memory. Raise an exception if allocation fails
+  m_stack_colinx = new int[m_stack_max_entries];
+  m_stack_start  = new int[m_stack_max_entries+1];
+  m_stack_data   = new double[m_stack_size];
+  m_stack_rowinx = new int[m_stack_size];
+  if (m_stack_colinx == NULL)
+	throw mem_alloc("GSparseMatrix::stack_init", m_stack_max_entries);
+  if (m_stack_start == NULL)
+	throw mem_alloc("GSparseMatrix::stack_init", m_stack_max_entries+1);
+  if (m_stack_data == NULL || m_stack_rowinx == NULL)
+	throw mem_alloc("GSparseMatrix::stack_init", m_stack_size);
+	
+  // Initialise next free stack location to the first stack element
+  m_stack_start[0] = 0;
+
+  // Return
+  return;
+}
+
+
+/***************************************************************************
+ *                     Push a column on the matrix stack                   *
+ * ----------------------------------------------------------------------- *
+ * The return value indicates if the vector was successfully pushed onto   *
+ * the stack (=0) or if the elements need to be added directly (in this    *
+ * case the return value equals to the number of non-zero elements in the  *
+ * vector.                                                                 *
+ ***************************************************************************/
+int GSparseMatrix::stack_push_column(const GVector& v, int col)
+{
+  // Raise an exception if the column index is invalid
+  #if defined(G_RANGE_CHECK)
+  if (col >= m_cols)
+    throw out_of_range("GSparseMatrix::stack_push_column(const GVector&, int)", 
+	                   0, col, m_rows, m_cols);
+  #endif
+
+  // Raise an exception if the matrix and vector dimensions are incompatible
+  if (m_rows != v.m_num)
+    throw matrix_vector_mismatch("GSparseMatrix::stack_push_column(const GVector&, int)", 
+	                             v.m_num, m_rows, m_cols);
+
+  // Determines how many elements should be stored in stack
+  int elements = 0;
+  for (int i = 0; i < v.m_num; ++i) {
+    if (v.m_data[i] != 0.0)
+	  elements++;
+  }
+  
+  // If the vector is empty there is nothing to do
+  if (elements == 0)
+    return 0;
+	
+  // If there is no stack or the stack can not hold the requested number of 
+  // elements then report number of vector elements to the caller
+  if (m_stack_data == NULL || elements > m_stack_size)
+    return elements;
+
+  // If the specified column is already in the stack then flush stack (we don't
+  // allow the same column twice)
+  for (int i = 0; i < m_stack_entries; ++i) {
+    if (col == m_stack_colinx[i]) {
+      stack_flush();
+	  break;
+	}
+  }
+	
+  // If the stack is full then flush it before pushing the vector onto it
+  if ((m_stack_entries >= m_stack_max_entries) ||
+      (elements        >= (m_stack_size - m_stack_start[m_stack_entries])))
+	stack_flush();
+
+  // Push the vector on the stack
+  int inx = m_stack_start[m_stack_entries];
+  for (int i = 0; i < v.m_num; ++i) {
+    if (v.m_data[i] != 0.0) {
+	  m_stack_data[inx]   = v.m_data[i];
+	  m_stack_rowinx[inx] = i;
+	  inx++;
+	}
+  }
+  
+  // Store entry information and initialise start of next entry
+  m_stack_colinx[m_stack_entries] = col;     // Store column index for entry
+  m_stack_entries++;                         // Increase the # of entries
+  m_stack_start[m_stack_entries] = inx;      // Set start pointer for next entry
+
+  // Debug: show stack information  
+  #if defined(G_DEBUG_SPARSE_STACK_PUSH)
+  cout << "GSparseMatrix::stack_push_column(v, col=" << col << ")" << endl;
+  cout << " Number of entries: " << m_stack_entries << endl;
+  cout << " Entry columns ...:";
+  for (int i = 0; i < m_stack_entries; ++i)
+    cout << " " << m_stack_colinx[i];
+  cout << endl;
+  cout << " Entry starts ....:";
+  for (int i = 0; i < m_stack_entries; ++i)
+    cout << " " << m_stack_start[i];
+  cout << " (next at " << m_stack_start[m_stack_entries] << ")" << endl;
+  cout << " Stack data ......:";
+  for (int i = 0; i < m_stack_start[m_stack_entries]; ++i)
+    cout << " " << m_stack_data[i];
+  cout << endl;
+  cout << " Stack rows ......:";
+  for (int i = 0; i < m_stack_start[m_stack_entries]; ++i)
+    cout << " " << m_stack_rowinx[i];
+  cout << endl;
+  #endif
+	  
+  // Return success
+  return 0;
+}
+
+
+/***************************************************************************
+ *                              Flush matrix stack                         *
+ * ----------------------------------------------------------------------- *
+ * Adds the stack to the actual matrix. First, the total number of matrix  *
+ * elements is determined. The new memory is allocated to hold all         *
+ * elements. Finally, all elements are filled into a working array that is *
+ * then compressed into the matrix.                                        *
+ * ----------------------------------------------------------------------- *
+ * NOTE: The present algorithm assumes that each column occurs only once   *
+ * in the stack!                                                           *
+ ***************************************************************************/
+void GSparseMatrix::stack_flush(void)
+{
+  // Do nothing if there is no stack
+  if (m_stack_data == NULL)
+    return;
+
+  // Fill pending value
+  fill_pending();
+    
+  // Debug header 
+  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
+  cout << "GSparseMatrix::stack_flush" << endl;
+  cout << " Stack entries ...: " << m_stack_entries << endl;
+  cout << " Stack elements ..: " << m_stack_start[m_stack_entries] << endl;
+  cout << " Matrix elements .: " << m_elements << endl;
+  #endif
+ 
+  // Allocate workspace for column flags
+  int* wrk_int = new int[m_cols];
+  if (wrk_int == NULL)
+	throw mem_alloc("GSparseMatrix::stack_flush", m_cols);
+
+  // First gather all columns from existing matrix and flag all existing
+  // columns by 1 and all non-existing by 0
+  for (int col = 0; col < m_cols; ++col) {
+    if (m_colstart[col] < m_colstart[col+1])
+	  wrk_int[col] = 1;
+	else
+      wrk_int[col] = 0;
+  }
+  
+  // Loop over all stack entries and set column flags. They have the 
+  // following definition:
+  //  -(entry+2): Column exists and is also found in stack=entry
+  //  +(entry+2): Column exists in stack=entry only
+  //           1: Column exists in matrix only
+  //           0: Column does neither exist in matrix or stack  
+  int new_elements = 0;
+  for (int entry = 0; entry < m_stack_entries; ++entry) {
+  
+    // Get column for actual entry
+    int col = m_stack_colinx[entry];
+
+	// If column exists already in matrix then determine total number
+	// of elements
+	if (wrk_int[col] == 1) {
+	  int i_start = m_colstart[col];
+	  int i_stop  = m_colstart[col+1];
+	  int k_start = m_stack_start[entry];
+	  int k_stop  = m_stack_start[entry+1];
+	  int row_m   = m_rowinx[i_start];
+	  int row_s   = m_stack_rowinx[k_start];
+	  if (k_start < k_stop && i_start < i_stop) {
+	  	wrk_int[col] = -(entry+2);                 // Flag that we have to mix
+	    while (i_start < i_stop || k_start < k_stop) {
+	      if (row_m == row_s) {                    // Case A: element in matrix and stack
+	        i_start++;
+	        k_start++;
+	        row_m = (i_start < i_stop) ? m_rowinx[i_start]       : m_rows+1;
+	        row_s = (k_start < k_stop) ? m_stack_rowinx[k_start] : m_rows+1;
+		  }
+		  else if (row_m < row_s) {                // Case B: element in matrix only
+	        i_start++;
+	        row_m = (i_start < i_stop) ? m_rowinx[i_start] : m_rows+1;
+		  }
+		  else {                                   // Case C: element in stack only
+	        k_start++;
+	        row_s = (k_start < k_stop) ? m_stack_rowinx[k_start] : m_rows+1;
+		    new_elements++;
+		  }
+        } // endwhile: mixing
+	  } // endif: there were elements in stack
+	} // endif: column existed in the matrix
+
+	// Column did not exists in the matrix, so it's newand we add all elements 
+	else {
+	  wrk_int[col]  = (entry+2);
+	  new_elements += (m_stack_start[entry+1] - m_stack_start[entry]);
+	}
+  }
+  int elements = m_elements + new_elements;
+
+  // Dump number of elements in new matrix 
+  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
+  cout << " New elements ....: " << new_elements << endl;
+  #endif
+  
+  // Allocate memory for new matrix (always keep some elbow room)
+  m_alloc = elements + m_mem_block;
+  double* new_data   = new double[m_alloc];
+  int*    new_rowinx = new int[m_alloc];
+  if (new_data == NULL || new_rowinx == NULL)
+	throw mem_alloc("GSparseMatrix::stack_flush", m_alloc);
+  
+  // Fill new matrix. For this purpose we loop over all matrix columns
+  // and perform the operation that was identified in the previous scan
+  int index = 0;
+  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
+  int num_orig  = 0;
+  int num_stack = 0;
+  int num_mix   = 0;
+  #endif
+  for (int col = 0; col < m_cols; ++col) {
+
+    // If column does not exist then skip
+	if (wrk_int[col] == 0) {
+	  m_colstart[col] = index;
+	  continue;
+	}
+
+    // If column exists in the matrix but not in the stack we copy the
+	// existing column into the new matrix
+	if (wrk_int[col] == 1) {
+	  for (int i = m_colstart[col]; i < m_colstart[col+1]; ++i) {
+	    new_data[index]   = m_data[i];
+		new_rowinx[index] = m_rowinx[i];
+		index++;
+        #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
+		num_orig++;
+		#endif
+	  }
+	}
+	
+	// If column is new (i.e. it did not exist in the matrix before) we 
+	// copy the stack into the new matrix
+	else if (wrk_int[col] > 1) {
+	  int entry = wrk_int[col] - 2;
+	  for (int i = m_stack_start[entry]; i < m_stack_start[entry+1]; ++i) {
+	    new_data[index]   = m_stack_data[i];
+		new_rowinx[index] = m_stack_rowinx[i];
+		index++;
+        #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
+		num_stack++;
+		#endif
+	  }
+	}
+	
+	// If column exists in the matrix and in the stack then we have the
+	// most tricky case since we have to mix both ...
+	else {
+	  int entry    = -(wrk_int[col] + 2);
+	  int i_start  = m_colstart[col];
+	  int i_stop   = m_colstart[col+1];
+	  int k_start  = m_stack_start[entry];
+	  int k_stop   = m_stack_start[entry+1];
+	  int row_m    = m_rowinx[i_start];
+	  int row_s    = m_stack_rowinx[k_start];
+	  if (k_start < k_stop && i_start < i_stop) {
+	    while (i_start < i_stop || k_start < k_stop) {
+	      if (row_m == row_s) {                    // Case A: element in matrix and stack
+		    new_data[index]    = m_data[i_start] + m_stack_data[k_start];
+		    new_rowinx[index]  = row_m;
+		    i_start++;
+		    k_start++;
+		    row_m = (i_start < i_stop) ? m_rowinx[i_start]       : m_rows+1;
+		    row_s = (k_start < k_stop) ? m_stack_rowinx[k_start] : m_rows+1;
+		  }
+		  else if (row_m < row_s) {                // Case B: element in matrix only
+		    new_data[index]   = m_data[i_start];
+		    new_rowinx[index] = row_m;
+		    i_start++;
+		    row_m = (i_start < i_stop) ? m_rowinx[i_start] : m_rows+1;
+		  }
+		  else {                                   // Case C: element in stack only
+		    new_data[index]   = m_stack_data[k_start];
+		    new_rowinx[index] = row_s;
+		    k_start++;
+		    row_s = (k_start < k_stop) ? m_stack_rowinx[k_start] : m_rows+1;
+		  } 
+		  index++;
+          #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
+		  num_mix++;
+		  #endif
+        } // endwhile: mixing
+	  } // endif: there were elements in stack
+	  
+	} // endelse: we had to mix the entries
+
+	// Store actual index in column start array
+	m_colstart[col] = index;
+	
+  } // endfor: looped over all columns
+
+  // Dump number of elements in new matrix after addition
+  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
+  cout << " Added elements ..: " << index << " (sould be " << elements << ")" << endl;
+  cout << " - Matrix only ...: " << num_orig << endl;
+  cout << " - Stack only ....: " << num_stack << endl;
+  cout << " - Matrix & Stack : " << num_mix << endl;
+  #endif
+
+  // Correct columns start array
+  for (int col = m_cols; col > 0; --col)
+    m_colstart[col] = m_colstart[col-1];
+  m_colstart[0] = 0;
+  
+  // Delete old matrix memory
+  if (m_data   != NULL) delete [] m_data;
+  if (m_rowinx != NULL) delete [] m_rowinx;
+	
+  // Update pointers to new memory and update element counter
+  m_data     = new_data;
+  m_rowinx   = new_rowinx;
+  m_elements = elements;
+
+  // Stack is empty now, so reset stack counters
+  m_stack_entries  = 0;
+  m_stack_start[0] = 0;
+
+  // Free workspace
+  delete [] wrk_int;
+  
+  // Return
+  return;
+}
+
+
+/***************************************************************************
+ *                             Destroy matrix stack                        *
+ * ----------------------------------------------------------------------- *
+ * Flush and destroy matrix stack                                          *
+ ***************************************************************************/
+void GSparseMatrix::stack_destroy(void)
+{
+  // Flush stack first
+  stack_flush();
+  
+  // Free stack members
+  if (m_stack_colinx != NULL) delete [] m_stack_colinx;
+  if (m_stack_start  != NULL) delete [] m_stack_start;
+  if (m_stack_data   != NULL) delete [] m_stack_data;
+  if (m_stack_rowinx != NULL) delete [] m_stack_rowinx;
+
+  // Properly mark members as free
+  m_stack_colinx = NULL;
+  m_stack_start  = NULL;
+  m_stack_data   = NULL;
+  m_stack_rowinx = NULL;
+
+  // Initialise stack members
+  m_stack_max_entries = 0;
+  m_stack_size        = 0;
+  m_stack_entries     = 0;
+  
+  // Return
+  return;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                     GSparseMatrix private functions                     =
@@ -24,6 +24,8 @@ using namespace std;
 /* __ Definitions ________________________________________________________ */
 #define G_SPARSE_MATRIX_DEFAULT_MEM_BLOCK 512    // Default Memory block size
+#define G_SPARSE_MATRIX_DEFAULT_STACK_ENTRIES 100   // Max # of stack entries
+#define G_SPARSE_MATRIX_DEFAULT_STACK_SIZE 512     // Max # of stack elements
 /* __ Macros (to be moved to GammaLib general header) ____________________ */
 #define G_MIN(a,b) (((a) < (b)) ? (a) : (b))
@@ -104,25 +106,43 @@ public:
   virtual void    cholesky_decompose(int compress = 1);
   virtual GVector cholesky_solver(const GVector& v, int compress = 1);
   virtual void    cholesky_invert(int compress = 1);
+  
+  // Stack functions
+  virtual void stack_init(int size = 0, int entries = 0);
+  virtual int  stack_push_column(const GVector& v, int col);
+  virtual void stack_flush(void);
+  virtual void stack_destroy(void);
 private:
   // Private functions
-  int           get_index(int row, int col) const;
-  void          fill_pending(void);
-  void          alloc_elements(int start, int num);
-  void          free_elements(int start, int num);
-  void          remove_zero_row_col(void);
-  void          insert_zero_row_col(int rows, int cols);
+  void init_members(void);
+  void copy_members(const GSparseMatrix& m);
+  void free_members(void);
+  int  get_index(int row, int col) const;
+  void fill_pending(void);
+  void alloc_elements(int start, int num);
+  void free_elements(int start, int num);
+  void remove_zero_row_col(void);
+  void insert_zero_row_col(int rows, int cols);
   // Private data members
-  int*   m_rowinx;      // Row-indices of all elements
-  int    m_mem_block;   // Memory block to be allocated at once
-  double m_zero;        // The zero element (needed for data access)
-  double m_fill_val;    // Element to be filled
-  int    m_fill_row;    // Row into which element needs to be filled
-  int    m_fill_col;    // Column into which element needs to be filled
-  void*  m_symbolic;    // Holds GSparseSymbolic object after decomposition
-  void*  m_numeric;     // Holds GSparseNumeric object after decomposition
+  int*    m_rowinx;             // Row-indices of all elements
+  int     m_mem_block;          // Memory block to be allocated at once
+  double  m_zero;               // The zero element (needed for data access)
+  double  m_fill_val;           // Element to be filled
+  int     m_fill_row;           // Row into which element needs to be filled
+  int     m_fill_col;           // Column into which element needs to be filled
+  void*   m_symbolic;           // Holds GSparseSymbolic object after decomposition
+  void*   m_numeric;            // Holds GSparseNumeric object after decomposition
+
+  // Fill-stack
+  int     m_stack_max_entries;  // Maximum number of entries in the stack
+  int     m_stack_size;         // Maximum number of elements in the stack
+  int     m_stack_entries;      // Number of entries in the stack
+  int*    m_stack_colinx;       // Column index for each entry [m_stack_entries]
+  int*    m_stack_start;        // Start in stack for each entry [m_stack_entries+1]
+  double* m_stack_data;         // Stack data [m_stack_size]
+  int*    m_stack_rowinx;       // Srack row indices [m_stack_size]
 };
@@ -205,6 +225,7 @@ GSparseMatrix&amp; GSparseMatrix::operator*= (const double&amp; d)
 {
   fill_pending();
   this->GMatrix::operator*=(d);
+  return *this;
 }
 // Unary scalar division
@@ -213,6 +234,7 @@ GSparseMatrix&amp; GSparseMatrix::operator/= (const double&amp; d)
 {
   fill_pending();
   this->GMatrix::operator/=(d);
+  return *this;
 }
@@ -488,7 +488,7 @@ int* GSparseSymbolic::cs_amd(int order, const GSparseMatrix* A)
   #endif
   // Drop diagonal entries from C
-  int cnz_diag = cs_fkeep(&C, &cs_diag, NULL);
+  cs_fkeep(&C, &cs_diag, NULL);
   // Debug
   #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
@@ -389,6 +389,7 @@ void GSymMatrix::cholesky_decompose(int compress)
   if (no_zeros) {
     // Loop over upper triangle (col >= row)
+	double diag = 0.0;
     for (int row = 0; row < m_rows; ++row) {
 	  double* ptr = m_data + m_colstart[row];
       for (int col = row; col < m_cols; ++col, ++ptr) {
@@ -397,7 +398,6 @@ void GSymMatrix::cholesky_decompose(int compress)
 	      int offset = m_colstart[k] - k;            // is always positive
 	      sum -= m_data[offset+row] * m_data[offset+col]; // sum -= M(row,k)*M(col,k)
 	    }
-	    double diag;
 	    if (row == col) {
 	      if (sum <= 0.0)
 		    throw matrix_not_pos_definite("GSymMatrix::cholesky_decompose(int)", 
@@ -423,6 +423,7 @@ void GSymMatrix::cholesky_decompose(int compress)
 	int* k_ptr;
     // Loop over upper triangle (col >= row)
+	double diag = 0.0;
     for (row = 0, row_ptr = m_inx; row < m_num_inx; ++row, ++row_ptr) {
 	  double* ptr_0 = m_data + m_colstart[*row_ptr] - *row_ptr;
       for (col = row, col_ptr = m_inx + row; col < m_num_inx; ++col, ++col_ptr) {
@@ -432,7 +433,6 @@ void GSymMatrix::cholesky_decompose(int compress)
 	      int offset = m_colstart[*k_ptr] - *k_ptr;            // is always positive
 	      sum -= m_data[offset+*row_ptr] * m_data[offset+*col_ptr]; // sum -= M(row,k)*M(col,k)
 	    }
-	    double diag;
 	    if (*row_ptr == *col_ptr) {
 	      if (sum <= 0.0)
 		    throw matrix_not_pos_definite("GSymMatrix::cholesky_decompose(int)", 
@@ -159,6 +159,7 @@ inline
 GSymMatrix& GSymMatrix::operator*= (const double& d)
 {
   this->GMatrix::operator*=(d);
+  return *this;
 }
 // Unary scalar division
@@ -166,6 +167,7 @@ inline
 GSymMatrix& GSymMatrix::operator/= (const double& d)
 {
   this->GMatrix::operator/=(d);
+  return *this;
 }
+/***************************************************************************
+ *         GammaLib.hpp  -  Gamma-Ray Astronomy Library Header file        *
+ * ----------------------------------------------------------------------- *
+ *  copyright            : (C) 2006 by Jurgen Knodlseder                   *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef GAMMALIB_HPP
+#define GAMMALIB_HPP
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include "GException.hpp"
+#include "GVector.hpp"
+#include "GMatrix.hpp"
+#include "GSymMatrix.hpp"
+#include "GSparseMatrix.hpp"
+#include "GSparseSymbolic.hpp"
+#include "GSparseNumeric.hpp"
+
+/* __ Globals ____________________________________________________________ */
+
+/* __ Prototypes _________________________________________________________ */
+
+#endif /* GAMMALIB_HPP */
 ## Process this file with automake to produce Makefile.in
-files = GException.hpp GException.cpp \
-        GVector.hpp GVector.cpp \
-	GMatrix.hpp GMatrix.cpp \
-	GSymMatrix.hpp GSymMatrix.cpp \
-	GSparseMatrix.hpp GSparseMatrix.cpp \
-	GSparseSymbolic.hpp GSparseSymbolic.cpp \
-	GSparseNumeric.hpp GSparseNumeric.cpp
+SUBDIRS = src test
-# Build dynamic library
-#lib_LTLIBRARIES = libgamma.la
-#libgamma_la_SOURCES = $(files)
-#libgamma_la_LDFLAGS =
-
-# Build static library
-noinst_LIBRARIES = libgamma.a
-libgamma_a_SOURCES = $(files)
-#libgamma_a_LDFLAGS =
 ## Process this file with automake to produce Makefile.in
-check_PROGRAMS = test_GVector \
-                 test_GMatrix \
-                 test_GSymMatrix \
-		 test_GSparseMatrix
-TESTS = test_GVector \
-        test_GMatrix \
-        test_GSymMatrix \
-	test_GSparseMatrix
-test_GVector_SOURCES = test_GVector.cpp
-#test_GVector_LDADD = ../src/libgamma.la
-test_GVector_LDADD = ../src/libgamma.a
-test_GMatrix_SOURCES = test_GMatrix.cpp
-#test_GMatrix_LDADD = ../src/libgamma.la
-test_GMatrix_LDADD = ../src/libgamma.a
-test_GSymMatrix_SOURCES = test_GSymMatrix.cpp
-#test_GSymMatrix_LDADD = ../src/libgamma.la
-test_GSymMatrix_LDADD = ../src/libgamma.a
-test_GSparseMatrix_SOURCES = test_GSparseMatrix.cpp
-#test_GSparseMatrix_LDADD = ../src/libgamma.la
-test_GSparseMatrix_LDADD = ../src/libgamma.a
+SUBDIRS = src test
+
@@ -183,7 +183,7 @@ test_GSymMatrix_SOURCES = test_GSymMatrix.cpp
 test_GSymMatrix_LDADD = ../src/libgamma.a
 test_GSparseMatrix_SOURCES = test_GSparseMatrix.cpp
 #test_GSparseMatrix_LDADD = ../src/libgamma.la
-test_GSparseMatrix_LDADD = ../src/libgamma.a
+test_GSparseMatrix_LDADD = ../src/libgamma.a -lc_p
 all: all-am
 .SUFFIXES:
@@ -22,8 +22,8 @@ double   g_matrix[] = {1.0, 2.0, 3.0, 4.0,
                        5.0, 6.0, 7.0, 8.0, 
 					   9.0, 10., 11., 12.};
 double   g_vector[] = {1.0, 2.0, 3.0, 4.0};
-unsigned g_rows    = 3;
-unsigned g_cols    = 4;
+int g_rows    = 3;
+int g_cols    = 4;
 /***************************************************************************
@@ -33,8 +33,8 @@ GMatrix set_matrix(void)
 {
   try {
     GMatrix matrix(g_rows,g_cols);
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col)
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col)
 	    matrix(row,col) = g_matrix[col+row*g_cols];
 	}
     return matrix;
@@ -53,7 +53,7 @@ GVector set_vector(void)
 {
   try {
     GVector vector(g_cols);
-	for (unsigned col = 0; col < g_cols; ++col)
+	for (int col = 0; col < g_cols; ++col)
 	  vector(col) = g_vector[col];
     return vector;
   }
@@ -71,8 +71,8 @@ int check_matrix(const GMatrix&amp; m, const double scale, const double add)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col) {
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col) {
 	    double value = g_matrix[col+row*g_cols] * scale + add;
 	    if (fabs(m(row,col)-value) > 1.0e-15) {
 		  result = 0;
@@ -96,8 +96,8 @@ int check_transpose_matrix(const GMatrix&amp; m, const double scale, const double ad
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col) {
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col) {
 	    double value = g_matrix[col+row*g_cols] * scale + add;
 	    if (m(col,row) != value) {
 		  result = 0;
@@ -121,8 +121,8 @@ int check_matrix_lt(const GMatrix&amp; m, const double scale, const double add)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col) {
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col) {
 	    double value = (col <= row) ? g_matrix[col+row*g_cols] * scale + add : 0.0;
 	    if (m(row,col) != value) {
 		  result = 0;
@@ -146,8 +146,8 @@ int check_matrix_ut(const GMatrix&amp; m, const double scale, const double add)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col) {
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col) {
 	    double value = (col >= row) ? g_matrix[col+row*g_cols] * scale + add : 0.0;
 	    if (m(row,col) != value) {
 		  result = 0;
@@ -171,9 +171,9 @@ int check_matrix_vector(const GVector&amp; v)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < v.size(); ++row) {
+    for (int row = 0; row < v.size(); ++row) {
       double value = 0.0;
-	  for (unsigned col = 0; col < g_cols; ++col)
+	  for (int col = 0; col < g_cols; ++col)
 	    value += g_matrix[col+row*g_cols] * g_vector[col];
 	  if (v(row) != value) {
 	    result = 0;
@@ -198,10 +198,10 @@ int check_matrix_matrix(const GMatrix&amp; m)
     return 0;
   int result = 1;
   try {
-    for (unsigned row = 0; row < m.rows(); ++row) {
-	  for (unsigned col = 0; col < m.cols(); ++col) {
+    for (int row = 0; row < m.rows(); ++row) {
+	  for (int col = 0; col < m.cols(); ++col) {
         double value = 0.0;
-		for (unsigned i = 0; i < g_cols; ++i)
+		for (int i = 0; i < g_cols; ++i)
 	      value += g_matrix[i+row*g_cols] * g_matrix[i+col*g_cols];
 	    if (m(row,col) != value) {
 	      result = 0;
@@ -224,8 +224,8 @@ int check_matrix_matrix(const GMatrix&amp; m)
 int check_matrix_min(const double min)
 {
   double value = g_matrix[0];
-  for (unsigned row = 0; row < g_rows; ++row) {
-    for (unsigned col = 0; col < g_cols; ++col) {
+  for (int row = 0; row < g_rows; ++row) {
+    for (int col = 0; col < g_cols; ++col) {
 	  if (g_matrix[col+row*g_cols] < value)
 		value = g_matrix[col+row*g_cols];
 	}
@@ -240,8 +240,8 @@ int check_matrix_min(const double min)
 int check_matrix_max(const double max)
 {
   double value = g_matrix[0];
-  for (unsigned row = 0; row < g_rows; ++row) {
-    for (unsigned col = 0; col < g_cols; ++col) {
+  for (int row = 0; row < g_rows; ++row) {
+    for (int col = 0; col < g_cols; ++col) {
 	  if (g_matrix[col+row*g_cols] > value)
 		value = g_matrix[col+row*g_cols];
 	}
@@ -256,8 +256,8 @@ int check_matrix_max(const double max)
 int check_matrix_sum(const double sum)
 {
   double value = 0.0;
-  for (unsigned row = 0; row < g_rows; ++row) {
-    for (unsigned col = 0; col < g_cols; ++col)
+  for (int row = 0; row < g_rows; ++row) {
+    for (int col = 0; col < g_cols; ++col)
       value += g_matrix[col+row*g_cols];
   }
   return (sum == value);
@@ -674,4 +674,6 @@ int main(void)
   }
   cout << "ok." << endl;
+  // Return
+  return 0;
 }
@@ -23,12 +23,12 @@ using namespace std;
 /* __ Globals ____________________________________________________________ */
 double   g_matrix[] = {1.0, 7.0, 3.0, 2.0, 4.0, 8.0, 5.0, 6.0, 9.0};
-unsigned g_row[]    = {  0,   0,   1,   2,   2,   2,   3,   3,   3};
-unsigned g_col[]    = {  0,   4,   1,   0,   2,   4,   2,   3,   4};
+int g_row[]    = {  0,   0,   1,   2,   2,   2,   3,   3,   3};
+int g_col[]    = {  0,   4,   1,   0,   2,   4,   2,   3,   4};
 double   g_vector[] = {1.0, 2.0, 3.0, 4.0, 5.0};
-unsigned g_elements = 9;
-unsigned g_rows     = 4;
-unsigned g_cols     = 5;
+int g_elements = 9;
+int g_rows     = 4;
+int g_cols     = 5;
 /***************************************************************************
@@ -38,7 +38,7 @@ GSparseMatrix set_matrix(void)
 {
   try {
     GSparseMatrix matrix(g_rows,g_cols);
-    for (unsigned i = 0; i < g_elements; ++i)
+    for (int i = 0; i < g_elements; ++i)
 	  matrix(g_row[i],g_col[i]) = g_matrix[i];
     return matrix;
   }
@@ -56,7 +56,7 @@ GVector set_vector(void)
 {
   try {
     GVector vector(g_cols);
-	for (unsigned col = 0; col < g_cols; ++col)
+	for (int col = 0; col < g_cols; ++col)
 	  vector(col) = g_vector[col];
     return vector;
   }
@@ -74,9 +74,9 @@ int check_matrix(const GSparseMatrix&amp; m, const double scale, const double add)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col) {
-	    unsigned i;
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col) {
+	    int i;
 		double   ref_value = add;
 		for (i = 0; i < g_elements; ++i) {
 		  if (g_row[i] == row && g_col[i] == col) {
@@ -107,10 +107,10 @@ int check_matrix_vector(const GVector&amp; v)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < v.size(); ++row) {
+    for (int row = 0; row < v.size(); ++row) {
       double value = 0.0;
-	  for (unsigned col = 0; col < g_cols; ++col) {
-	    unsigned i;
+	  for (int col = 0; col < g_cols; ++col) {
+	    int i;
 		double   ref_value = 0.0;
 		for (i = 0; i < g_elements; ++i) {
 		  if (g_row[i] == row && g_col[i] == col) {
@@ -143,11 +143,11 @@ int check_matrix_matrix(const GSparseMatrix&amp; m)
     return 0;
   int result = 1;
   try {
-    for (unsigned row = 0; row < m.rows(); ++row) {
-	  for (unsigned col = 0; col < m.cols(); ++col) {
+    for (int row = 0; row < m.rows(); ++row) {
+	  for (int col = 0; col < m.cols(); ++col) {
         double value = 0.0;
-		for (unsigned i = 0; i < g_cols; ++i) {
-	      unsigned k;
+		for (int i = 0; i < g_cols; ++i) {
+	      int k;
 		  double   ref_value_1 = 0.0;
 		  double   ref_value_2 = 0.0;
 		  for (k = 0; k < g_elements; ++k) {
@@ -187,9 +187,9 @@ int check_transpose_matrix(const GSparseMatrix&amp; m, const double scale, const dou
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col) {
-	    unsigned i;
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col) {
+	    int i;
 		double   ref_value = 0.0;
 		for (i = 0; i < g_elements; ++i) {
 		  if (g_row[i] == row && g_col[i] == col) {
@@ -218,7 +218,7 @@ int check_transpose_matrix(const GSparseMatrix&amp; m, const double scale, const dou
 int check_matrix_min(const double min)
 {
   double value = g_matrix[0];
-  for (unsigned i = 0; i < g_elements; ++i) {
+  for (int i = 0; i < g_elements; ++i) {
 	if (g_matrix[i] < value)
 	  value = g_matrix[i];
   }
@@ -233,7 +233,7 @@ int check_matrix_min(const double min)
 int check_matrix_max(const double max)
 {
   double value = g_matrix[0];
-  for (unsigned i = 0; i < g_elements; ++i) {
+  for (int i = 0; i < g_elements; ++i) {
 	if (g_matrix[i] > value)
 	  value = g_matrix[i];
   }
@@ -247,7 +247,7 @@ int check_matrix_max(const double max)
 int check_matrix_sum(const double sum)
 {
   double value = 0.0;
-  for (unsigned i = 0; i < g_elements; ++i)
+  for (int i = 0; i < g_elements; ++i)
 	  value += g_matrix[i];
   return (sum == value);
 }
@@ -375,7 +375,31 @@ void test_assign_values(const GSparseMatrix&amp; m_test)
 	throw;
   }
   cout << ".";
-  #endif  
+  #endif
+  //
+  // Insert column into 10 x 10 matrix using matrix stack
+  try {
+	GSparseMatrix stack(10,10);
+	stack.clear();
+	for (int i = 3; i < 5; ++i)
+	  stack(i,i) = 5.0;
+	GVector column(10);
+	stack.stack_init(100,50);
+	for (int j = 0; j < 10; ++j) {
+	  int i_min = (j < 2) ?  0 : j-2;
+	  int i_max = (j > 8) ? 10 : j+2;
+	  column = 0.0;
+	  for (int i = i_min; i < i_max; ++i)
+	    column(i) = (i+1)*1;
+	  stack.add_col(column, j);
+    }
+	stack.stack_destroy();
+//	cout << stack << endl;
+  }
+  catch (exception &e) {
+    cout << e.what() << endl;
+	throw;
+  }
   cout << " ok." << endl;
 }
@@ -1346,16 +1370,18 @@ void test_heavy(int block = 512)
   clock_t t_start;
   double  t_elapse;
   int     number;
+  int     num_cols;
+  int     i_min, i_max;
   // Perform test
   cout << "Test GSparseMatrix: Heavy calculus: ";
   try {
     //
 	// Fill 10000 x 10000 matrix with 10000 values
+    t_start = clock();
 	number  = 10000;
 	GSparseMatrix large(number,number);
 	large.set_mem_block(block);
-    t_start = clock();
 	for (int i = 0; i < number; ++i)
 	  large(i,i) = (i+1)*3.14;
 	t_elapse = double(clock() - t_start)/double(CLOCKS_PER_SEC);
@@ -1375,32 +1401,67 @@ void test_heavy(int block = 512)
 	#endif   
     //
 	// Insert columns into 10000 x 10000 matrix
-	number  = 10000;
+	number   = 10000;
+	num_cols = 1000;
 	large.clear();
 	large.set_mem_block(block);
+	for (int j = 200; j < 400; ++j)
+	  large(j,j) = 3.14;
 	GVector column(number);
     t_start = clock();
-	for (int j = 0; j < 100; ++j) {
+	for (int j = 0; j < num_cols; ++j) {
+	  i_min = (j < 2) ? 0 : j-2;
+	  i_max = (j > num_cols-2) ? num_cols : j+2;
 	  column = 0.0;
-	  for (int i = j; i < number; ++i)
-	    column(i) = (i+1)*3.14;
+	  for (int i = i_min; i < i_max; ++i)
+	    column(i) = (i+1)*0.01;
 	  large.add_col(column, j);
     }
 	t_elapse = double(clock() - t_start)/double(CLOCKS_PER_SEC);
     #if defined(DUMP_TIMING)
-	cout << endl << " - 100 columns adding needed " << t_elapse << 
-	                " sec (reference ~ 2.6 sec)";
+	cout << endl << " - " << num_cols << " columns adding needed " << t_elapse << 
+	                " sec (reference ~ 0.4 sec)";
 	#endif   
     //
+	// Insert columns into 10000 x 10000 matrix using matrix stack
+	GSparseMatrix stack(number,number);
+	stack.clear();
+	stack.set_mem_block(block);
+	for (int j = 200; j < 400; ++j)
+	  stack(j,j) = 3.14;
+    t_start = clock();
+	column = GVector(number);
+	stack.stack_init(20000,100);
+	for (int j = 0; j < num_cols; ++j) {
+	  i_min = (j < 2)          ? 0 : j-2;
+	  i_max = (j > num_cols-2) ? num_cols : j+2;
+	  column = 0.0;
+	  for (int i = i_min; i < i_max; ++i)
+	    column(i) = (i+1)*0.01;
+	  stack.add_col(column, j);
+    }
+	stack.stack_destroy();
+	t_elapse = double(clock() - t_start)/double(CLOCKS_PER_SEC);
+    #if defined(DUMP_TIMING)
+	cout << endl << " - " << num_cols << " columns stack-adding needed " << t_elapse << 
+	                " sec (reference ~ 0.1 sec)";
+	#endif
+	//
+	// Check that both are identical
+	if (large != stack) {
+      cout << endl << "TEST ERROR: Stack-based insertion corrupted." << endl;
+	  throw;
+	}
+    //
 	// Subsequent addition and subtractions
 	number  = 300;
+    t_start = clock();
 	GSparseMatrix sum(number,number);
 	GSparseMatrix add(number,number);
 	GSparseMatrix sub(number,number);
 	sum.set_mem_block(block);
 	add.set_mem_block(block);
 	sub.set_mem_block(block);
-    t_start = clock();
 	for (int i = 0; i < number; ++i) {
 	  add(i,number-i-1) = (i*i+1)*0.01;
 	  sum += add;
@@ -1425,7 +1486,7 @@ void test_heavy(int block = 512)
 	t_elapse = double(clock() - t_start)/double(CLOCKS_PER_SEC);
     #if defined(DUMP_TIMING)
 	cout << endl << " - Matrix multiplication 1 times " << t_elapse << 
-	                " sec (reference ~ 3.5 sec)";
+	                " sec (reference ~ 3.1 sec)";
 	#endif   
   }
   catch (exception &e) {
@@ -18,10 +18,10 @@
 using namespace std;
 /* __ Globals ____________________________________________________________ */
-double   g_matrix[] = {4.0, 1.0, 2.0, 1.0, 5.0, 3.0, 2.0, 3.0, 6.0};
-double   g_vector[] = {1.0, 2.0, 3.0};
-unsigned g_rows    = 3;
-unsigned g_cols    = 3;
+double g_matrix[] = {4.0, 1.0, 2.0, 1.0, 5.0, 3.0, 2.0, 3.0, 6.0};
+double g_vector[] = {1.0, 2.0, 3.0};
+int    g_rows    = 3;
+int    g_cols    = 3;
 /***************************************************************************
@@ -31,8 +31,8 @@ GSymMatrix set_matrix(void)
 {
   try {
     GSymMatrix matrix(g_rows,g_cols);
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col)
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col)
 	    matrix(row,col) = g_matrix[col+row*g_cols];
 	}
     return matrix;
@@ -51,11 +51,11 @@ GSymMatrix set_matrix_zero(void)
 {
   try {
     GSymMatrix matrix(g_rows+1,g_cols+1);
-	unsigned i = 0;
-    for (unsigned row = 0; row < g_rows+1; ++row) {
+	int i = 0;
+    for (int row = 0; row < g_rows+1; ++row) {
 	  if (row == 2) continue;
-	  unsigned j = 0;
-      for (unsigned col = 0; col < g_cols+1; ++col) {
+	  int j = 0;
+      for (int col = 0; col < g_cols+1; ++col) {
 	    if (col == 2) continue;
 	    matrix(row,col) = g_matrix[j+i*g_cols];
 		j++;
@@ -78,7 +78,7 @@ GVector set_vector(void)
 {
   try {
     GVector vector(g_cols);
-	for (unsigned col = 0; col < g_cols; ++col)
+	for (int col = 0; col < g_cols; ++col)
 	  vector(col) = g_vector[col];
     return vector;
   }
@@ -96,8 +96,8 @@ int check_matrix(const GSymMatrix&amp; m, const double scale, const double add)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col) {
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col) {
 	    double value = g_matrix[col+row*g_cols] * scale + add;
 	    if (fabs(m(row,col)-value) > 1.0e-15) {
 		  result = 0;
@@ -121,8 +121,8 @@ int check_matrix(const GMatrix&amp; m, const double scale, const double add)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col) {
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col) {
 	    double value = g_matrix[col+row*g_cols] * scale + add;
 	    if (m(row,col) != value) {
 		  result = 0;
@@ -146,8 +146,8 @@ int check_matrix_lt(const GMatrix&amp; m, const double scale, const double add)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col) {
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col) {
 	    double value = (col <= row) ? g_matrix[col+row*g_cols] * scale + add : 0.0;
 	    if (m(row,col) != value) {
 		  result = 0;
@@ -171,8 +171,8 @@ int check_matrix_ut(const GMatrix&amp; m, const double scale, const double add)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < g_rows; ++row) {
-      for (unsigned col = 0; col < g_cols; ++col) {
+    for (int row = 0; row < g_rows; ++row) {
+      for (int col = 0; col < g_cols; ++col) {
 	    double value = (col >= row) ? g_matrix[col+row*g_cols] * scale + add : 0.0;
 	    if (m(row,col) != value) {
 		  result = 0;
@@ -196,9 +196,9 @@ int check_matrix_vector(const GVector&amp; v)
 {
   int result = 1;
   try {
-    for (unsigned row = 0; row < v.size(); ++row) {
+    for (int row = 0; row < v.size(); ++row) {
       double value = 0.0;
-	  for (unsigned col = 0; col < g_cols; ++col)
+	  for (int col = 0; col < g_cols; ++col)
 	    value += g_matrix[col+row*g_cols] * g_vector[col];
 	  if (v(row) != value) {
 	    result = 0;
@@ -223,10 +223,10 @@ int check_matrix_matrix(const GSymMatrix&amp; m)
     return 0;
   int result = 1;
   try {
-    for (unsigned row = 0; row < m.rows(); ++row) {
-	  for (unsigned col = 0; col < m.cols(); ++col) {
+    for (int row = 0; row < m.rows(); ++row) {
+	  for (int col = 0; col < m.cols(); ++col) {
         double value = 0.0;
-		for (unsigned i = 0; i < g_cols; ++i)
+		for (int i = 0; i < g_cols; ++i)
 	      value += g_matrix[i+row*g_cols] * g_matrix[i+col*g_cols];
 	    if (m(row,col) != value) {
 	      result = 0;
@@ -249,8 +249,8 @@ int check_matrix_matrix(const GSymMatrix&amp; m)
 int check_matrix_min(const double min)
 {
   double value = g_matrix[0];
-  for (unsigned row = 0; row < g_rows; ++row) {
-    for (unsigned col = 0; col < g_cols; ++col) {
+  for (int row = 0; row < g_rows; ++row) {
+    for (int col = 0; col < g_cols; ++col) {
 	  if (g_matrix[col+row*g_cols] < value)
 		value = g_matrix[col+row*g_cols];
 	}
@@ -265,8 +265,8 @@ int check_matrix_min(const double min)
 int check_matrix_max(const double max)
 {
   double value = g_matrix[0];
-  for (unsigned row = 0; row < g_rows; ++row) {
-    for (unsigned col = 0; col < g_cols; ++col) {
+  for (int row = 0; row < g_rows; ++row) {
+    for (int col = 0; col < g_cols; ++col) {
 	  if (g_matrix[col+row*g_cols] > value)
 		value = g_matrix[col+row*g_cols];
 	}
@@ -281,8 +281,8 @@ int check_matrix_max(const double max)
 int check_matrix_sum(const double sum)
 {
   double value = 0.0;
-  for (unsigned row = 0; row < g_rows; ++row) {
-    for (unsigned col = 0; col < g_cols; ++col)
+  for (int row = 0; row < g_rows; ++row) {
+    for (int col = 0; col < g_cols; ++col)
       value += g_matrix[col+row*g_cols];
   }
   return (sum == value);
@@ -928,4 +928,7 @@ int main(void)
   }
   cout << "ok." << endl;
+  // Return
+  return 0;
+
 }
@@ -25,18 +25,18 @@ int main(void)
   cout << "*************************" << endl;
   // Set parameters
-  unsigned num = 5;
+  int num = 5;
   // Define vectors
   GVector test(num);
   GVector result(num);
   GVector smaller(num+1);
   GVector bigger(num+1);
-  for (unsigned i = 0; i < num; ++i)
+  for (int i = 0; i < num; ++i)
     test(i) = (i+1) * 1.1;
-  for (unsigned i = 0; i < num-1; ++i)
+  for (int i = 0; i < num-1; ++i)
     smaller(i) = (i+1) * 1.1;
-  for (unsigned i = 0; i < num+1; ++i)
+  for (int i = 0; i < num+1; ++i)
     bigger(i) = (i+1) * 1.1;
   // Test 1: Allocate zero vector
@@ -329,4 +329,7 @@ int main(void)
     cout << e.what() << endl;
   }
+
+  // Return
+  return 0;
 }