diff --git a/src/GMatrix.cpp b/src/GMatrix.cpp index af86b75..38aaf5c 100644 --- a/src/GMatrix.cpp +++ b/src/GMatrix.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& 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& 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& 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& GMatrix::operator= (const GMatrix& 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& GMatrix::operator= (const GMatrix& 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& GMatrix::operator= (const GMatrix& 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& operator<< (ostream& os, const GMatrix& 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; diff --git a/src/GMatrix.hpp b/src/GMatrix.hpp index e6c31c8..e1adbbf 100644 --- a/src/GMatrix.hpp +++ b/src/GMatrix.hpp @@ -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) diff --git a/src/GSparseMatrix.cpp b/src/GSparseMatrix.cpp index 52df670..4331711 100644 --- a/src/GSparseMatrix.cpp +++ b/src/GSparseMatrix.cpp @@ -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& 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& GSparseMatrix::operator= (const GSparseMatrix& 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& GSparseMatrix::operator*= (const GSparseMatrix& 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& 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& 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& 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& operator<< (ostream& os, const GSparseMatrix& 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 diff --git a/src/GSparseMatrix.hpp b/src/GSparseMatrix.hpp index b5376d4..3f031df 100644 --- a/src/GSparseMatrix.hpp +++ b/src/GSparseMatrix.hpp @@ -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& a, const double& 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& 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) diff --git a/src/GSymMatrix.cpp b/src/GSymMatrix.cpp index d3e6727..6807f90 100644 --- a/src/GSymMatrix.cpp +++ b/src/GSymMatrix.cpp @@ -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; diff --git a/src/GVector.cpp b/src/GVector.cpp index 5d574b3..0078723 100644 --- a/src/GVector.cpp +++ b/src/GVector.cpp @@ -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 = diff --git a/src/GVector.hpp b/src/GVector.hpp index c95bc90..5aa0b2c 100644 --- a/src/GVector.hpp +++ b/src/GVector.hpp @@ -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 { diff --git a/src/Makefile.am b/src/Makefile.am index 20481c0..1ae3f5f 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,12 +1,4 @@ ## 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) diff --git a/test/Makefile.am b/test/Makefile.am index a8b0a4c..1ae3f5f 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -1,18 +1,4 @@ ## 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 + diff --git a/test/test_GMatrix.cpp b/test/test_GMatrix.cpp index 9ba71b9..26b11f1 100644 --- a/test/test_GMatrix.cpp +++ b/test/test_GMatrix.cpp @@ -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: "; diff --git a/test/test_GSparseMatrix.cpp b/test/test_GSparseMatrix.cpp index 05b0879..bf6027a 100644 --- a/test/test_GSparseMatrix.cpp +++ b/test/test_GSparseMatrix.cpp @@ -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; diff --git a/test/test_GSymMatrix.cpp b/test/test_GSymMatrix.cpp index 507391e..ed361a5 100644 --- a/test/test_GSymMatrix.cpp +++ b/test/test_GSymMatrix.cpp @@ -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: ";