Implement block memory allocating scheme

knodlseder
1 parent 53893431
Showing 12 changed files with 431 additions and 162 deletions
src/GMatrix.cpp
src/GMatrix.hpp
src/GSparseMatrix.cpp
src/GSparseMatrix.hpp
src/GSymMatrix.cpp
src/GVector.cpp
src/GVector.hpp
src/Makefile.am
test/Makefile.am
test/test_GMatrix.cpp
test/test_GSparseMatrix.cpp
test/test_GSymMatrix.cpp
@@ -33,6 +33,7 @@ GMatrix::GMatrix(void)
   m_rows       = 0;
   m_cols       = 0;
   m_elements   = 0;
+  m_alloc      = 0;
   m_num_rowsel = 0;
   m_num_colsel = 0;
   m_data       = NULL;
@@ -54,6 +55,7 @@ GMatrix::GMatrix(int rows, int cols)
   m_rows       = 0;
   m_cols       = 0;
   m_elements   = 0;
+  m_alloc      = 0;
   m_num_rowsel = 0;
   m_num_colsel = 0;
   m_data       = NULL;
@@ -61,7 +63,7 @@ GMatrix::GMatrix(int rows, int cols)
   m_rowsel     = NULL;
   m_colsel     = NULL;
  
-  // Determine number of physical elements in matrix
+  // Determine number of elements to store in matrix
   int elements = rows*cols;
  
   // Throw exception if requested matrix size is zero
@@ -75,10 +77,11 @@ GMatrix::GMatrix(int rows, int cols)
   if (m_data == NULL || m_colstart == NULL)
 	throw mem_alloc("GMatrix constructor", elements);
  
-  // Store matrix size (logical and physical)
+  // Store matrix size (logical, storage, allocated)
   m_rows     = rows;
   m_cols     = cols;
   m_elements = elements;
+  m_alloc    = elements;
  
   // Set-up column start indices
   m_colstart[0] = 0;
@@ -106,6 +109,7 @@ GMatrix::GMatrix(const GMatrix&amp; m)
   m_rows       = 0;
   m_cols       = 0;
   m_elements   = 0;
+  m_alloc      = 0;
   m_num_rowsel = 0;
   m_num_colsel = 0;
   m_data       = NULL;
@@ -121,10 +125,10 @@ GMatrix::GMatrix(const GMatrix&amp; m)
     m_colstart[i] = m.m_colstart[i];
  
   // Allocate memory for elements and copy them (only if there are elements)
-  if (m.m_elements > 0) {
-    m_data = new double[m.m_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_elements);
+	  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];
   }
@@ -153,6 +157,7 @@ GMatrix::GMatrix(const GMatrix&amp; m)
   m_rows     = m.m_rows;
   m_cols     = m.m_cols;
   m_elements = m.m_elements;
+  m_alloc    = m.m_alloc;
  
   // Return
   return;
@@ -199,6 +204,7 @@ GMatrix&amp; GMatrix::operator= (const GMatrix&amp; m)
     m_rows       = 0;
     m_cols       = 0;
     m_elements   = 0;
+	m_alloc      = 0;
     m_num_rowsel = 0;
     m_num_colsel = 0;
     m_data       = NULL;
@@ -214,10 +220,10 @@ GMatrix&amp; GMatrix::operator= (const GMatrix&amp; m)
       m_colstart[i] = m.m_colstart[i];
  
     // Allocate memory for elements and copy them (only if there are elements)
-    if (m.m_elements > 0) {
-      m_data = new double[m.m_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_elements);
+	    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];
     }
@@ -246,6 +252,8 @@ GMatrix&amp; GMatrix::operator= (const GMatrix&amp; m)
     m_rows     = m.m_rows;
     m_cols     = m.m_cols;
     m_elements = m.m_elements;
+	m_alloc    = m.m_alloc;
+	
   } // endif: object was not identical
  
   // Return this object
@@ -709,6 +717,7 @@ void GMatrix::select_non_zero(void)
   return;
 }
  
+
 /*==========================================================================
  =                                                                         =
  =                             GMatrix friends                             =
@@ -720,8 +729,24 @@ void GMatrix::select_non_zero(void)
  ***************************************************************************/
 ostream& operator<< (ostream& os, const GMatrix& m)
 {
+  // Put header in stream
+  os << "=== GMatrix ===" << endl;
+  if (m.m_rowsel != NULL)
+    os << " Number of rows ............: " << m.m_rows << " (compressed " <<
+	      m.m_num_rowsel << ")" << endl;
+  else
+    os << " Number of rows ............: " << m.m_rows << endl;
+  if (m.m_colsel != NULL)
+    os << " Number of columns .........: " << m.m_cols << " (compressed " <<
+	      m.m_num_colsel << ")" << endl;
+  else
+    os << " Number of columns .........: " << m.m_cols << endl;
+  os << " Number of elements ........: " << m.m_elements << endl;
+  os << " Number of allocated cells .: " << m.m_alloc << endl;
+
   // Loop over all matrix elements
   for (int row = 0; row < m.m_rows; ++row) {
+    os << " ";
     for (int col = 0; col < m.m_cols; ++col) {
       os << m(row,col);
 	  if (col != m.m_cols-1)
@@ -730,6 +755,20 @@ ostream&amp; operator&lt;&lt; (ostream&amp; os, const GMatrix&amp; m)
 	if (row != m.m_rows-1)
 	  os << endl;
   }
+
+  // If there is a row compression the show scheme
+  if (m.m_rowsel != NULL) {
+    os << endl << " Row selection ..:";
+    for (int row = 0; row < m.m_num_rowsel; ++row)
+      os << " " << m.m_rowsel[row];
+  }
+
+  // If there is a column compression the show scheme
+  if (m.m_colsel != NULL) {
+    os << endl << " Column selection:";
+    for (int col = 0; col < m.m_num_colsel; ++col)
+      os << " " << m.m_colsel[col];
+  }
  
   // Return output stream
   return os;
@@ -142,7 +142,8 @@ protected:
   // Private data area
   int     m_rows;         // Number of rows
   int     m_cols;         // Number of columns
-  int     m_elements;     // Number of (physical) elements in matrix
+  int     m_elements;     // Number of elements stored in matrix
+  int     m_alloc;        // Size of allocated matrix memory
   int     m_num_rowsel;   // Number of selected rows (for compressed decomposition)
   int     m_num_colsel;   // Number of selected columns (for compressed decomposition)
   int*    m_colstart;     // Column start indices (m_cols+1)
@@ -31,7 +31,8 @@ using namespace std;
 //#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_USE_MEMCPY
  
 /*==========================================================================
  =                                                                         =
@@ -45,18 +46,20 @@ using namespace std;
 GSparseMatrix::GSparseMatrix(int rows, int cols, int elements) : GMatrix()
 {
   // Initialise private members for clean destruction
-  m_rows     = 0;
-  m_cols     = 0;
-  m_elements = 0;
-  m_data     = NULL;
-  m_colstart = NULL;
-  m_rowinx   = NULL;
-  m_zero     = 0.0;
-  m_fill_val = 0.0;
-  m_fill_row = 0;
-  m_fill_col = 0;
-  m_symbolic = NULL;
-  m_numeric  = NULL;
+  m_rows      = 0;
+  m_cols      = 0;
+  m_elements  = 0;
+  m_alloc     = 0;
+  m_data      = NULL;
+  m_colstart  = NULL;
+  m_rowinx    = NULL;
+  m_mem_block = G_SPARSE_MATRIX_DEFAULT_MEM_BLOCK;
+  m_zero      = 0.0;
+  m_fill_val  = 0.0;
+  m_fill_row  = 0;
+  m_fill_col  = 0;
+  m_symbolic  = NULL;
+  m_numeric   = NULL;
  
   // Allocate column start array. This is the only array that we can
   // allocate at this time. The other arrays can only be allocated during
@@ -98,16 +101,17 @@ GSparseMatrix::GSparseMatrix(const GSparseMatrix&amp; m) : GMatrix(m)
   m_numeric  = NULL;
  
   // Copy data members
-  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;
+  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_elements > 0) {
-    m_rowinx = new int[m_elements];
+  if (m_alloc > 0) {
+    m_rowinx = new int[m_alloc];
     if (m_rowinx == NULL)
-	  throw mem_alloc("GSparseMatrix copy constructor", m_elements);
+	  throw mem_alloc("GSparseMatrix copy constructor", m_alloc);
     for (int i = 0; i < m_elements; ++i)
       m_rowinx[i] = m.m_rowinx[i];
   }
@@ -181,17 +185,18 @@ GSparseMatrix&amp; GSparseMatrix::operator= (const GSparseMatrix&amp; m)
 	m_numeric  = NULL;
  
     // Copy data members
-    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;
+	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_elements > 0) {
-      m_rowinx = new int[m_elements];
+    if (m_alloc > 0) {
+      m_rowinx = new int[m_alloc];
       if (m_rowinx == NULL)
 	    throw mem_alloc("GSparseMatrix::operator= (const GSparseMatrix&)", 
-		                m_elements);
+		                m_alloc);
       for (int i = 0; i < m_elements; ++i)
         m_rowinx[i] = m.m_rowinx[i];
     }
@@ -444,6 +449,23 @@ GSparseMatrix&amp; GSparseMatrix::operator*= (const GSparseMatrix&amp; m)
  ==========================================================================*/
  
 /***************************************************************************
+ *                       Set all matrix elements to 0                      *
+ ***************************************************************************/
+void GSparseMatrix::clear()
+{
+  // Free memory
+  free_elements(0, m_elements);
+
+  // Initialise column start indices to 0
+  for (int col = 0; col <= m_cols; ++col)
+    m_colstart[col] = 0;
+  
+  // Return
+  return;
+}
+
+
+/***************************************************************************
  *                        Get minimum matrix element                       *
  ***************************************************************************/
 double GSparseMatrix::min() const
@@ -593,7 +615,6 @@ GVector GSparseMatrix::extract_col(int col) const
  
   // Return vector
   return result;
-
 }
  
  
@@ -645,7 +666,7 @@ void GSparseMatrix::insert_col(const GVector&amp; v, int col)
   // Determine the number of non-zero elements in the vector
   int n_vector = 0;
   for (int row = 0; row < m_rows; ++row) {
-    if (v(row) != 0.0)
+    if (v.m_data[row] != 0.0)
 	  n_vector++;
   }
  
@@ -678,8 +699,8 @@ void GSparseMatrix::insert_col(const GVector&amp; v, int col)
   // Insert the vector elements in the matrix
   if (n_vector > 0) {
     for (int row = 0, i = i_start; row < m_rows; ++row) {
-	  if (v(row) != 0.0) {
-	    m_data[i]   = v(row);
+	  if (v.m_data[row] != 0.0) {
+	    m_data[i]   = v.m_data[row];
 	    m_rowinx[i] = row;
 	    i++;
 	  }
@@ -726,11 +747,26 @@ void GSparseMatrix::add_col(const GVector&amp; v, int col)
     throw matrix_vector_mismatch("GSparseMatrix::add_col(const GVector&, int)", 
 	                             v.m_num, m_rows, m_cols);
  
-  // Extract vector for column, add elements, and re-insert vector
-  GVector column = extract_col(col);
-  column += v;
-  insert_col(column, col);
+  // Extract vector for column, add elements, and re-insert vector (only if
+  // vector to insert has non-zeros)
+  if (v.non_zeros() > 0) {
  
+    // Copy input vector
+    GVector column = v;
+
+    // Add elements to vector
+    for (int i = m_colstart[col]; i < m_colstart[col+1]; ++i)
+      column.m_data[m_rowinx[i]] += m_data[i];
+
+    // If there is a pending element then put it in the vector
+    if (m_fill_val != 0.0 && m_fill_col == col)
+      column.m_data[m_fill_row] += m_fill_val;
+
+	// Insert vector into matrix
+    insert_col(column, col);
+	  
+  }
+
   // Return
   return;
 }
@@ -1225,49 +1261,115 @@ void GSparseMatrix::fill_pending(void)
  ***************************************************************************/
 void GSparseMatrix::alloc_elements(int start, int num)
 {
+  // Dump header
+  #if defined(G_DEBUG_SPARSE_MALLOC)
+  cout << "GSparseMatrix::alloc_elements(start=" << start << ", num=" << 
+          num << ")" << endl;
+  cout << " Before allocation : " << m_elements << " " << m_alloc << endl;
+  #endif
+  
   // Continue only if we need memory
   if (num > 0) {
  
     // If start is after the end the append memory
 	if (start > m_elements)
 	  start = m_elements;
+
+    // Determine the requested new logical size of the matrix
+	int new_size = m_elements + num;
+
+    // Case A: the requested memory is already available, so just move the
+	// data to make space for new elements and initialise the new cells
+	if (new_size <= m_alloc) {
+
+	  // Move up all elements after index to insert
+	  #if defined(G_USE_MEMCPY)
+	  int n_copy = m_elements - start;
+	  memmove(&(m_data[start+num]), &(m_data[start]), sizeof(double)*n_copy);
+	  memmove(&(m_rowinx[start+num]), &(m_rowinx[start]), sizeof(long)*n_copy);
+	  #else
+	  for (int i = m_elements - 1; i >= start; --i) {
+	    m_data[i+num]   = m_data[i];
+	    m_rowinx[i+num] = m_rowinx[i];
+	  }
+	  #endif
+
+	  // Clear new elements (zero content for row 0)
+	  for (int i = start; i < start+num; ++i) {
+	    m_data[i]   = 0.0;
+	    m_rowinx[i] = 0;
+	  }
+
+	  // Update element counter
+	  m_elements += num;
+
+	} // endif: Case A: memory already existed
+	
+	// Case B: more memory is needed, so allocate it, copy the existing
+	// content, and initialise new cells
+	else {
+	
+	  // Make sure that enough memory is allocated
+	  int new_propose = m_alloc + m_mem_block;
+	  m_alloc = (new_size > new_propose) ? new_size : new_propose;
  
-	// Allocate memory for new elements
-	double* new_data   = new double[m_elements+num];
-	int*    new_rowinx = new int[m_elements+num];
-	if (new_data == NULL || new_rowinx == NULL)
-	  throw mem_alloc("GSparseMatrix::alloc_elements(int, int)", 
-	                  m_elements+num);
-
-    // Copy all elements before index to insert
-	for (int i = 0; i < start; ++i) {
-	  new_data[i]   = m_data[i];
-	  new_rowinx[i] = m_rowinx[i];
-	}
+	  // Allocate memory for new elements
+	  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::alloc_elements(int, int)", 
+	                    m_alloc);
+
+      // Copy all elements before index to insert
+	  #if defined(G_USE_MEMCPY)
+	  int n_copy = start;
+	  memcpy(new_data,   m_data,   sizeof(double)*n_copy);
+	  memcpy(new_rowinx, m_rowinx, sizeof(long)*n_copy);
+	  #else
+	  for (int i = 0; i < start; ++i) {
+	    new_data[i]   = m_data[i];
+	    new_rowinx[i] = m_rowinx[i];
+	  }
+	  #endif
  
-	// Clear new elements (zero content for row 0)
-	for (int i = start; i < start+num; ++i) {
-	  new_data[i]   = 0.0;
-	  new_rowinx[i] = 0;
-	}
+	  // Clear new elements (zero content for row 0)
+	  for (int i = start; i < start+num; ++i) {
+	    new_data[i]   = 0.0;
+	    new_rowinx[i] = 0;
+	  }
  
-	// Copy all elements after index to insert
-	for (int i = start; i < m_elements; ++i) {
-	  new_data[i+num]   = m_data[i];
-	  new_rowinx[i+num] = m_rowinx[i];
-	}
+	  // Copy all elements after index to insert
+	  #if defined(G_USE_MEMCPY)
+	  n_copy = m_elements - start;
+	  memcpy(&(new_data[start+num]),   &(m_data[start]),   sizeof(double)*n_copy);
+	  memcpy(&(new_rowinx[start+num]), &(m_rowinx[start]), sizeof(long)*n_copy);
+	  #else
+	  for (int i = start; i < m_elements; ++i) {
+	    new_data[i+num]   = m_data[i];
+	    new_rowinx[i+num] = m_rowinx[i];
+	  }
+	  #endif
  
-    // Delete old memory
-    if (m_data   != NULL) delete [] m_data;
-    if (m_rowinx != NULL) delete [] m_rowinx;
+      // Delete old 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 += num;
+	  // Update pointers to new memory and update element counter
+	  m_data      = new_data;
+	  m_rowinx    = new_rowinx;
+	  m_elements += num;
+	  
+	} // endelse: Case B: more memory was needed
  
   } // endif: needed new memory
  
+  // Dump new memory size
+  #if defined(G_DEBUG_SPARSE_MALLOC)
+  cout << " After allocation .: " << m_elements << " " << m_alloc << endl;
+  #endif
+
+  // Return
+  return;
 }
  
  
@@ -1280,55 +1382,105 @@ void GSparseMatrix::alloc_elements(int start, int num)
  ***************************************************************************/
 void GSparseMatrix::free_elements(int start, int num)
 {
+  // Dump header
+  #if defined(G_DEBUG_SPARSE_MALLOC)
+  cout << "GSparseMatrix::free_elements(start=" << start << ", num=" << 
+          num << ")" << endl;
+  #endif
+
   // Continue only if we need to free memory and if start is within the
   // range
   if (num > 0 && start < m_elements) {
  
-    // Determine the new size of the matrix. If there are no elements then
-	// simply delete all matrix elements ...
+    // Determine the requested new logical size of the matrix
 	int new_size = m_elements - num;
+
+    // If there are no elements then simply delete all matrix elements ...
 	if (new_size < 1) {
       if (m_data   != NULL) delete [] m_data;
       if (m_rowinx != NULL) delete [] m_rowinx;
 	  m_data     = NULL;
 	  m_rowinx   = NULL;
 	  m_elements = 0;
+	  m_alloc    = 0;
 	}
  
 	// ... otherwise shrink the array
 	else {
  
-	  // Allocate new memory
-	  double* new_data   = new double[new_size];
-	  int*    new_rowinx = new int[new_size];
-	  if (new_data == NULL || new_rowinx == NULL)
-	    throw mem_alloc("GSparseMatrix::free_elements(int, int)", 
-	                    new_size);
+	  // Case A: If at least one entire memory block has been liberated then 
+	  // physically shrink the matrix array
+	  if (m_alloc - new_size > m_mem_block) {
+	  
+	    // Shrink, but leave a memory block for possible future filling
+	    m_alloc = new_size + m_mem_block;
+	  
+	    // Allocate new memory
+	    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::free_elements(int, int)", m_alloc);
+
+	    // Copy all elements before the starting index
+	    #if defined(G_USE_MEMCPY)
+	    int n_copy = start;
+	    memcpy(new_data,   m_data,   sizeof(double)*n_copy);
+	    memcpy(new_rowinx, m_rowinx, sizeof(long)*n_copy);
+		#else
+	    for (int i = 0; i < start; ++i) {
+	      new_data[i]   = m_data[i];
+	      new_rowinx[i] = m_rowinx[i];
+	    }
+		#endif
+
+	    // Copy all elements after the starting index
+	    #if defined(G_USE_MEMCPY)
+	    n_copy = new_size - start;
+	    memcpy(&(new_data[start]),   &(m_data[start+num]),   sizeof(double)*n_copy);
+	    memcpy(&(new_rowinx[start]), &(m_rowinx[start+num]), sizeof(long)*n_copy);
+	    #else
+	    for (int i = start; i < new_size; ++i) {
+	      new_data[i]   = m_data[i+num];
+	      new_rowinx[i] = m_rowinx[i+num];
+	    }
+		#endif
+		
+        // Delete old 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 = new_size;
+		
+	  } // endif: Case A: memory shrinkage performed
+	  
+	  // Case B: we keep the memory and just move the elements
+	  else {
+
+	    // Move all elements after the starting index
+	    #if defined(G_USE_MEMCPY)
+	    int n_copy = new_size - start;
+	    memmove(&(m_data[start]),   &(m_data[start+num]),   sizeof(double)*n_copy);
+	    memmove(&(m_rowinx[start]), &(m_rowinx[start+num]), sizeof(long)*n_copy);
+		#else
+	    for (int i = start; i < new_size; ++i) {
+	      m_data[i]   = m_data[i+num];
+	      m_rowinx[i] = m_rowinx[i+num];
+	    }
+		#endif
  
-	  // Copy all elements before the starting index
-	  for (int i = 0; i < start; ++i) {
-	    new_data[i]   = m_data[i];
-	    new_rowinx[i] = m_rowinx[i];
-	  }
+	    // Update element counter
+	    m_elements = new_size;
  
-	  // Copy all elements after the starting index
-	  for (int i = start; i < new_size; ++i) {
-	    new_data[i]   = m_data[i+num];
-	    new_rowinx[i] = m_rowinx[i+num];
-	  }
-	
-      // Delete old 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 = new_size;
+	  } // endelse: Case B
  
 	} // endif: array shrinkage needed	
   } // endif: needed new memory
  
+  // Return
+  return;
 }
  
  
@@ -1374,8 +1526,8 @@ void GSparseMatrix::remove_zero_row_col(void)
     row_map[m_rowsel[c_row]] = c_row;
  
   // Initialise pointers to compressed array
-  double* d_data     = m_data;
-  int*    d_rowinx   = m_rowinx;
+  double* d_data   = m_data;
+  int*    d_rowinx = m_rowinx;
  
   // Initialise column start of first column to zero
   m_colstart[0] = 0;
@@ -1581,32 +1733,41 @@ ostream&amp; operator&lt;&lt; (ostream&amp; os, const GSparseMatrix&amp; m)
   else
     os << " Number of columns .........: " << m.m_cols << endl;
   os << " Number of non-zero elements: " << m.m_elements << endl;
+  os << " Number of allocated cells .: " << m.m_alloc << endl;
+  os << " Memory block size .........: " << m.m_mem_block << endl;
   os << " Sparse matrix fill ........: " << m.fill() << endl;
  
-  // Loop over all matrix elements and put them into the stream
-  for (int row = 0; row < m.m_rows; ++row) {
-    os << " ";
-    for (int col = 0; col < m.m_cols; ++col) {
-      os << m(row,col);
-	  if (col != m.m_cols-1)
-	    os << ", ";
-	}
-	if (row != m.m_rows-1)
-	  os << endl;
+  // Loop over all matrix elements and put them into the stream (only if
+  // their number is not too large)
+  if (m.m_rows < 20 && m.m_cols < 20) {
+    for (int row = 0; row < m.m_rows; ++row) {
+      os << " ";
+      for (int col = 0; col < m.m_cols; ++col) {
+        os << m(row,col);
+	    if (col != m.m_cols-1)
+	      os << ", ";
+	  }
+	  if (row != m.m_rows-1)
+	    os << endl;
+    }
   }
  
   // If there is a row compression the show scheme
   if (m.m_rowsel != NULL) {
     os << endl << " Row selection ..:";
-    for (int row = 0; row < m.m_num_rowsel; ++row)
-      os << " " << m.m_rowsel[row];
+    if (m.m_num_rowsel < 20) {
+      for (int row = 0; row < m.m_num_rowsel; ++row)
+        os << " " << m.m_rowsel[row];
+	}
   }
  
   // If there is a column compression the show scheme
   if (m.m_colsel != NULL) {
     os << endl << " Column selection:";
-    for (int col = 0; col < m.m_num_colsel; ++col)
-      os << " " << m.m_colsel[col];
+    if (m.m_num_colsel < 20) {
+      for (int col = 0; col < m.m_num_colsel; ++col)
+        os << " " << m.m_colsel[col];
+	}
   }
  
   // If there is a symbolic decomposition then put it also in the stream
@@ -22,6 +22,9 @@
 /* __ Namespaces _________________________________________________________ */
 using namespace std;
  
+/* __ Definitions ________________________________________________________ */
+#define G_SPARSE_MATRIX_DEFAULT_MEM_BLOCK 512    // Default Memory block size
+
 /* __ Enumerators ________________________________________________________ */
  
 /* __ Structures _________________________________________________________ */
@@ -94,6 +97,7 @@ public:
   virtual GVector extract_col(int col) const;
   virtual void    insert_col(const GVector& v, int col);
   virtual void    add_col(const GVector& v, int col);
+  virtual void    set_mem_block(int block);
   virtual GMatrix convert_to_full() const;
   //TBD virtual GMatrix extract_lower_triangle() const;
   //TBD virtual GMatrix extract_upper_triangle() const;
@@ -112,13 +116,14 @@ private:
   GSparseMatrix cs_transpose(GSparseMatrix *m, int values);
  
   // Private data members
-  int*   m_rowinx;    // Row-indices of all elements
-  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
 };
  
  
@@ -245,14 +250,6 @@ GSparseMatrix operator/ (const GSparseMatrix&amp; a, const double&amp; b)
   return result;
 }
  
-// Matrix clear function
-inline
-void GSparseMatrix::clear()
-{
-  free_elements(0, m_elements);
-  return;
-}
-
 // Matrix fill function
 inline
 double GSparseMatrix::fill() const
@@ -269,6 +266,14 @@ GSparseMatrix transpose(const GSparseMatrix&amp; m)
   return result;
 }
  
+// Set memory block size
+inline
+void GSparseMatrix::set_mem_block(int block)
+{
+  m_mem_block = (block > 0) ? block : 1;
+  return;
+}
+
 // Cholesky decomposition
 inline
 GSparseMatrix cholesky_decompose(const GSparseMatrix& m, int compress)
@@ -33,6 +33,7 @@ GSymMatrix::GSymMatrix(int rows, int cols) : GMatrix()
   m_rows       = 0;
   m_cols       = 0;
   m_elements   = 0;
+  m_alloc      = 0;
   m_num_rowsel = 0;
   m_num_colsel = 0;
   m_data       = NULL;
@@ -63,6 +64,7 @@ GSymMatrix::GSymMatrix(int rows, int cols) : GMatrix()
   m_rows     = rows;
   m_cols     = cols;
   m_elements = elements;
+  m_alloc    = elements;
  
   // Set-up column start indices
   m_colstart[0]   = 0;
@@ -73,6 +73,25 @@ GVector::~GVector()
  =                                                                         =
  ==========================================================================*/
  
+/***************************************************************************
+ *                   Returns number of non-zeros in vector                 *
+ ***************************************************************************/
+int GVector::non_zeros() const
+{
+  // Initialise number of non-zeros
+  int non_zeros = 0;
+  
+  // Gather all non-zero elements
+  for (int i = 0; i < m_num; ++i) {
+    if (m_data[i] != 0.0)
+	  non_zeros++;
+  }
+  
+  // Return number of non-zero elements
+  return non_zeros;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                         GVector private functions                       =
@@ -104,7 +104,8 @@ public:
   GVector  operator- () const;
  
   // Vector functions
-  int size() const; // Return dimension of vector
+  int size() const;       // Return dimension of vector
+  int non_zeros() const;  // Returns # of non zeros in vector
  
   // Exception: Vector index out of range
   class out_of_range : public GException {
 ## 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
  
-lib_LTLIBRARIES = libgamma.la
-libgamma_la_SOURCES = $(files)
 ## 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_GMatrix_SOURCES = test_GMatrix.cpp
-test_GMatrix_LDADD = ../src/libgamma.la
-test_GSymMatrix_SOURCES = test_GSymMatrix.cpp
-test_GSymMatrix_LDADD = ../src/libgamma.la
-test_GSparseMatrix_SOURCES = test_GSparseMatrix.cpp
-test_GSparseMatrix_LDADD = ../src/libgamma.la
+SUBDIRS = src test
+
@@ -265,6 +265,23 @@ int check_matrix_sum(const double sum)
  
  
 /***************************************************************************
+ *                                Test: Output                             *
+ ***************************************************************************/
+void test_output(const GMatrix& m_test)
+{
+  cout << "Test GMatrix: Output test matrix: " << endl;
+  try {
+    cout << m_test << endl;
+  }
+  catch (exception &e) {
+	cout << endl << "TEST ERROR: Unable to output matrix." << endl;
+    cout << e.what() << endl;
+	throw;
+  }
+}
+
+
+/***************************************************************************
  *                            Main test function                           *
  ***************************************************************************/
 int main(void)
@@ -285,6 +302,9 @@ int main(void)
   // Prepare result matrix
   GMatrix result = m_test;
  
+  // Execute the tests
+  test_output(m_test);
+
   // Test 1: Allocate zero matrix
   try {
     cout << "GMatrix - Test 1: Allocate zero matrix: ";
@@ -1340,7 +1340,7 @@ void test_cholesky(void)
 /***************************************************************************
  *                            Test: Heavy calculus                         *
  ***************************************************************************/
-void test_heavy()
+void test_heavy(int block = 512)
 {
   // Allocate some variables
   clock_t t_start;
@@ -1353,14 +1353,15 @@ void test_heavy()
     //
 	// Fill 10000 x 10000 matrix with 10000 values
 	number  = 10000;
-    t_start = clock();
 	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);
     #if defined(DUMP_TIMING)
 	cout << endl << " - Fill of 10000 values needed " << t_elapse << 
-	                " sec (reference ~ 1.9 sec)";
+	                " sec (reference ~ 0.2 sec)";
 	#endif   
     //
 	// Modify 10000 x 10000 matrix with 10000 (exisiting) values
@@ -1373,12 +1374,33 @@ void test_heavy()
 	                " sec (reference ~ 0 sec)";
 	#endif   
     //
+	// Insert columns into 10000 x 10000 matrix
+	number  = 10000;
+	large.clear();
+	large.set_mem_block(block);
+	GVector column(number);
+    t_start = clock();
+	for (int j = 0; j < 100; ++j) {
+	  column = 0.0;
+	  for (int i = j; i < number; ++i)
+	    column(i) = (i+1)*3.14;
+	  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)";
+	#endif   
+    //
 	// 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;
@@ -1392,8 +1414,9 @@ void test_heavy()
 	#endif   
     //
 	// Subsequent multiplications
-    t_start = clock();
 	GSparseMatrix fac(number,number);
+	fac.set_mem_block(block);
+    t_start = clock();
 	for (int i = 0; i < 1; ++i) {
 	  for (int j = 0; j < number; ++j)
 	    fac(j,j) = (i+1)*1.0;
@@ -290,6 +290,23 @@ int check_matrix_sum(const double sum)
  
  
 /***************************************************************************
+ *                                Test: Output                             *
+ ***************************************************************************/
+void test_output(const GSymMatrix& m_test)
+{
+  cout << "Test GSymMatrix: Output test matrix: " << endl;
+  try {
+    cout << m_test << endl;
+  }
+  catch (exception &e) {
+	cout << endl << "TEST ERROR: Unable to output matrix." << endl;
+    cout << e.what() << endl;
+	throw;
+  }
+}
+
+
+/***************************************************************************
  *                            Main test function                           *
  ***************************************************************************/
 int main(void)
@@ -310,6 +327,9 @@ int main(void)
   // Prepare result matrix
   GSymMatrix result = m_test;
  
+  // Execute the tests
+  test_output(m_test);
+
   // Test 1: Allocate zero matrix
   try {
     cout << "GSymMatrix - Test 1: Allocate zero matrix: ";