diff --git a/src/GMatrix.cpp b/src/GMatrix.cpp index 38aaf5c..7b5e962 100644 --- a/src/GMatrix.cpp +++ b/src/GMatrix.cpp @@ -170,9 +170,9 @@ GMatrix::GMatrix(const GMatrix& m) GMatrix::~GMatrix() { // De-allocate only if memory 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_colstart != NULL) delete [] m_colstart; if (m_data != NULL) delete [] m_data; // Return diff --git a/src/GSparseMatrix.cpp b/src/GSparseMatrix.cpp index 4331711..953f81a 100644 --- a/src/GSparseMatrix.cpp +++ b/src/GSparseMatrix.cpp @@ -48,8 +48,8 @@ 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_alloc = 0; + m_elements = 0; // Logically used # of elements + m_alloc = 0; // Allocated # of elements (>= m_elements) m_data = NULL; m_colstart = NULL; m_rowinx = NULL; @@ -537,7 +537,14 @@ double GSparseMatrix::sum() const ***************************************************************************/ void GSparseMatrix::transpose() { - *this = cs_transpose(this, 1); + // Fill pending value in input matrix + fill_pending(); + + // Transpose matrix + *this = cs_transpose(*this, 1); + + // Return + return; } @@ -803,8 +810,12 @@ void GSparseMatrix::cholesky_decompose(int compress) int matrix_rows = m_rows; int matrix_cols = m_cols; - // Delete any existing symbolic analysis object + // Delete any existing symbolic and numeric analysis object and reset + // pointers if (m_symbolic != NULL) delete (GSparseSymbolic*)m_symbolic; + if (m_numeric != NULL) delete (GSparseNumeric*)m_numeric; + m_symbolic = NULL; + m_numeric = NULL; // Allocate symbolic analysis object GSparseSymbolic* symbolic = new GSparseSymbolic(); @@ -827,44 +838,26 @@ void GSparseMatrix::cholesky_decompose(int compress) // Ordering an symbolic analysis of matrix. This sets up an array 'pinv' // which contains the fill-in reducing permutations symbolic->cholesky_symbolic_analysis(1, *this); + + // Store symbolic pointer in sparse matrix object + m_symbolic = (void*)symbolic; // Perform numeric Cholesky decomposition numeric.cholesky_numeric_analysis(*this, *symbolic); - - // Insert zero rows and columns if they have been removed previously. - // Since no selection information is present in numeric.m_L we have to - // copy it from this in order to allow for insertion. - if (compress) { - - // If there is a row selection then copy it - if (m_rowsel != NULL && m_num_rowsel > 0) { - numeric.m_L->m_rowsel = new int[m_num_rowsel]; - if (numeric.m_L->m_rowsel == NULL) - throw mem_alloc("GSparseMatrix::cholesky_decompose(int)", m_num_rowsel); - for (int i = 0; i < m_num_rowsel; ++i) - numeric.m_L->m_rowsel[i] = m_rowsel[i]; - numeric.m_L->m_num_rowsel = m_num_rowsel; - } - - // If there is a column selection then copy it - if (m_colsel != NULL && m_num_colsel > 0) { - numeric.m_L->m_colsel = new int[m_num_colsel]; - if (numeric.m_L->m_colsel == NULL) - throw mem_alloc("GSparseMatrix::cholesky_decompose(int)", m_num_colsel); - for (int i = 0; i < m_num_colsel; ++i) - numeric.m_L->m_colsel[i] = m_colsel[i]; - numeric.m_L->m_num_colsel = m_num_colsel; - } - - // Now we can insert zero rows and columns - numeric.m_L->insert_zero_row_col(matrix_rows, matrix_cols); + + // Copy L matrix into this object + free_elements(0, m_elements); + alloc_elements(0, numeric.m_L->m_elements); + for (int i = 0; i < m_elements; ++i) { + m_data[i] = numeric.m_L->m_data[i]; + m_rowinx[i] = numeric.m_L->m_rowinx[i]; } + for (int col = 0; col <= m_cols; ++col) + m_colstart[col] = numeric.m_L->m_colstart[col]; - // Store Cholesky decomposition in matrix - *this = *(numeric.m_L); - - // Store symbolic pointer in sparse matrix object - m_symbolic = (void*)symbolic; + // Insert zero rows and columns if they have been removed previously. + if (compress) + insert_zero_row_col(matrix_rows, matrix_cols); // Return return; @@ -1534,7 +1527,7 @@ void GSparseMatrix::remove_zero_row_col(void) // Dump column start array #if defined(G_DEBUG_SPARSE_COMPRESSION) - cout << " before compression:"; + cout << " before compression (" << m_colstart[m_cols] << "):"; for (int col = 0; col <= m_cols; ++col) cout << " " << m_colstart[col]; cout << endl; @@ -1567,7 +1560,7 @@ void GSparseMatrix::remove_zero_row_col(void) // Dump column start array #if defined(G_DEBUG_SPARSE_COMPRESSION) - cout << " after compression :"; + cout << " after compression (" << m_colstart[m_num_colsel] << ") :"; for (int c_col = 0; c_col <= m_num_colsel; ++c_col) cout << " " << m_colstart[c_col]; cout << endl; @@ -1576,10 +1569,9 @@ void GSparseMatrix::remove_zero_row_col(void) // Free row mapping array delete [] row_map; - // Update matrix attributes - m_rows = m_num_rowsel; - m_cols = m_num_colsel; - m_elements = m_colstart[m_num_colsel]; + // Update matrix size + m_rows = m_num_rowsel; + m_cols = m_num_colsel; // Return return; @@ -1616,7 +1608,7 @@ void GSparseMatrix::insert_zero_row_col(int rows, int cols) // Dump column start array #if defined(G_DEBUG_SPARSE_COMPRESSION) - cout << " before restoration:"; + cout << " before restoration (" << m_colstart[m_num_colsel] << "):"; for (int c_col = 0; c_col <= m_num_colsel; ++c_col) cout << " " << m_colstart[c_col]; cout << endl; @@ -1624,7 +1616,11 @@ void GSparseMatrix::insert_zero_row_col(int rows, int cols) // If column selection exists then restore column counters if (m_colsel != NULL) { + + // Start insertion from the last original column int col_stop = cols - 1; + + // Loop over all compressed columns for (int c_col = m_num_colsel-1; c_col > 0; --c_col) { int col_start = m_colsel[c_col-1] + 1; int col_value = m_colstart[c_col]; @@ -1632,18 +1628,24 @@ void GSparseMatrix::insert_zero_row_col(int rows, int cols) m_colstart[col] = col_value; col_stop = col_start - 1; } + + // Set first columns if they are not yet set + for (int col = 0; col <= col_stop; ++col) + m_colstart[col] = 0; + + // Restore the number of elements m_colstart[cols] = m_elements; } // Dump column start array #if defined(G_DEBUG_SPARSE_COMPRESSION) - cout << " after restoration :"; + cout << " after restoration (" << m_colstart[cols] << ") :"; for (int col = 0; col <= cols; ++col) cout << " " << m_colstart[col]; cout << endl; #endif - // Update matrix attributes + // Update matrix size m_rows = rows; m_cols = cols; @@ -1652,63 +1654,6 @@ void GSparseMatrix::insert_zero_row_col(int rows, int cols) } -/*************************************************************************** - * cs_transpose * - * ----------------------------------------------------------------------- * - * Transpose matrix. The flag 'values' allows to avoid copying the actual * - * data values. This allows to perform a logical matrix transposition, as * - * needed by the symbolic matrix analysis class. * - ***************************************************************************/ -GSparseMatrix GSparseMatrix::cs_transpose(GSparseMatrix* m, int values) -{ - // Fill pending value in input matrix - m->fill_pending(); - - // Declare and allocate result matrix - GSparseMatrix result(m->m_cols, m->m_rows, m->m_elements); - - // Allocate and initialise workspace - int wrk_size = m->m_rows; - int* wrk_int = new int[wrk_size]; - if (wrk_int == NULL) - throw mem_alloc("GSparseMatrix::cs_transpose(GSparseMatrix*, int)", - wrk_size); - for (int i = 0; i < wrk_size; i++) - wrk_int[i] = 0; - - // Setup the number of non-zero elements in each row - for (int p = 0; p < m->m_elements; p++) - wrk_int[m->m_rowinx[p]]++; - - // Set row pointers. To use a GSparseSymbolic function we have to - // allocate and object (but this does not take memory) - cs_cumsum(result.m_colstart, wrk_int, m->m_rows); - - // Case A: Normal transponse, including assignment of values - if (values) { - for (int col = 0; col < m->m_cols; col++) { - for (int p = m->m_colstart[col] ; p < m->m_colstart[col+1] ; p++) { - int i = wrk_int[m->m_rowinx[p]]++; - result.m_rowinx[i] = col; - result.m_data[i] = m->m_data[p] ; - } - } - } - - // Case B: Logical transponse, no assignment of values is performed - else { - for (int col = 0; col < m->m_cols; col++) { - for (int p = m->m_colstart[col] ; p < m->m_colstart[col+1] ; p++) - result.m_rowinx[wrk_int[m->m_rowinx[p]]++] = col; - } - } - - // Return transponse matrix - return result; - -} - - /*========================================================================== = = = GSparseMatrix friends = @@ -1732,42 +1677,36 @@ ostream& operator<< (ostream& os, const GSparseMatrix& m) m.m_num_colsel << ")" << endl; else os << " Number of columns .........: " << m.m_cols << endl; - os << " Number of non-zero elements: " << m.m_elements << endl; + os << " Number of non-zero elements: " << m.m_colstart[m.m_cols] << " (" << + 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 (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; - } + // 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; } // If there is a row compression the show scheme if (m.m_rowsel != NULL) { os << endl << " Row selection ..:"; - if (m.m_num_rowsel < 20) { - for (int row = 0; row < m.m_num_rowsel; ++row) - os << " " << m.m_rowsel[row]; - } + 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:"; - if (m.m_num_colsel < 20) { - for (int col = 0; col < m.m_num_colsel; ++col) - os << " " << m.m_colsel[col]; - } + 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 @@ -1805,6 +1744,145 @@ GSparseMatrix fabs(const GSparseMatrix& m) } +/*************************************************************************** + * cs_symperm * + * ----------------------------------------------------------------------- * + * C = A(p,p) where A and C are symmetric the upper part stored. * + * ----------------------------------------------------------------------- * + * Input: A Sparse matrix * + * pinv pinv[0..n-1] * + * Output: C Sparse matrix * + ***************************************************************************/ +GSparseMatrix cs_symperm(const GSparseMatrix& m, const int* pinv) +{ + // Declare loop variables + int i, j, p, q, i2, j2; + + // Assign matrix attributes + int n = m.m_cols; + int* Ap = m.m_colstart; + int* Ai = m.m_rowinx; + double* Ax = m.m_data; + + // Allocate result matrix + GSparseMatrix C(n, n, Ap[n]); + + // Allocate and initialise workspace + int wrk_size = n; + int* wrk_int = new int[wrk_size]; + if (wrk_int == NULL) + throw mem_alloc("cs_symperm(GSparseMatrix*, const int*, int", wrk_size); + for (i = 0; i < wrk_size; ++i) + wrk_int[i] = 0; + + // Assign result matrix attributes + int* Cp = C.m_colstart; + int* Ci = C.m_rowinx; + double* Cx = C.m_data; + + // Count entries in each column of C + for (j = 0; j < n; j++) { + + // Column j of A is column j2 of C + j2 = pinv ? pinv[j] : j; + + // Loop over entries in column j + for (p = Ap[j]; p < Ap[j+1]; p++) { + i = Ai [p]; + if (i > j) continue; // skip lower triangular part of A + i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C + wrk_int[G_MAX(i2, j2)]++; // column count of C + } + } + + // Compute column pointers of C + cs_cumsum(Cp, wrk_int, n); + + // Loop over all columns of A + for (j = 0 ; j < n ; j++) { + + // Column j of A is column j2 of C + j2 = pinv ? pinv[j] : j; + + // Loop over entries in column j + for (p = Ap[j]; p < Ap[j+1]; p++) { + i = Ai [p] ; + if (i > j) continue; // skip lower triangular part of A + i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C + Ci[q = wrk_int[G_MAX(i2,j2)]++] = G_MIN(i2,j2); + if (Cx) Cx[q] = Ax[p]; + } + } + + // Free workspace + delete [] wrk_int; + + // Rectify the number of elements in matrix C + C.free_elements(Cp[n], (C.m_elements-Cp[n])); + + // Return result + return C; + +} + + +/*************************************************************************** + * cs_transpose * + * ----------------------------------------------------------------------- * + * Transpose matrix. The flag 'values' allows to avoid copying the actual * + * data values. This allows to perform a logical matrix transposition, as * + * needed by the symbolic matrix analysis class. * + * ----------------------------------------------------------------------- * + * NOTE: This routine DOES NOT support pending elements (they have to be * + * filled before. * + ***************************************************************************/ +GSparseMatrix cs_transpose(const GSparseMatrix& m, int values) +{ + // Declare and allocate result matrix + GSparseMatrix result(m.m_cols, m.m_rows, m.m_elements); + + // Allocate and initialise workspace + int wrk_size = m.m_rows; + int* wrk_int = new int[wrk_size]; + if (wrk_int == NULL) + throw mem_alloc("GSparseMatrix::cs_transpose(GSparseMatrix*, int)", + wrk_size); + for (int i = 0; i < wrk_size; i++) + wrk_int[i] = 0; + + // Setup the number of non-zero elements in each row + for (int p = 0; p < m.m_elements; p++) + wrk_int[m.m_rowinx[p]]++; + + // Set row pointers. To use a GSparseSymbolic function we have to + // allocate and object (but this does not take memory) + cs_cumsum(result.m_colstart, wrk_int, m.m_rows); + + // Case A: Normal transponse, including assignment of values + if (values) { + for (int col = 0; col < m.m_cols; col++) { + for (int p = m.m_colstart[col] ; p < m.m_colstart[col+1] ; p++) { + int i = wrk_int[m.m_rowinx[p]]++; + result.m_rowinx[i] = col; + result.m_data[i] = m.m_data[p] ; + } + } + } + + // Case B: Logical transponse, no assignment of values is performed + else { + for (int col = 0; col < m.m_cols; col++) { + for (int p = m.m_colstart[col] ; p < m.m_colstart[col+1] ; p++) + result.m_rowinx[wrk_int[m.m_rowinx[p]]++] = col; + } + } + + // Return transponse matrix + return result; + +} + + /*========================================================================== = = = GSparseMatrix exception classes = @@ -1849,4 +1927,3 @@ double cs_cumsum(int* p, int* c, int n) // Return cumulative sum of c[0..n-1] return nz2; } - diff --git a/src/GSparseMatrix.hpp b/src/GSparseMatrix.hpp index 3f031df..e702f04 100644 --- a/src/GSparseMatrix.hpp +++ b/src/GSparseMatrix.hpp @@ -25,12 +25,9 @@ using namespace std; /* __ Definitions ________________________________________________________ */ #define G_SPARSE_MATRIX_DEFAULT_MEM_BLOCK 512 // Default Memory block size -/* __ Enumerators ________________________________________________________ */ - -/* __ Structures _________________________________________________________ */ - -/* __ Prototypes _________________________________________________________ */ -double cs_cumsum(int* p, int* c, int n); +/* __ Macros (to be moved to GammaLib general header) ____________________ */ +#define G_MIN(a,b) (((a) < (b)) ? (a) : (b)) +#define G_MAX(a,b) (((a) > (b)) ? (a) : (b)) /*************************************************************************** @@ -54,6 +51,9 @@ class GSparseMatrix : public GMatrix { friend GSparseMatrix fabs(const GSparseMatrix& m); friend GSparseMatrix cholesky_decompose(const GSparseMatrix& m, int compress = 1); friend GSparseMatrix cholesky_invert(const GSparseMatrix& m, int compress = 1); + // + friend GSparseMatrix cs_symperm(const GSparseMatrix& m, const int* pinv); + friend GSparseMatrix cs_transpose(const GSparseMatrix& m, int values); public: // Constructors and destructors (not inherited) @@ -113,7 +113,6 @@ private: void free_elements(int start, int num); void remove_zero_row_col(void); void insert_zero_row_col(int rows, int cols); - GSparseMatrix cs_transpose(GSparseMatrix *m, int values); // Private data members int* m_rowinx; // Row-indices of all elements @@ -292,4 +291,10 @@ GSparseMatrix cholesky_invert(const GSparseMatrix& m, int compress) return result; } + +/*************************************************************************** + * Prototypes for other functions * + ***************************************************************************/ +double cs_cumsum(int* p, int* c, int n); + #endif /* GSPARSEMATRIX_HPP */ diff --git a/src/GSparseNumeric.cpp b/src/GSparseNumeric.cpp index 6710dd9..5eb0905 100644 --- a/src/GSparseNumeric.cpp +++ b/src/GSparseNumeric.cpp @@ -78,6 +78,12 @@ GSparseNumeric& GSparseNumeric::operator= (const GSparseNumeric& n) // Execute only if object is not identical if (this != &n) { + // De-allocate only if memory has indeed been allocated + if (m_L != NULL) delete m_L; + if (m_U != NULL) delete m_U; + if (m_pinv != NULL) delete [] m_pinv; + if (m_B != NULL) delete [] m_B; + // Initialise private members for clean destruction m_L = NULL; m_U = NULL; @@ -87,15 +93,23 @@ GSparseNumeric& GSparseNumeric::operator= (const GSparseNumeric& n) m_n_B = 0; // Copy m_L if it exists - if (n.m_L != NULL) - m_L = new GSparseMatrix(*n.m_L); + if (n.m_L != NULL) { + m_L = new GSparseMatrix(*n.m_L); + if (m_L == NULL) + throw mem_alloc("GSparseNumeric::operator= (const GSparseNumeric&)", + 1); + } // Copy m_U if it exists - if (n.m_U != NULL) - m_U = new GSparseMatrix(*n.m_U); + if (n.m_U != NULL) { + m_U = new GSparseMatrix(*n.m_U); + if (m_U == NULL) + throw mem_alloc("GSparseNumeric::operator= (const GSparseNumeric&)", + 1); + } // Copy m_pinv array if it exists - if (n.m_pinv != NULL) { + if (n.m_pinv != NULL && n.m_n_pinv > 0) { m_pinv = new int[n.m_n_pinv]; if (m_pinv == NULL) throw mem_alloc("GSparseNumeric::operator= (const GSparseNumeric&)", @@ -106,7 +120,7 @@ GSparseNumeric& GSparseNumeric::operator= (const GSparseNumeric& n) } // Copy m_B array if it exists - if (n.m_B != NULL) { + if (n.m_B != NULL && n.m_n_B > 0) { m_B = new double[n.m_n_B]; if (m_B == NULL) throw mem_alloc("GSparseNumeric::operator= (const GSparseNumeric&)", @@ -138,9 +152,23 @@ GSparseNumeric& GSparseNumeric::operator= (const GSparseNumeric& n) * Input: A Sparse matrix * * S Symbolic analysis of sparse matrix * ***************************************************************************/ -void GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix& A, +void GSparseNumeric::cholesky_numeric_analysis(const GSparseMatrix& A, const GSparseSymbolic& S) { + // De-allocate memory that has indeed been previously allocated + if (m_L != NULL) delete m_L; + if (m_U != NULL) delete m_U; + if (m_pinv != NULL) delete [] m_pinv; + if (m_B != NULL) delete [] m_B; + + // Initialise members + m_L = NULL; + m_U = NULL; + m_pinv = NULL; + m_B = NULL; + m_n_pinv = 0; + m_n_B = 0; + // Return if arrays in the symbolic analysis have not been allocated if (!S.m_cp || !S.m_parent) return; @@ -173,7 +201,7 @@ void GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix& A, int* parent = S.m_parent; // Assign C = A(p,p) where A and C are symmetric and the upper part stored - GSparseMatrix C = (pinv) ? S.cs_symperm(&A, pinv, 1) : (A); + GSparseMatrix C = (pinv) ? cs_symperm(A, pinv) : (A); // Assign workspace pointer int* c = wrk_int; @@ -186,7 +214,7 @@ void GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix& A, double* Cx = C.m_data; // Allocate L matrix - m_L = new GSparseMatrix(n,n,cp[n]); + m_L = new GSparseMatrix(n, n, cp[n]); if (m_L == NULL) throw mem_alloc( "GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix&, const GSparseSymbolic&)", @@ -204,7 +232,8 @@ void GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix& A, // Compute L(:,k) for L*L' = C for (k = 0; k < n; k++) { - // Nonzero pattern of L(k,:) + // Nonzero pattern of L(k,:). + // Returns -1 if parent = NULL, s = NULL or c = NULL top = cs_ereach(&C, k, parent, s, c); // find pattern of L(k,:) x[k] = 0; // x (0:k) is now zero @@ -222,14 +251,14 @@ void GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix& A, // Triangular solve: Solve L(0:k-1,0:k-1) * x = C(:,k) for ( ; top < n; top++) { - i = s[top]; // s [top..n-1] is pattern of L(k,:) - lki = x[i]/Lx[Lp[i]]; // L(k,i) = x (i) / L(i,i) - x[i] = 0; // clear x for k+1st iteration + i = s[top]; // s [top..n-1] is pattern of L(k,:) + lki = x[i]/Lx[Lp[i]]; // L(k,i) = x (i) / L(i,i) + x[i] = 0; // clear x for k+1st iteration for (p = Lp[i]+1; p < c[i]; p++) x[Li[p]] -= Lx[p] * lki; - d -= lki * lki; // d = d - L(k,i)*L(k,i) - p = c [i]++ ; - Li[p] = k; // store L(k,i) in column i + d -= lki * lki; // d = d - L(k,i)*L(k,i) + p = c[i]++; + Li[p] = k; // store L(k,i) in column i Lx[p] = lki; } @@ -319,7 +348,7 @@ int GSparseNumeric::cs_ereach(const GSparseMatrix* A, int k, CS_MARK(w, k); // Return index top, where s[top..n-1] contains pattern of L(k,:) - return (top); + return top; } diff --git a/src/GSparseNumeric.hpp b/src/GSparseNumeric.hpp index 5b7be6a..11f13cf 100644 --- a/src/GSparseNumeric.hpp +++ b/src/GSparseNumeric.hpp @@ -51,7 +51,7 @@ public: GSparseNumeric& operator= (const GSparseNumeric& n); // Functions - void cholesky_numeric_analysis(GSparseMatrix& m, const GSparseSymbolic& s); + void cholesky_numeric_analysis(const GSparseMatrix& m, const GSparseSymbolic& s); private: // Functions diff --git a/src/GSparseSymbolic.cpp b/src/GSparseSymbolic.cpp index 79605c6..3aecf30 100644 --- a/src/GSparseSymbolic.cpp +++ b/src/GSparseSymbolic.cpp @@ -24,6 +24,7 @@ using namespace std; /* __ Definitions ________________________________________________________ */ //#define G_DEBUG_SPARSE_SYM_AMD_ORDERING // Debug AMD ordering +//#define G_DEBUG_SPARSE_CHOLESKY // Analyse Cholesky decomposition /*========================================================================== @@ -63,11 +64,11 @@ GSparseSymbolic::GSparseSymbolic() GSparseSymbolic::~GSparseSymbolic() { // De-allocate only if memory has indeed been allocated - if (m_pinv != NULL) delete[] m_pinv; - if (m_q != NULL) delete[] m_q; - if (m_parent != NULL) delete[] m_parent; - if (m_cp != NULL) delete[] m_cp; - if (m_leftmost != NULL) delete[] m_leftmost; + if (m_pinv != NULL) delete [] m_pinv; + if (m_q != NULL) delete [] m_q; + if (m_parent != NULL) delete [] m_parent; + if (m_cp != NULL) delete [] m_cp; + if (m_leftmost != NULL) delete [] m_leftmost; // Return return; @@ -88,6 +89,13 @@ GSparseSymbolic& GSparseSymbolic::operator= (const GSparseSymbolic& s) // Execute only if object is not identical if (this != &s) { + // De-allocate only if memory has indeed been allocated + if (m_pinv != NULL) delete [] m_pinv; + if (m_q != NULL) delete [] m_q; + if (m_parent != NULL) delete [] m_parent; + if (m_cp != NULL) delete [] m_cp; + if (m_leftmost != NULL) delete [] m_leftmost; + // Initialise private members for clean destruction m_n_pinv = 0; m_n_q = 0; @@ -109,7 +117,7 @@ GSparseSymbolic& GSparseSymbolic::operator= (const GSparseSymbolic& s) m_unz = s.m_unz; // Copy m_pinv array if it exists - if (s.m_pinv != NULL) { + if (s.m_pinv != NULL && s.m_n_pinv > 0) { m_pinv = new int[s.m_n_pinv]; if (m_pinv == NULL) throw mem_alloc("GSparseSymbolic::operator= (const GSparseSymbolic&)", @@ -120,7 +128,7 @@ GSparseSymbolic& GSparseSymbolic::operator= (const GSparseSymbolic& s) } // Copy m_q array if it exists - if (s.m_q != NULL) { + if (s.m_q != NULL && s.m_n_q > 0) { m_q = new int[s.m_n_q]; if (m_q == NULL) throw mem_alloc("GSparseSymbolic::operator= (const GSparseSymbolic&)", @@ -131,7 +139,7 @@ GSparseSymbolic& GSparseSymbolic::operator= (const GSparseSymbolic& s) } // Copy m_parent array if it exists - if (s.m_parent != NULL) { + if (s.m_parent != NULL && s.m_n_parent > 0) { m_parent = new int[s.m_n_parent]; if (m_parent == NULL) throw mem_alloc("GSparseSymbolic::operator= (const GSparseSymbolic&)", @@ -142,7 +150,7 @@ GSparseSymbolic& GSparseSymbolic::operator= (const GSparseSymbolic& s) } // Copy m_cp array if it exists - if (s.m_cp != NULL) { + if (s.m_cp != NULL && s.m_n_cp > 0) { m_cp = new int[s.m_n_cp]; if (m_cp == NULL) throw mem_alloc("GSparseSymbolic::operator= (const GSparseSymbolic&)", @@ -153,20 +161,20 @@ GSparseSymbolic& GSparseSymbolic::operator= (const GSparseSymbolic& s) } // Copy m_leftmost array if it exists - if (s.m_leftmost != NULL) { + if (s.m_leftmost != NULL && s.m_n_leftmost > 0) { m_leftmost = new int[s.m_n_leftmost]; if (m_leftmost == NULL) throw mem_alloc("GSparseSymbolic::operator= (const GSparseSymbolic&)", s.m_n_leftmost); for (int i = 0; i < s.m_n_leftmost; ++i) m_leftmost[i] = s.m_leftmost[i]; - m_leftmost = s.m_leftmost; + m_n_leftmost = s.m_n_leftmost; } } // endif: object was not identical // Return - return (*this); + return *this; } @@ -186,44 +194,120 @@ GSparseSymbolic& GSparseSymbolic::operator= (const GSparseSymbolic& s) * 1: Cholesky decomposition * * m Sparse matrix * ***************************************************************************/ -void GSparseSymbolic::cholesky_symbolic_analysis(int order, GSparseMatrix& m) +void GSparseSymbolic::cholesky_symbolic_analysis(int order, + const GSparseMatrix& m) { + // Debug + #if defined(G_DEBUG_SPARSE_CHOLESKY) + cout << endl << "GSparseSymbolic::cholesky_symbolic_analysis entered" << endl; + cout << " order ......................: " << order << endl; + cout << " number of rows .............: " << m.m_rows << endl; + cout << " number of columns ..........: " << m.m_cols << endl; + cout << " number of non-zero elements : " << m.m_elements << endl; + #endif + + // De-allocate memory that has indeed been previously allocated + if (m_pinv != NULL) delete [] m_pinv; + if (m_q != NULL) delete [] m_q; + if (m_parent != NULL) delete [] m_parent; + if (m_cp != NULL) delete [] m_cp; + if (m_leftmost != NULL) delete [] m_leftmost; + + // Initialise members + m_n_pinv = 0; + m_n_q = 0; + m_n_parent = 0; + m_n_cp = 0; + m_n_leftmost = 0; + m_pinv = NULL; + m_q = NULL; + m_parent = NULL; + m_cp = NULL; + m_leftmost = NULL; + m_m2 = 0; + m_lnz = 0.0; + m_unz = 0.0; + // Check if order type is valid if (order < 0 || order > 1) - throw invalid_order("GSparseSymbolic::cholesky_symbolic_analysis(int, GSparseMatrix*)", - order, 0, 1); + throw invalid_order( + "GSparseSymbolic::cholesky_symbolic_analysis(int, GSparseMatrix*)", + order, 0, 1); - // Fill pending element - m.fill_pending(); - // Assign input matrix attributes int n = m.m_cols; // P = amd(A+A'), or natural int* P = cs_amd(order, &m); + #if defined(G_DEBUG_SPARSE_CHOLESKY) + cout << " AMD ordering permutation (" << P << ")" << endl; + if (P != NULL) { + for (int i = 0; i < n; ++i) + cout << " " << P[i]; + cout << endl; + } + #endif - // Find inverse permutation and store it in class member 'm_pinv' + // Find inverse permutation and store it in class member 'm_pinv'. + // Note that if P = NULL or n < 1 this function returns NULL. m_pinv = cs_pinv(P, n); m_n_pinv = n; + #if defined(G_DEBUG_SPARSE_CHOLESKY) + cout << " Inverse permutation (" << m_pinv << ")" << endl; + if (m_pinv != NULL) { + for (int i = 0; i < n; ++i) + cout << " " << m_pinv[i]; + cout << endl; + } + #endif // Delete workspace - delete [] P; + if (P != NULL) delete [] P; // C = spones(triu(A(P,P))) - GSparseMatrix C = cs_symperm(&m, m_pinv, 0); + GSparseMatrix C = cs_symperm(m, m_pinv); + #if defined(G_DEBUG_SPARSE_CHOLESKY) + cout << " C = spones(triu(A(P,P))) " << C << endl; + #endif // Find etree of C and store it in class member 'm_parent' m_parent = cs_etree(&C, 0); m_n_parent = n; + #if defined(G_DEBUG_SPARSE_CHOLESKY) + cout << " Elimination tree (" << m_parent << ")" << endl; + if (m_parent != NULL) { + for (int i = 0; i < n; ++i) + cout << " " << m_parent[i]; + cout << endl; + } + #endif // Post order the etree + // Note that if m_parent = NULL or n < 1 this function returns NULL. int* post = cs_post(m_parent, n); + #if defined(G_DEBUG_SPARSE_CHOLESKY) + cout << " Post ordered elimination tree (" << post << ")" << endl; + if (post != NULL) { + for (int i = 0; i < n; ++i) + cout << " " << post[i]; + cout << endl; + } + #endif // Find column counts of chol(C) + // Note that if m_parent = NULL or post = NULL this function returns NULL. int* c = cs_counts(&C, m_parent, post, 0); + #if defined(G_DEBUG_SPARSE_CHOLESKY) + cout << " Column counts (" << c << ")" << endl; + if (c != NULL) { + for (int i = 0; i < n; ++i) + cout << " " << c[i]; + cout << endl; + } + #endif // Delete workspace - delete [] post; + if (post != NULL) delete [] post; // Allocate column pointers for Cholesky decomposition m_n_cp = n+1; @@ -233,11 +317,52 @@ void GSparseSymbolic::cholesky_symbolic_analysis(int order, GSparseMatrix& m) // Find column pointers for L m_unz = m_lnz = cs_cumsum(m_cp, c, n); + #if defined(G_DEBUG_SPARSE_CHOLESKY) + cout << " Column pointers for L (" << m_cp << ")" << endl; + if (m_cp != NULL) { + for (int i = 0; i < m_n_cp; ++i) + cout << " " << m_cp[i]; + cout << endl; + } + cout << " Number of non-zero elements in L: " << m_lnz << endl; + #endif // Delete workspace - delete [] c; + if (c != NULL) delete [] c; + + // Delete symbolic analysis if it is not valid + if (m_lnz < 0) { + + // De-allocate memory that has indeed been previously allocated + if (m_pinv != NULL) delete [] m_pinv; + if (m_q != NULL) delete [] m_q; + if (m_parent != NULL) delete [] m_parent; + if (m_cp != NULL) delete [] m_cp; + if (m_leftmost != NULL) delete [] m_leftmost; -// return ((S->lnz >= 0) ? S : cs_sfree (S)) ; + // Initialise members + m_n_pinv = 0; + m_n_q = 0; + m_n_parent = 0; + m_n_cp = 0; + m_n_leftmost = 0; + m_pinv = NULL; + m_q = NULL; + m_parent = NULL; + m_cp = NULL; + m_leftmost = NULL; + m_m2 = 0; + m_lnz = 0.0; + m_unz = 0.0; + } + + // Debug + #if defined(G_DEBUG_SPARSE_CHOLESKY) + cout << "GSparseSymbolic::cholesky_symbolic_analysis finished" << endl; + #endif + + // Return + return; } @@ -261,11 +386,11 @@ void GSparseSymbolic::cholesky_symbolic_analysis(int order, GSparseMatrix& m) * A Sparse matrix * * Output: p[0..n] integer array of n+1 elements * ***************************************************************************/ -int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) +int* GSparseSymbolic::cs_amd(int order, const GSparseMatrix* A) { // Throw exception if order is not a decomposition if (order < 0 || order > 3) - throw invalid_order("GSparseSymbolic::cs_amd(int, GSparseMatrix*)", + throw invalid_order("GSparseSymbolic::cs_amd(int, const GSparseMatrix*)", order, 0, 3); // Debug @@ -305,17 +430,20 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) GSparseMatrix C(m,n); // Get (logical) transpose of A: AT = A' - GSparseMatrix AT = A->cs_transpose(A,0); +// NOTE: WE ONLY NEED THE LOGICAL TRANSPOSE HERE. HOWEVER, WE HAVE NO +// LOGICAL ADDITION SO FAR ... +// GSparseMatrix AT = cs_transpose(*A, 0); + GSparseMatrix AT = cs_transpose(*A, 1); // Find dense threshold dense = (int)CS_MAX(16, 10 * sqrt((double)n)); dense = CS_MIN(n-2, dense); // Case A: Cholesky decomposition of a symmetric matrix: C = A + A' - if (order == 1 && n == m) { + if (order == 1 && n == m) +// NOTE: HERE WE ONLY NEED LOGICAL ADDITION // C = cs_add (A, AT, 0, 0); C = *A + AT; - } // Case B: LU decomposition: C = A' * A with no dense rows else if (order == 2) { @@ -337,24 +465,35 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) ATp[m] = p2; // Get (logical) transpose of AT: A2 = AT' - GSparseMatrix A2 = AT.cs_transpose(&AT,0); +// NOTE: WE ONLY NEED THE LOGICAL TRANSPOSE HERE. HOWEVER, WE HAVE NO +// LOGICAL MULTIPLICATION SO FAR ... +// GSparseMatrix A2 = cs_transpose(AT, 0); + GSparseMatrix A2 = cs_transpose(AT, 1); // C = A' * A with no dense rows +// NOTE: cs_multiply NOT YET IMPLEMENTED // C = A2 ? cs_multiply(AT, A2) : NULL; C = AT * A2; } // Case C: QR decomposition: C = A' * A - else { + else +// NOTE: cs_multiply NOT YET IMPLEMENTED // C = cs_multiply(AT, A); C = AT * *A; - } - - // Stop in case of an error -// if (!C) return (NULL) ; + // Debug + #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING) + cout << " Matrix C " << C << endl; + #endif + // Drop diagonal entries from C - cs_fkeep(&C, &cs_diag, NULL); + int cnz_diag = cs_fkeep(&C, &cs_diag, NULL); + + // Debug + #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING) + cout << " Dropped diagonal entries from matrix C " << C << endl; + #endif // Assign C matrix attributes int* Cp = C.m_colstart; @@ -372,10 +511,19 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) throw mem_alloc("GSparseSymbolic::cs_amd (int, GSparseMatrix*)", wrk_size); // Add elbow room to C -// t = cnz + cnz/5 + 2*n; -// if (!P || !W || !cs_sprealloc (C, t)) return (cs_idone (P, C, W, 0)) ; - int t = cnz/5 + 2*n; - C.alloc_elements(cnz, t); + int elbow_room = cnz/5 + 2*n; // Request additional # of elements + C.alloc_elements(cnz, elbow_room); // Appand elements + int nzmax = C.m_elements; // Save the maximum # of entries for garbage coll. + C.m_elements -= elbow_room; // Reverse element increase of alloc_elements + + // (Re-)assign C matrix attributes + Cp = C.m_colstart; + int* Ci = C.m_rowinx; + + // Debug + #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING) + cout << " Added elbow room to matrix C " << C << endl; + #endif // Set workspace pointers int* len = wrk_int; @@ -386,7 +534,7 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) int* degree = wrk_int + 5*(n+1); int* w = wrk_int + 6*(n+1); int* hhead = wrk_int + 7*(n+1); - int* last = P; // use P as workspace for last + int* last = P; // use P as workspace for last // ====================================================================== @@ -398,35 +546,38 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) cout << " Step 2: Initialize quotient graph" << endl; #endif - // (Re-)assign C matrix attributes - int nzmax = C.m_elements; // Maximum number of entries - Cp = C.m_colstart; - int* Ci = C.m_rowinx; - // Setup array that contains for each column the number of non-zero elements for (k = 0 ; k < n ; k++) len[k] = Cp[k+1] - Cp[k]; len[n] = 0; + + // Debug + #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING) + cout << " Non-zero elements per column: "; + for (k = 0 ; k <= n ; k++) + cout << " " << len[k]; + cout << endl; + #endif - // Loop over all columns + // Loop over all nodes (i.e. matrix columns) for (i = 0 ; i <= n ; i++) { - head[i] = -1; // degree list i is empty + head[i] = -1; // degree list i is empty last[i] = -1; next[i] = -1; - hhead[i] = -1; // hash list i is empty - nv[i] = 1; // node i is just one node - w[i] = 1; // node i is alive - elen[i] = 0; // Ek of node i is empty - degree[i] = len[i]; // degree of node i + hhead[i] = -1; // hash list i is empty + nv[i] = 1; // node i is just one node + w[i] = 1; // node i is alive + elen[i] = 0; // Ek of node i is empty + degree[i] = len[i]; // degree of node i } // Clear w mark = cs_wclear(0, 0, w, n); // End-point initialisations - elen[n] = -2; // n is a dead element - Cp[n] = -1; // n is a root of assembly tree - w[n] = 0; // n is a dead element + elen[n] = -2; // n is a dead element + Cp[n] = -1; // n is a root of assembly tree + w[n] = 0; // n is a dead element // ====================================================================== @@ -438,7 +589,7 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) cout << " Step 3: Initialize degree lists" << endl; #endif - // Loop over all columns + // Loop over all nodes for (i = 0 ; i < n ; i++) { // Get degree of node @@ -446,25 +597,25 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) // Case A: Node i is empty if (d == 0) { - elen[i] = -2; // element i is dead + elen[i] = -2; // element i is dead nel++ ; - Cp[i] = -1; // i is a root of assembly tree - w[i] = 0; // i is a dead element + Cp[i] = -1; // i is a root of assembly tree + w[i] = 0; // i is a dead element } // Case B: Node is dense else if (d > dense) { - nv[i] = 0; // absorb i into element n - elen[i] = -1; // node i is dead - nel++ ; - Cp[i] = CS_FLIP(n); - nv[n]++ ; + nv[i] = 0; // absorb i into element n + elen[i] = -1; // node i is dead + nel++; + Cp[i] = CS_FLIP(n); + nv[n]++; } // Case C: Node is sparse else { if (head[d] != -1) last[head[d]] = i; - next[i] = head[d]; // put node i in degree list d + next[i] = head[d]; // put node i in degree list d head[d] = i; } } // endfor: looped over all columns @@ -485,65 +636,65 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) // Select node of minimum approximate degree for (k = -1; mindeg < n && (k = head [mindeg]) == -1; mindeg++) ; if (next[k] != -1) last [next[k]] = -1; - head[mindeg] = next[k]; // remove k from degree list - elenk = elen[k]; // elenk = |Ek| - nvk = nv[k]; // # of nodes k represents - nel += nvk; // nv[k] nodes of A eliminated + head[mindeg] = next[k]; // remove k from degree list + elenk = elen[k]; // elenk = |Ek| + nvk = nv[k]; // # of nodes k represents + nel += nvk; // nv[k] nodes of A eliminated - // Garbage collection + // Garbage collection. This frees the elbow rooms in matrix C if (elenk > 0 && cnz + mindeg >= nzmax) { // Loop over all columns - for (j = 0 ; j < n ; j++) { + for (j = 0; j < n; j++) { // If j is a live node or element, then ... if ((p = Cp[j]) >= 0) { - Cp[j] = Ci[p]; // save first entry of object - Ci[p] = CS_FLIP(j); // first entry is now CS_FLIP(j) + Cp[j] = Ci[p]; // save first entry of object + Ci[p] = CS_FLIP(j); // first entry is now CS_FLIP(j) } } // Scan all memory - for (q = 0, p = 0; p < cnz ; ) { + for (q = 0, p = 0; p < cnz; ) { // If we found object j, then ... if ((j = CS_FLIP(Ci[p++])) >= 0) { - Ci[q] = Cp[j]; // restore first entry of object - Cp[j] = q++; // new pointer to object j - for (k3 = 0; k3 < len [j]-1; k3++) + Ci[q] = Cp[j]; // restore first entry of object + Cp[j] = q++; // new pointer to object j + for (k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++]; } } - cnz = q; // Ci[cnz...nzmax-1] now free + cnz = q; // Ci[cnz...nzmax-1] now free } // endif: garbage collection done // Construct new element ---------------------------------------- */ dk = 0; - nv[k] = -nvk; // flag k as in Lk + nv[k] = -nvk; // flag k as in Lk p = Cp[k]; - pk1 = (elenk == 0) ? p : cnz; // do in place if elen[k] == 0 + pk1 = (elenk == 0) ? p : cnz; // do in place if elen[k] == 0 pk2 = pk1; for (k1 = 1; k1 <= elenk + 1; k1++) { // ... if (k1 > elenk) { - e = k; // search the nodes in k - pj = p; // list of nodes starts at Ci[pj] - ln = len [k] - elenk; // length of list of nodes in k + e = k; // search the nodes in k + pj = p; // list of nodes starts at Ci[pj] + ln = len[k] - elenk; // length of list of nodes in k } else { - e = Ci[p++]; // search the nodes in e + e = Ci[p++]; // search the nodes in e pj = Cp[e]; - ln = len[e]; // length of list of nodes in e + ln = len[e]; // length of list of nodes in e } // ... for (k2 = 1 ; k2 <= ln ; k2++) { i = Ci[pj++]; - if ((nvi = nv [i]) <= 0) // Skip of node i is dead or seen + if ((nvi = nv [i]) <= 0) // Skip of node i is dead or seen continue; - dk += nvi; // degree[Lk] += size of node i - nv[i] = -nvi; // negate nv[i] to denote i in Lk - Ci[pk2++] = i; // place i in Lk + dk += nvi; // degree[Lk] += size of node i + nv[i] = -nvi; // negate nv[i] to denote i in Lk + Ci[pk2++] = i; // place i in Lk if (next[i] != -1) last[next[i]] = last[i]; @@ -556,16 +707,16 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) // ... if (e != k) { - Cp[e] = CS_FLIP (k); // absorb e into k - w[e] = 0; // e is now a dead element + Cp[e] = CS_FLIP(k); // absorb e into k + w[e] = 0; // e is now a dead element } } // endfor: k1 - if (elenk != 0) cnz = pk2; // Ci [cnz...nzmax] is free - degree[k] = dk; // external degree of k - |Lk\i| - Cp[k] = pk1; // element k is in Ci[pk1..pk2-1] + if (elenk != 0) cnz = pk2; // Ci [cnz...nzmax] is free + degree[k] = dk; // external degree of k - |Lk\i| + Cp[k] = pk1; // element k is in Ci[pk1..pk2-1] len[k] = pk2 - pk1; - elen[k] = -2; // k is now an element + elen[k] = -2; // k is now an element // Clear w if necessary mark = cs_wclear(mark, lemax, w, n); @@ -573,85 +724,86 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) // Scan 1: Find set differences |Le\Lk| for (pk = pk1 ; pk < pk2 ; pk++) { i = Ci [pk]; - if ((eln = elen [i]) <= 0) // skip if elen[i] empty + if ((eln = elen[i]) <= 0) // skip if elen[i] empty continue; - nvi = -nv[i]; // nv [i] was negated + nvi = -nv[i]; // nv [i] was negated wnvi = mark - nvi; // Scan Ei for (p = Cp[i]; p <= Cp[i] + eln - 1; p++) { - e = Ci [p] ; + e = Ci[p]; if (w[e] >= mark) - w[e] -= nvi; // decrement |Le\Lk| - else if (w[e] != 0) // ensure e is a live element - w[e] = degree[e] + wnvi; // 1st time e seen in scan 1 + w[e] -= nvi; // decrement |Le\Lk| + else if (w[e] != 0) // ensure e is a live element + w[e] = degree[e] + wnvi; // 1st time e seen in scan 1 } } // endfor: scan 1 finished // Scan 2: Degree update for (pk = pk1; pk < pk2; pk++) { - i = Ci [pk] ; /* consider node i in Lk */ - p1 = Cp [i] ; - p2 = p1 + elen [i] - 1 ; - pn = p1 ; + i = Ci[pk]; // consider node i in Lk + p1 = Cp[i]; + p2 = p1 + elen [i] - 1; + pn = p1; // Scan Ei for (h = 0, d = 0, p = p1 ; p <= p2 ; p++) { - e = Ci [p] ; + e = Ci[p]; // If e is an unabsorbed element, then ... - if (w [e] != 0) { - dext = w [e] - mark ; /* dext = |Le\Lk| */ + if (w[e] != 0) { + dext = w[e] - mark; // dext = |Le\Lk| if (dext > 0) { - d += dext ; /* sum up the set differences */ - Ci [pn++] = e ; /* keep e in Ei */ - h += e ; /* compute the hash of node i */ + d += dext; // sum up the set differences + Ci[pn++] = e; // keep e in Ei + h += e; // compute the hash of node i } else { - Cp [e] = CS_FLIP (k) ; /* aggressive absorb. e->k */ - w [e] = 0 ; /* e is a dead element */ + Cp[e] = CS_FLIP(k); // aggressive absorb. e->k + w[e] = 0; // e is a dead element } } } // endfor: scanned Ei // elen[i] = |Ei| - elen [i] = pn - p1 + 1 ; + elen[i] = pn - p1 + 1; p3 = pn ; - p4 = p1 + len [i] ; + p4 = p1 + len[i]; // Prune edges in Ai for (p = p2 + 1 ; p < p4 ; p++) { - j = Ci [p] ; - if ((nvj = nv [j]) <= 0) continue ; /* node j dead or in Lk */ - d += nvj ; /* degree(i) += |j| */ - Ci [pn++] = j ; /* place j in node list of i */ - h += j ; /* compute hash for node i */ + j = Ci[p]; + if ((nvj = nv[j]) <= 0) // node j dead or in Lk + continue; + d += nvj; // degree(i) += |j| + Ci[pn++] = j; // place j in node list of i + h += j; // compute hash for node i } // Check for mass elimination if (d == 0) { - Cp [i] = CS_FLIP (k) ; /* absorb i into k */ - nvi = -nv [i] ; - dk -= nvi ; /* |Lk| -= |i| */ - nvk += nvi ; /* |k| += nv[i] */ - nel += nvi ; - nv [i] = 0 ; - elen [i] = -1 ; /* node i is dead */ + Cp[i] = CS_FLIP(k); // absorb i into k + nvi = -nv[i]; + dk -= nvi; // |Lk| -= |i| + nvk += nvi; // |k| += nv[i] + nel += nvi; + nv[i] = 0; + elen[i] = -1; // node i is dead } else { - degree [i] = CS_MIN (degree [i], d) ; /* update degree(i) */ - Ci [pn] = Ci [p3] ; /* move first node to end */ - Ci [p3] = Ci [p1] ; /* move 1st el. to end of Ei */ - Ci [p1] = k ; /* add k as 1st element in of Ei */ - len [i] = pn - p1 + 1 ; /* new len of adj. list of node i */ - h %= n ; /* finalize hash of i */ - next [i] = hhead [h] ; /* place i in hash bucket */ - hhead [h] = i ; - last [i] = h ; /* save hash of i in last[i] */ + degree[i] = CS_MIN(degree[i], d); // update degree(i) + Ci[pn] = Ci[p3]; // move first node to end + Ci[p3] = Ci[p1]; // move 1st el. to end of Ei + Ci[p1] = k; // add k as 1st element in of Ei + len[i] = pn - p1 + 1; // new len of adj. list of node i + h %= n; // finalize hash of i + next[i] = hhead[h]; // place i in hash bucket + hhead[h] = i; + last[i] = h; // save hash of i in last[i] } } // endfor: scan 2 is done // finalize |Lk| - degree[k] = dk ; + degree[k] = dk; lemax = CS_MAX(lemax, dk); // Clear w @@ -659,33 +811,35 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) // Supernode detection for (pk = pk1 ; pk < pk2 ; pk++) { - i = Ci [pk] ; - if (nv [i] >= 0) continue ; /* skip if i is dead */ - h = last [i] ; /* scan hash bucket of node i */ - i = hhead [h] ; - hhead [h] = -1 ; /* hash bucket will be empty */ + i = Ci[pk]; + if (nv[i] >= 0) continue; // skip if i is dead + h = last[i]; // scan hash bucket of node i + i = hhead[h]; + hhead[h] = -1; // hash bucket will be empty // ... - for ( ; i != -1 && next [i] != -1 ; i = next [i], mark++) { - ln = len [i] ; - eln = elen [i] ; - for (p = Cp [i]+1 ; p <= Cp [i] + ln-1 ; p++) w [Ci [p]] = mark; - jlast = i ; + for ( ; i != -1 && next[i] != -1 ; i = next[i], mark++) { + ln = len[i]; + eln = elen[i]; + for (p = Cp[i]+1 ; p <= Cp[i] + ln-1 ; p++) + w[Ci[p]] = mark; + jlast = i; // Compare i with all j - for (j = next [i] ; j != -1 ; ) { - ok = (len [j] == ln) && (elen [j] == eln) ; - for (p = Cp [j] + 1 ; ok && p <= Cp [j] + ln - 1 ; p++) { - if (w [Ci [p]] != mark) ok = 0 ; /* compare i and j*/ + for (j = next[i]; j != -1; ) { + ok = (len[j] == ln) && (elen[j] == eln) ; + for (p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++) { + if (w[Ci[p]] != mark) // compare i and j + ok = 0; } // i and j are identical if (ok) { - Cp[j] = CS_FLIP(i); /* absorb j into i */ + Cp[j] = CS_FLIP(i); // absorb j into i nv[i] += nv[j]; nv[j] = 0; - elen[j] = -1; /* node j is dead */ - j = next[j]; /* delete j from hash bucket */ + elen[j] = -1; // node j is dead + j = next[j]; // delete j from hash bucket next[jlast] = j; } // j and i are different @@ -697,57 +851,75 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) } // } // endfor: supernode detection - // Finalize new element - for (p = pk1, pk = pk1 ; pk < pk2 ; pk++) /* finalize Lk */ - { - i = Ci [pk] ; - if ((nvi = -nv [i]) <= 0) continue ;/* skip if i is dead */ - nv [i] = nvi ; /* restore nv[i] */ - d = degree [i] + dk - nvi ; /* compute external degree(i) */ - d = CS_MIN (d, n - nel - nvi) ; - if (head [d] != -1) last [head [d]] = i ; - next [i] = head [d] ; /* put i back in degree list */ - last [i] = -1 ; - head [d] = i ; - mindeg = CS_MIN (mindeg, d) ; /* find new minimum degree */ - degree [i] = d ; - Ci [p++] = i ; /* place i in Lk */ + // Finalize new element Lk + for (p = pk1, pk = pk1; pk < pk2; pk++) { + i = Ci[pk]; + if ((nvi = -nv[i]) <= 0) // skip if i is dead + continue; + nv[i] = nvi; // restore nv[i] + d = degree[i] + dk - nvi; // compute external degree(i) + d = CS_MIN(d, n - nel - nvi); + if (head[d] != -1) last[head[d]] = i; + next [i] = head [d]; // put i back in degree list + last [i] = -1 ; + head [d] = i ; + mindeg = CS_MIN (mindeg, d); // find new minimum degree + degree [i] = d ; + Ci [p++] = i ; // place i in Lk } - nv [k] = nvk ; /* # nodes absorbed into k */ - if ((len [k] = p-pk1) == 0) /* length of adj list of element k*/ - { - Cp [k] = -1 ; /* k is a root of the tree */ - w [k] = 0 ; /* k is now a dead element */ + + // Number of nodes absorbed into k + nv[k] = nvk; + + // ... + if ((len [k] = p-pk1) == 0) { // length of adj list of element k + Cp [k] = -1 ; // k is a root of the tree + w [k] = 0 ; // k is now a dead element } - if (elenk != 0) cnz = p ; /* free unused space in Lk */ - } + if (elenk != 0) cnz = p ; // free unused space in Lk + + } // endwhile + + + // ====================================================================== + // Step 5: Postordering + // ====================================================================== + + // Debug + #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING) + cout << " Step 5: Postordering" << endl; + #endif - // Postordering - for (i = 0 ; i < n ; i++) Cp [i] = CS_FLIP (Cp [i]) ;/* fix assembly tree */ - for (j = 0 ; j <= n ; j++) head [j] = -1 ; - for (j = n ; j >= 0 ; j--) /* place unordered nodes in lists */ - { - if (nv [j] > 0) continue ; /* skip if j is an element */ - next [j] = head [Cp [j]] ; /* place j in list of its parent */ - head [Cp [j]] = j ; - } - for (e = n ; e >= 0 ; e--) /* place elements in lists */ - { - if (nv [e] <= 0) continue ; /* skip unless e is an element */ - if (Cp [e] != -1) - { - next [e] = head [Cp [e]] ; /* place e in list of its parent */ - head [Cp [e]] = e ; + // Fix assembly tree + for (i = 0; i < n; i++) + Cp[i] = CS_FLIP(Cp[i]); + for (j = 0 ; j <= n ; j++) + head[j] = -1; + + // Place unordered nodes in lists + for (j = n ; j >= 0 ; j--) { + if (nv[j] > 0) continue; // skip if j is an element + next[j] = head [Cp[j]]; // place j in list of its parent + head[Cp[j]] = j; + } + + // Place elements in lists + for (e = n ; e >= 0 ; e--) { + if (nv[e] <= 0) continue; // skip unless e is an element + if (Cp[e] != -1) { + next[e] = head[Cp[e]]; // place e in list of its parent + head[Cp[e]] = e; } } + // ====================================================================== - // Step 5: Postorder the assembly tree + // Step 6: Postorder the assembly tree // ====================================================================== // Debug #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING) - cout << " Step 5: Postorder the assembly tree" << endl; + cout << " Step 6: Postorder the assembly tree" << endl; #endif // Postorder the assembly tree @@ -769,7 +941,7 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) #endif // Return result - return (P); + return P; } @@ -784,7 +956,9 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A) * ata 0: count LL'=A, 1: count LL'=A'A * * Output: p[0..n] integer array of n+1 elements * ***************************************************************************/ -int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent, +#define HEAD(k,j) (ata ? head [k] : j) +#define NEXT(J) (ata ? next [J] : -1) +int* GSparseSymbolic::cs_counts(const GSparseMatrix* A, const int* parent, const int* post, int ata) { // Declare loop variables @@ -792,7 +966,7 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent, // Return NULL pointer if input arrays are NULL or sparse matrix is // empty - if (!parent || !post) return (NULL); + if (!parent || !post) return NULL; // Assign input matrix attributes int m = A->m_rows; @@ -801,8 +975,9 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent, // Allocate result int* colcount = new int[n]; if (colcount == NULL) - throw mem_alloc("GSparseSymbolic::cs_counts(GSparseMatrix*, const int*, const int*, int)", - n); + throw mem_alloc( + "GSparseSymbolic::cs_counts(GSparseMatrix*, const int*, const int*, int)", + n); int* delta = colcount; // Allocate workspace @@ -812,8 +987,8 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent, throw mem_alloc("GSparseSymbolic::cs_counts(GSparseMatrix*, const int*, const int*, int)", wrk_size); - // Get transpose of A: AT = A' - GSparseMatrix AT = A->cs_transpose(A,0); + // Get (logical) transpose of A: AT = A' + GSparseMatrix AT = cs_transpose(*A, 0); // Set-up pointers to workspace int* ancestor = wrk_int; @@ -836,7 +1011,9 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent, int* ATp = AT.m_colstart; int* ATi = AT.m_rowinx; - // Initialise + // Initialise. Note that init_ata wires the working array into the *head + // and *next pointers, so both arrays are just aliases of the wrk_int + // array. int* head = NULL; int* next = NULL; if (ata) init_ata(&AT, post, wrk_int, &head, &next); @@ -845,8 +1022,6 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent, for (i = 0; i < n; i++) ancestor[i] = i; // Loop over all columns of matrix A (or rows of matrix AT) - #define HEAD(k,j) (ata ? head [k] : j) - #define NEXT(J) (ata ? next [J] : -1) for (k = 0; k < n; k++) { j = post[k]; // j is the kth node in postordered etree if (parent[j] != -1) delta[parent[j]]--; // j is not a root @@ -872,7 +1047,7 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent, delete [] wrk_int; // Return result - return (colcount); + return colcount; } @@ -885,7 +1060,7 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent, * ata flag * * Output: etree(A) Evaluation tree of matrix A * ***************************************************************************/ -int* GSparseSymbolic::cs_etree(GSparseMatrix* A, int ata) +int* GSparseSymbolic::cs_etree(const GSparseMatrix* A, int ata) { // Declare loop variables int i, k, p; @@ -910,7 +1085,10 @@ int* GSparseSymbolic::cs_etree(GSparseMatrix* A, int ata) int* prev = wrk_int + n; // If 'prev' is requested, initialise array with -1 - if (ata) for (i = 0; i < m ; ++i) prev[i] = -1; + if (ata) { + for (i = 0; i < m ; ++i) + prev[i] = -1; + } // Loop over all nodes (columns) for (k = 0; k < n; ++k) { @@ -940,7 +1118,7 @@ int* GSparseSymbolic::cs_etree(GSparseMatrix* A, int ata) delete [] wrk_int; // Return result - return (parent); + return parent; } @@ -956,7 +1134,8 @@ int* GSparseSymbolic::cs_etree(GSparseMatrix* A, int ata) * other Other function * * Output: result 0..n-1 if ok, -1 in case of error * ***************************************************************************/ -int GSparseSymbolic::cs_fkeep(GSparseMatrix* A, int(*fkeep)(int, int, double, void*), +int GSparseSymbolic::cs_fkeep(GSparseMatrix* A, + int(*fkeep)(int, int, double, void*), void* other) { // Return error if some of the input pointers is invalid @@ -984,13 +1163,15 @@ int GSparseSymbolic::cs_fkeep(GSparseMatrix* A, int(*fkeep)(int, int, double, vo } // Finalise A - Ap[n] = nz; + Ap[n] = nz; + A->m_elements = nz; // Remove extra space from A -// cs_sprealloc (A, 0) ; /* remove extra space from A */ +// cs_sprealloc(A, 0); + A->free_elements(nz, (A->m_elements-nz)); // Return number of non-zero elements - return (nz) ; + return nz; } @@ -1061,7 +1242,7 @@ int GSparseSymbolic::cs_leaf(int i, int j, const int *first, int *maxfirst, int* GSparseSymbolic::cs_pinv(int const *p, int n) { // Return NULL pointer if input pointer is NULL. This denotes identity - if (!p || n < 1) return (NULL); + if (!p || n < 1) return NULL; // Allocate result array int* pinv = new int[n]; @@ -1073,7 +1254,7 @@ int* GSparseSymbolic::cs_pinv(int const *p, int n) pinv[p[k]] = k; // Return result - return (pinv); + return pinv; } @@ -1137,90 +1318,7 @@ int* GSparseSymbolic::cs_post(const int* parent, int n) delete [] wrk_int; // Return result - return (post); - -} - - -/*************************************************************************** - * cs_symperm * - * ----------------------------------------------------------------------- * - * C = A(p,p) where A and C are symmetric the upper part stored. * - * pinv not p * - * ----------------------------------------------------------------------- * - * Input: A Sparse matrix * - * pinv pinv[0..n-1] * - * values =1 allocate memory for values (NOT USED) * - * Output: C Sparse matrix * - ***************************************************************************/ -GSparseMatrix GSparseSymbolic::cs_symperm(GSparseMatrix* A, const int* pinv, - int values) const -{ - // Declare loop variables - int i, j, p, q, i2, j2; - - // Assign matrix attributes - int n = A->m_cols; - int* Ap = A->m_colstart; - int* Ai = A->m_rowinx; - double* Ax = A->m_data; - - // Allocate result matrix - GSparseMatrix C(n, n, Ap[n]); - - // Allocate and initialise workspace - int wrk_size = n; - int* wrk_int = new int[wrk_size]; - if (wrk_int == NULL) - throw mem_alloc("GSparseSymbolic::cs_symperm(GSparseMatrix*, const int*, int", - wrk_size); - for (i = 0; i < wrk_size; ++i) - wrk_int[i] = 0; - - // Assign result matriix attributes - int* Cp = C.m_colstart; - int* Ci = C.m_rowinx; - double* Cx = C.m_data; - - // Count entries in each column of C - for (j = 0; j < n; j++) { - - // Column j of A is column j2 of C - j2 = pinv ? pinv[j] : j; - - // Loop over entries in column j - for (p = Ap[j]; p < Ap[j+1]; p++) { - i = Ai [p]; - if (i > j) continue; // skip lower triangular part of A - i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C - wrk_int[CS_MAX(i2, j2)]++; // column count of C - } - } - - // Compute column pointers of C - cs_cumsum(Cp, wrk_int, n); - - // Loop over all columns of A - for (j = 0 ; j < n ; j++) { - - // Column j of A is column j2 of C - j2 = pinv ? pinv[j] : j; - - // Loop over entries in column j - for (p = Ap[j]; p < Ap[j+1]; p++) { - i = Ai [p] ; - if (i > j) continue; // skip lower triangular part of A - i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C - Ci[q = wrk_int[CS_MAX(i2,j2)]++] = CS_MIN(i2,j2); - Cx[q] = Ax[p]; - } - } - - // Free workspace - delete [] wrk_int; - - // Return result - return (C); + return post; } @@ -1266,7 +1364,7 @@ int GSparseSymbolic::cs_tdfs(int j, int k, int* head, const int* next, } // Return result - return (k) ; + return k; } @@ -1290,7 +1388,7 @@ int GSparseSymbolic::cs_tdfs(int j, int k, int* head, const int* next, * next pointer to next[0..n-1] linked list info * * Output: - - * ***************************************************************************/ -void GSparseSymbolic::init_ata(GSparseMatrix* AT, const int* post, +void GSparseSymbolic::init_ata(const GSparseMatrix* AT, const int* post, int* wrk_int, int** head, int** next) { // Declare loop variables @@ -1321,6 +1419,9 @@ void GSparseSymbolic::init_ata(GSparseMatrix* AT, const int* post, (*next)[i] = (*head)[k]; (*head)[k] = i; } + + // Return + return; } diff --git a/src/GSparseSymbolic.hpp b/src/GSparseSymbolic.hpp index d83e9c2..43707f0 100644 --- a/src/GSparseSymbolic.hpp +++ b/src/GSparseSymbolic.hpp @@ -61,7 +61,7 @@ public: // Unary operators // Functions - void cholesky_symbolic_analysis(int order, GSparseMatrix& m); + void cholesky_symbolic_analysis(int order, const GSparseMatrix& m); // Exception: Invalid ordering scheme class invalid_order : public GException { @@ -71,18 +71,17 @@ public: private: // Functions - int* cs_amd (int order, GSparseMatrix* A); - int* cs_counts(GSparseMatrix* A, const int* parent, const int* post, int ata); - int* cs_etree(GSparseMatrix* A, int ata); - int cs_fkeep(GSparseMatrix* A, int(*fkeep)(int, int, double, void*), void* other); - int cs_leaf(int i, int j, const int* first, int* maxfirst, int* prevleaf, int* ancestor, int* jleaf); - int* cs_pinv(int const* p, int n); - int* cs_post(const int* parent, int n); - GSparseMatrix cs_symperm(GSparseMatrix* A, const int* pinv, int values) const; - int cs_tdfs(int j, int k, int* head, const int* next, int* post, int* stack); - static void init_ata(GSparseMatrix* AT, const int* post, int* wrk_int, int** head, int** next); - static int cs_diag(int i, int j, double aij, void* other); - static int cs_wclear(int mark, int lemax, int* w, int n); + int* cs_amd (int order, const GSparseMatrix* A); + int* cs_counts(const GSparseMatrix* A, const int* parent, const int* post, int ata); + int* cs_etree(const GSparseMatrix* A, int ata); + int cs_fkeep(GSparseMatrix* A, int(*fkeep)(int, int, double, void*), void* other); + int cs_leaf(int i, int j, const int* first, int* maxfirst, int* prevleaf, int* ancestor, int* jleaf); + int* cs_pinv(int const* p, int n); + int* cs_post(const int* parent, int n); + int cs_tdfs(int j, int k, int* head, const int* next, int* post, int* stack); + static void init_ata(const GSparseMatrix* AT, const int* post, int* wrk_int, int** head, int** next); + static int cs_diag(int i, int j, double aij, void* other); + static int cs_wclear(int mark, int lemax, int* w, int n); // Data int* m_pinv; // Inverse row permutation for QR, fill reduce permutation for Cholesky diff --git a/src/GVector.cpp b/src/GVector.cpp index 0078723..9e90125 100644 --- a/src/GVector.cpp +++ b/src/GVector.cpp @@ -46,15 +46,32 @@ GVector::GVector(int num) /*************************************************************************** + * GVector copy constructor * + ***************************************************************************/ +GVector::GVector(const GVector& v) +{ + // Allocate vector array. Throw exception if allocation failed + m_data = new double[v.m_num]; + if (m_data == NULL) + throw mem_alloc("GVector copy constructor", v.m_num); + + // Store vector size and copy elements + m_num = v.m_num; + for (int i = 0; i < m_num; ++i) + m_data[i] = v.m_data[i]; + + // Return + return; +} + + +/*************************************************************************** * GVector destructor * ***************************************************************************/ GVector::~GVector() { // Deallocate only if memory has indeed been allocated - if (m_num > 0) { - m_num = 0; - delete[] m_data; - } + if (m_data != NULL) delete[] m_data; // Return return; @@ -67,6 +84,33 @@ GVector::~GVector() = = ==========================================================================*/ +/*************************************************************************** + * GVector assignment operator * + ***************************************************************************/ +GVector& GVector::operator= (const GVector& v) +{ + // Assign only if not identical + if (this != &v) { + + // First delete old data if it exists + if (m_data != NULL) delete [] m_data; + + // Allocate vector array. Throw exception if allocation failed + m_data = new double[v.m_num]; + if (m_data == NULL) + throw mem_alloc("GVector assignment operator", v.m_num); + + // Store vector size and copy elements + m_num = v.m_num; + for (int i = 0; i < m_num; ++i) + m_data[i] = v.m_data[i]; + } + + // Return vector + return *this; +} + + /*========================================================================== = = = GVector member functions = diff --git a/src/GVector.hpp b/src/GVector.hpp index 5aa0b2c..be2f783 100644 --- a/src/GVector.hpp +++ b/src/GVector.hpp @@ -134,18 +134,6 @@ private: * Inline members * ***************************************************************************/ -// Copy constructor -inline -GVector::GVector(const GVector& v) -{ - m_data = new double[v.m_num]; - if (m_data == NULL) - throw mem_alloc("GVector copy constructor", v.m_num); - m_num = v.m_num; - for (int i = 0; i < m_num; ++i) - m_data[i] = v.m_data[i]; -} - // Vector element access operator inline double& GVector::operator() (int inx) @@ -168,22 +156,6 @@ const double& GVector::operator() (int inx) const return m_data[inx]; } -// Vector assignment operator -inline -GVector& GVector::operator= (const GVector& v) -{ - if (this != &v) { - if (m_data != NULL) delete [] m_data; - m_data = new double[v.m_num]; - if (m_data == NULL) - throw mem_alloc("GVector assignment operator", v.m_num); - m_num = v.m_num; - for (int i = 0; i < m_num; ++i) - m_data[i] = v.m_data[i]; - } - return *this; -} - // Vector unary addition operator inline GVector& GVector::operator+= (const GVector& v)