Jürgen Knödlseder / gammalib

Commit cb738b250bd32e8d741c8ebe4928cf70a0175c3a

Authored by knodlseder
1 parent 9dc93639

Correct small bugs and make spi_toolslib compliant

Inline Side-by-side
Showing 9 changed files with 725 additions and 498 deletions
src/GMatrix.cpp
  View file @cb738b2
... ... @@ -170,9 +170,9 @@ GMatrix::GMatrix(const GMatrix& m)
170 170 GMatrix::~GMatrix()
171 171 {
172 172 // De-allocate only if memory has indeed been allocated
173   - if (m_colstart != NULL) delete [] m_colstart;
174 173 if (m_rowsel != NULL) delete [] m_rowsel;
175 174 if (m_colsel != NULL) delete [] m_colsel;
  175 + if (m_colstart != NULL) delete [] m_colstart;
176 176 if (m_data != NULL) delete [] m_data;
177 177  
178 178 // Return
... ...
src/GSparseMatrix.cpp
  View file @cb738b2
... ... @@ -48,8 +48,8 @@ GSparseMatrix::GSparseMatrix(int rows, int cols, int elements) : GMatrix()
48 48 // Initialise private members for clean destruction
49 49 m_rows = 0;
50 50 m_cols = 0;
51   - m_elements = 0;
52   - m_alloc = 0;
  51 + m_elements = 0; // Logically used # of elements
  52 + m_alloc = 0; // Allocated # of elements (>= m_elements)
53 53 m_data = NULL;
54 54 m_colstart = NULL;
55 55 m_rowinx = NULL;
... ... @@ -537,7 +537,14 @@ double GSparseMatrix::sum() const
537 537 ***************************************************************************/
538 538 void GSparseMatrix::transpose()
539 539 {
540   - *this = cs_transpose(this, 1);
  540 + // Fill pending value in input matrix
  541 + fill_pending();
  542 +
  543 + // Transpose matrix
  544 + *this = cs_transpose(*this, 1);
  545 +
  546 + // Return
  547 + return;
541 548 }
542 549  
543 550  
... ... @@ -803,8 +810,12 @@ void GSparseMatrix::cholesky_decompose(int compress)
803 810 int matrix_rows = m_rows;
804 811 int matrix_cols = m_cols;
805 812  
806   - // Delete any existing symbolic analysis object
  813 + // Delete any existing symbolic and numeric analysis object and reset
  814 + // pointers
807 815 if (m_symbolic != NULL) delete (GSparseSymbolic*)m_symbolic;
  816 + if (m_numeric != NULL) delete (GSparseNumeric*)m_numeric;
  817 + m_symbolic = NULL;
  818 + m_numeric = NULL;
808 819  
809 820 // Allocate symbolic analysis object
810 821 GSparseSymbolic* symbolic = new GSparseSymbolic();
... ... @@ -827,44 +838,26 @@ void GSparseMatrix::cholesky_decompose(int compress)
827 838 // Ordering an symbolic analysis of matrix. This sets up an array 'pinv'
828 839 // which contains the fill-in reducing permutations
829 840 symbolic->cholesky_symbolic_analysis(1, *this);
  841 +
  842 + // Store symbolic pointer in sparse matrix object
  843 + m_symbolic = (void*)symbolic;
830 844  
831 845 // Perform numeric Cholesky decomposition
832 846 numeric.cholesky_numeric_analysis(*this, *symbolic);
833   -
834   - // Insert zero rows and columns if they have been removed previously.
835   - // Since no selection information is present in numeric.m_L we have to
836   - // copy it from this in order to allow for insertion.
837   - if (compress) {
838   -
839   - // If there is a row selection then copy it
840   - if (m_rowsel != NULL && m_num_rowsel > 0) {
841   - numeric.m_L->m_rowsel = new int[m_num_rowsel];
842   - if (numeric.m_L->m_rowsel == NULL)
843   - throw mem_alloc("GSparseMatrix::cholesky_decompose(int)", m_num_rowsel);
844   - for (int i = 0; i < m_num_rowsel; ++i)
845   - numeric.m_L->m_rowsel[i] = m_rowsel[i];
846   - numeric.m_L->m_num_rowsel = m_num_rowsel;
847   - }
848   -
849   - // If there is a column selection then copy it
850   - if (m_colsel != NULL && m_num_colsel > 0) {
851   - numeric.m_L->m_colsel = new int[m_num_colsel];
852   - if (numeric.m_L->m_colsel == NULL)
853   - throw mem_alloc("GSparseMatrix::cholesky_decompose(int)", m_num_colsel);
854   - for (int i = 0; i < m_num_colsel; ++i)
855   - numeric.m_L->m_colsel[i] = m_colsel[i];
856   - numeric.m_L->m_num_colsel = m_num_colsel;
857   - }
858   -
859   - // Now we can insert zero rows and columns
860   - numeric.m_L->insert_zero_row_col(matrix_rows, matrix_cols);
  847 +
  848 + // Copy L matrix into this object
  849 + free_elements(0, m_elements);
  850 + alloc_elements(0, numeric.m_L->m_elements);
  851 + for (int i = 0; i < m_elements; ++i) {
  852 + m_data[i] = numeric.m_L->m_data[i];
  853 + m_rowinx[i] = numeric.m_L->m_rowinx[i];
861 854 }
  855 + for (int col = 0; col <= m_cols; ++col)
  856 + m_colstart[col] = numeric.m_L->m_colstart[col];
862 857  
863   - // Store Cholesky decomposition in matrix
864   - *this = *(numeric.m_L);
865   -
866   - // Store symbolic pointer in sparse matrix object
867   - m_symbolic = (void*)symbolic;
  858 + // Insert zero rows and columns if they have been removed previously.
  859 + if (compress)
  860 + insert_zero_row_col(matrix_rows, matrix_cols);
868 861  
869 862 // Return
870 863 return;
... ... @@ -1534,7 +1527,7 @@ void GSparseMatrix::remove_zero_row_col(void)
1534 1527  
1535 1528 // Dump column start array
1536 1529 #if defined(G_DEBUG_SPARSE_COMPRESSION)
1537   - cout << " before compression:";
  1530 + cout << " before compression (" << m_colstart[m_cols] << "):";
1538 1531 for (int col = 0; col <= m_cols; ++col)
1539 1532 cout << " " << m_colstart[col];
1540 1533 cout << endl;
... ... @@ -1567,7 +1560,7 @@ void GSparseMatrix::remove_zero_row_col(void)
1567 1560  
1568 1561 // Dump column start array
1569 1562 #if defined(G_DEBUG_SPARSE_COMPRESSION)
1570   - cout << " after compression :";
  1563 + cout << " after compression (" << m_colstart[m_num_colsel] << ") :";
1571 1564 for (int c_col = 0; c_col <= m_num_colsel; ++c_col)
1572 1565 cout << " " << m_colstart[c_col];
1573 1566 cout << endl;
... ... @@ -1576,10 +1569,9 @@ void GSparseMatrix::remove_zero_row_col(void)
1576 1569 // Free row mapping array
1577 1570 delete [] row_map;
1578 1571  
1579   - // Update matrix attributes
1580   - m_rows = m_num_rowsel;
1581   - m_cols = m_num_colsel;
1582   - m_elements = m_colstart[m_num_colsel];
  1572 + // Update matrix size
  1573 + m_rows = m_num_rowsel;
  1574 + m_cols = m_num_colsel;
1583 1575  
1584 1576 // Return
1585 1577 return;
... ... @@ -1616,7 +1608,7 @@ void GSparseMatrix::insert_zero_row_col(int rows, int cols)
1616 1608  
1617 1609 // Dump column start array
1618 1610 #if defined(G_DEBUG_SPARSE_COMPRESSION)
1619   - cout << " before restoration:";
  1611 + cout << " before restoration (" << m_colstart[m_num_colsel] << "):";
1620 1612 for (int c_col = 0; c_col <= m_num_colsel; ++c_col)
1621 1613 cout << " " << m_colstart[c_col];
1622 1614 cout << endl;
... ... @@ -1624,7 +1616,11 @@ void GSparseMatrix::insert_zero_row_col(int rows, int cols)
1624 1616  
1625 1617 // If column selection exists then restore column counters
1626 1618 if (m_colsel != NULL) {
  1619 +
  1620 + // Start insertion from the last original column
1627 1621 int col_stop = cols - 1;
  1622 +
  1623 + // Loop over all compressed columns
1628 1624 for (int c_col = m_num_colsel-1; c_col > 0; --c_col) {
1629 1625 int col_start = m_colsel[c_col-1] + 1;
1630 1626 int col_value = m_colstart[c_col];
... ... @@ -1632,18 +1628,24 @@ void GSparseMatrix::insert_zero_row_col(int rows, int cols)
1632 1628 m_colstart[col] = col_value;
1633 1629 col_stop = col_start - 1;
1634 1630 }
  1631 +
  1632 + // Set first columns if they are not yet set
  1633 + for (int col = 0; col <= col_stop; ++col)
  1634 + m_colstart[col] = 0;
  1635 +
  1636 + // Restore the number of elements
1635 1637 m_colstart[cols] = m_elements;
1636 1638 }
1637 1639  
1638 1640 // Dump column start array
1639 1641 #if defined(G_DEBUG_SPARSE_COMPRESSION)
1640   - cout << " after restoration :";
  1642 + cout << " after restoration (" << m_colstart[cols] << ") :";
1641 1643 for (int col = 0; col <= cols; ++col)
1642 1644 cout << " " << m_colstart[col];
1643 1645 cout << endl;
1644 1646 #endif
1645 1647  
1646   - // Update matrix attributes
  1648 + // Update matrix size
1647 1649 m_rows = rows;
1648 1650 m_cols = cols;
1649 1651  
... ... @@ -1652,63 +1654,6 @@ void GSparseMatrix::insert_zero_row_col(int rows, int cols)
1652 1654 }
1653 1655  
1654 1656  
1655   -/***************************************************************************
1656   - * cs_transpose *
1657   - * ----------------------------------------------------------------------- *
1658   - * Transpose matrix. The flag 'values' allows to avoid copying the actual *
1659   - * data values. This allows to perform a logical matrix transposition, as *
1660   - * needed by the symbolic matrix analysis class. *
1661   - ***************************************************************************/
1662   -GSparseMatrix GSparseMatrix::cs_transpose(GSparseMatrix* m, int values)
1663   -{
1664   - // Fill pending value in input matrix
1665   - m->fill_pending();
1666   -
1667   - // Declare and allocate result matrix
1668   - GSparseMatrix result(m->m_cols, m->m_rows, m->m_elements);
1669   -
1670   - // Allocate and initialise workspace
1671   - int wrk_size = m->m_rows;
1672   - int* wrk_int = new int[wrk_size];
1673   - if (wrk_int == NULL)
1674   - throw mem_alloc("GSparseMatrix::cs_transpose(GSparseMatrix*, int)",
1675   - wrk_size);
1676   - for (int i = 0; i < wrk_size; i++)
1677   - wrk_int[i] = 0;
1678   -
1679   - // Setup the number of non-zero elements in each row
1680   - for (int p = 0; p < m->m_elements; p++)
1681   - wrk_int[m->m_rowinx[p]]++;
1682   -
1683   - // Set row pointers. To use a GSparseSymbolic function we have to
1684   - // allocate and object (but this does not take memory)
1685   - cs_cumsum(result.m_colstart, wrk_int, m->m_rows);
1686   -
1687   - // Case A: Normal transponse, including assignment of values
1688   - if (values) {
1689   - for (int col = 0; col < m->m_cols; col++) {
1690   - for (int p = m->m_colstart[col] ; p < m->m_colstart[col+1] ; p++) {
1691   - int i = wrk_int[m->m_rowinx[p]]++;
1692   - result.m_rowinx[i] = col;
1693   - result.m_data[i] = m->m_data[p] ;
1694   - }
1695   - }
1696   - }
1697   -
1698   - // Case B: Logical transponse, no assignment of values is performed
1699   - else {
1700   - for (int col = 0; col < m->m_cols; col++) {
1701   - for (int p = m->m_colstart[col] ; p < m->m_colstart[col+1] ; p++)
1702   - result.m_rowinx[wrk_int[m->m_rowinx[p]]++] = col;
1703   - }
1704   - }
1705   -
1706   - // Return transponse matrix
1707   - return result;
1708   -
1709   -}
1710   -
1711   -
1712 1657 /*==========================================================================
1713 1658 = =
1714 1659 = GSparseMatrix friends =
... ... @@ -1732,42 +1677,36 @@ ostream&amp; operator&lt;&lt; (ostream&amp; os, const GSparseMatrix&amp; m)
1732 1677 m.m_num_colsel << ")" << endl;
1733 1678 else
1734 1679 os << " Number of columns .........: " << m.m_cols << endl;
1735   - os << " Number of non-zero elements: " << m.m_elements << endl;
  1680 + os << " Number of non-zero elements: " << m.m_colstart[m.m_cols] << " (" <<
  1681 + m.m_elements << ")" << endl;
1736 1682 os << " Number of allocated cells .: " << m.m_alloc << endl;
1737 1683 os << " Memory block size .........: " << m.m_mem_block << endl;
1738 1684 os << " Sparse matrix fill ........: " << m.fill() << endl;
1739 1685  
1740   - // Loop over all matrix elements and put them into the stream (only if
1741   - // their number is not too large)
1742   - if (m.m_rows < 20 && m.m_cols < 20) {
1743   - for (int row = 0; row < m.m_rows; ++row) {
1744   - os << " ";
1745   - for (int col = 0; col < m.m_cols; ++col) {
1746   - os << m(row,col);
1747   - if (col != m.m_cols-1)
1748   - os << ", ";
1749   - }
1750   - if (row != m.m_rows-1)
1751   - os << endl;
1752   - }
  1686 + // Loop over all matrix elements and put them into the stream
  1687 + for (int row = 0; row < m.m_rows; ++row) {
  1688 + os << " ";
  1689 + for (int col = 0; col < m.m_cols; ++col) {
  1690 + os << m(row,col);
  1691 + if (col != m.m_cols-1)
  1692 + os << ", ";
  1693 + }
  1694 + if (row != m.m_rows-1)
  1695 + os << endl;
1753 1696 }
1754 1697  
1755 1698 // If there is a row compression the show scheme
1756 1699 if (m.m_rowsel != NULL) {
1757 1700 os << endl << " Row selection ..:";
1758   - if (m.m_num_rowsel < 20) {
1759   - for (int row = 0; row < m.m_num_rowsel; ++row)
1760   - os << " " << m.m_rowsel[row];
1761   - }
  1701 + for (int row = 0; row < m.m_num_rowsel; ++row)
  1702 + os << " " << m.m_rowsel[row];
1762 1703 }
1763 1704  
1764 1705 // If there is a column compression the show scheme
1765 1706 if (m.m_colsel != NULL) {
1766 1707 os << endl << " Column selection:";
1767   - if (m.m_num_colsel < 20) {
1768   - for (int col = 0; col < m.m_num_colsel; ++col)
1769   - os << " " << m.m_colsel[col];
1770   - }
  1708 + for (int col = 0; col < m.m_num_colsel; ++col)
  1709 + os << " " << m.m_colsel[col];
1771 1710 }
1772 1711  
1773 1712 // If there is a symbolic decomposition then put it also in the stream
... ... @@ -1805,6 +1744,145 @@ GSparseMatrix fabs(const GSparseMatrix&amp; m)
1805 1744 }
1806 1745  
1807 1746  
  1747 +/***************************************************************************
  1748 + * cs_symperm *
  1749 + * ----------------------------------------------------------------------- *
  1750 + * C = A(p,p) where A and C are symmetric the upper part stored. *
  1751 + * ----------------------------------------------------------------------- *
  1752 + * Input: A Sparse matrix *
  1753 + * pinv pinv[0..n-1] *
  1754 + * Output: C Sparse matrix *
  1755 + ***************************************************************************/
  1756 +GSparseMatrix cs_symperm(const GSparseMatrix& m, const int* pinv)
  1757 +{
  1758 + // Declare loop variables
  1759 + int i, j, p, q, i2, j2;
  1760 +
  1761 + // Assign matrix attributes
  1762 + int n = m.m_cols;
  1763 + int* Ap = m.m_colstart;
  1764 + int* Ai = m.m_rowinx;
  1765 + double* Ax = m.m_data;
  1766 +
  1767 + // Allocate result matrix
  1768 + GSparseMatrix C(n, n, Ap[n]);
  1769 +
  1770 + // Allocate and initialise workspace
  1771 + int wrk_size = n;
  1772 + int* wrk_int = new int[wrk_size];
  1773 + if (wrk_int == NULL)
  1774 + throw mem_alloc("cs_symperm(GSparseMatrix*, const int*, int", wrk_size);
  1775 + for (i = 0; i < wrk_size; ++i)
  1776 + wrk_int[i] = 0;
  1777 +
  1778 + // Assign result matrix attributes
  1779 + int* Cp = C.m_colstart;
  1780 + int* Ci = C.m_rowinx;
  1781 + double* Cx = C.m_data;
  1782 +
  1783 + // Count entries in each column of C
  1784 + for (j = 0; j < n; j++) {
  1785 +
  1786 + // Column j of A is column j2 of C
  1787 + j2 = pinv ? pinv[j] : j;
  1788 +
  1789 + // Loop over entries in column j
  1790 + for (p = Ap[j]; p < Ap[j+1]; p++) {
  1791 + i = Ai [p];
  1792 + if (i > j) continue; // skip lower triangular part of A
  1793 + i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C
  1794 + wrk_int[G_MAX(i2, j2)]++; // column count of C
  1795 + }
  1796 + }
  1797 +
  1798 + // Compute column pointers of C
  1799 + cs_cumsum(Cp, wrk_int, n);
  1800 +
  1801 + // Loop over all columns of A
  1802 + for (j = 0 ; j < n ; j++) {
  1803 +
  1804 + // Column j of A is column j2 of C
  1805 + j2 = pinv ? pinv[j] : j;
  1806 +
  1807 + // Loop over entries in column j
  1808 + for (p = Ap[j]; p < Ap[j+1]; p++) {
  1809 + i = Ai [p] ;
  1810 + if (i > j) continue; // skip lower triangular part of A
  1811 + i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C
  1812 + Ci[q = wrk_int[G_MAX(i2,j2)]++] = G_MIN(i2,j2);
  1813 + if (Cx) Cx[q] = Ax[p];
  1814 + }
  1815 + }
  1816 +
  1817 + // Free workspace
  1818 + delete [] wrk_int;
  1819 +
  1820 + // Rectify the number of elements in matrix C
  1821 + C.free_elements(Cp[n], (C.m_elements-Cp[n]));
  1822 +
  1823 + // Return result
  1824 + return C;
  1825 +
  1826 +}
  1827 +
  1828 +
  1829 +/***************************************************************************
  1830 + * cs_transpose *
  1831 + * ----------------------------------------------------------------------- *
  1832 + * Transpose matrix. The flag 'values' allows to avoid copying the actual *
  1833 + * data values. This allows to perform a logical matrix transposition, as *
  1834 + * needed by the symbolic matrix analysis class. *
  1835 + * ----------------------------------------------------------------------- *
  1836 + * NOTE: This routine DOES NOT support pending elements (they have to be *
  1837 + * filled before. *
  1838 + ***************************************************************************/
  1839 +GSparseMatrix cs_transpose(const GSparseMatrix& m, int values)
  1840 +{
  1841 + // Declare and allocate result matrix
  1842 + GSparseMatrix result(m.m_cols, m.m_rows, m.m_elements);
  1843 +
  1844 + // Allocate and initialise workspace
  1845 + int wrk_size = m.m_rows;
  1846 + int* wrk_int = new int[wrk_size];
  1847 + if (wrk_int == NULL)
  1848 + throw mem_alloc("GSparseMatrix::cs_transpose(GSparseMatrix*, int)",
  1849 + wrk_size);
  1850 + for (int i = 0; i < wrk_size; i++)
  1851 + wrk_int[i] = 0;
  1852 +
  1853 + // Setup the number of non-zero elements in each row
  1854 + for (int p = 0; p < m.m_elements; p++)
  1855 + wrk_int[m.m_rowinx[p]]++;
  1856 +
  1857 + // Set row pointers. To use a GSparseSymbolic function we have to
  1858 + // allocate and object (but this does not take memory)
  1859 + cs_cumsum(result.m_colstart, wrk_int, m.m_rows);
  1860 +
  1861 + // Case A: Normal transponse, including assignment of values
  1862 + if (values) {
  1863 + for (int col = 0; col < m.m_cols; col++) {
  1864 + for (int p = m.m_colstart[col] ; p < m.m_colstart[col+1] ; p++) {
  1865 + int i = wrk_int[m.m_rowinx[p]]++;
  1866 + result.m_rowinx[i] = col;
  1867 + result.m_data[i] = m.m_data[p] ;
  1868 + }
  1869 + }
  1870 + }
  1871 +
  1872 + // Case B: Logical transponse, no assignment of values is performed
  1873 + else {
  1874 + for (int col = 0; col < m.m_cols; col++) {
  1875 + for (int p = m.m_colstart[col] ; p < m.m_colstart[col+1] ; p++)
  1876 + result.m_rowinx[wrk_int[m.m_rowinx[p]]++] = col;
  1877 + }
  1878 + }
  1879 +
  1880 + // Return transponse matrix
  1881 + return result;
  1882 +
  1883 +}
  1884 +
  1885 +
1808 1886 /*==========================================================================
1809 1887 = =
1810 1888 = GSparseMatrix exception classes =
... ... @@ -1849,4 +1927,3 @@ double cs_cumsum(int* p, int* c, int n)
1849 1927 // Return cumulative sum of c[0..n-1]
1850 1928 return nz2;
1851 1929 }
1852   -
... ...
src/GSparseMatrix.hpp
  View file @cb738b2
... ... @@ -25,12 +25,9 @@ using namespace std;
25 25 /* __ Definitions ________________________________________________________ */
26 26 #define G_SPARSE_MATRIX_DEFAULT_MEM_BLOCK 512 // Default Memory block size
27 27  
28   -/* __ Enumerators ________________________________________________________ */
29   -
30   -/* __ Structures _________________________________________________________ */
31   -
32   -/* __ Prototypes _________________________________________________________ */
33   -double cs_cumsum(int* p, int* c, int n);
  28 +/* __ Macros (to be moved to GammaLib general header) ____________________ */
  29 +#define G_MIN(a,b) (((a) < (b)) ? (a) : (b))
  30 +#define G_MAX(a,b) (((a) > (b)) ? (a) : (b))
34 31  
35 32  
36 33 /***************************************************************************
... ... @@ -54,6 +51,9 @@ class GSparseMatrix : public GMatrix {
54 51 friend GSparseMatrix fabs(const GSparseMatrix& m);
55 52 friend GSparseMatrix cholesky_decompose(const GSparseMatrix& m, int compress = 1);
56 53 friend GSparseMatrix cholesky_invert(const GSparseMatrix& m, int compress = 1);
  54 + //
  55 + friend GSparseMatrix cs_symperm(const GSparseMatrix& m, const int* pinv);
  56 + friend GSparseMatrix cs_transpose(const GSparseMatrix& m, int values);
57 57  
58 58 public:
59 59 // Constructors and destructors (not inherited)
... ... @@ -113,7 +113,6 @@ private:
113 113 void free_elements(int start, int num);
114 114 void remove_zero_row_col(void);
115 115 void insert_zero_row_col(int rows, int cols);
116   - GSparseMatrix cs_transpose(GSparseMatrix *m, int values);
117 116  
118 117 // Private data members
119 118 int* m_rowinx; // Row-indices of all elements
... ... @@ -292,4 +291,10 @@ GSparseMatrix cholesky_invert(const GSparseMatrix&amp; m, int compress)
292 291 return result;
293 292 }
294 293  
  294 +
  295 +/***************************************************************************
  296 + * Prototypes for other functions *
  297 + ***************************************************************************/
  298 +double cs_cumsum(int* p, int* c, int n);
  299 +
295 300 #endif /* GSPARSEMATRIX_HPP */
... ...
src/GSparseNumeric.cpp
  View file @cb738b2
... ... @@ -78,6 +78,12 @@ GSparseNumeric&amp; GSparseNumeric::operator= (const GSparseNumeric&amp; n)
78 78 // Execute only if object is not identical
79 79 if (this != &n) {
80 80  
  81 + // De-allocate only if memory has indeed been allocated
  82 + if (m_L != NULL) delete m_L;
  83 + if (m_U != NULL) delete m_U;
  84 + if (m_pinv != NULL) delete [] m_pinv;
  85 + if (m_B != NULL) delete [] m_B;
  86 +
81 87 // Initialise private members for clean destruction
82 88 m_L = NULL;
83 89 m_U = NULL;
... ... @@ -87,15 +93,23 @@ GSparseNumeric&amp; GSparseNumeric::operator= (const GSparseNumeric&amp; n)
87 93 m_n_B = 0;
88 94  
89 95 // Copy m_L if it exists
90   - if (n.m_L != NULL)
91   - m_L = new GSparseMatrix(*n.m_L);
  96 + if (n.m_L != NULL) {
  97 + m_L = new GSparseMatrix(*n.m_L);
  98 + if (m_L == NULL)
  99 + throw mem_alloc("GSparseNumeric::operator= (const GSparseNumeric&)",
  100 + 1);
  101 + }
92 102  
93 103 // Copy m_U if it exists
94   - if (n.m_U != NULL)
95   - m_U = new GSparseMatrix(*n.m_U);
  104 + if (n.m_U != NULL) {
  105 + m_U = new GSparseMatrix(*n.m_U);
  106 + if (m_U == NULL)
  107 + throw mem_alloc("GSparseNumeric::operator= (const GSparseNumeric&)",
  108 + 1);
  109 + }
96 110  
97 111 // Copy m_pinv array if it exists
98   - if (n.m_pinv != NULL) {
  112 + if (n.m_pinv != NULL && n.m_n_pinv > 0) {
99 113 m_pinv = new int[n.m_n_pinv];
100 114 if (m_pinv == NULL)
101 115 throw mem_alloc("GSparseNumeric::operator= (const GSparseNumeric&)",
... ... @@ -106,7 +120,7 @@ GSparseNumeric&amp; GSparseNumeric::operator= (const GSparseNumeric&amp; n)
106 120 }
107 121  
108 122 // Copy m_B array if it exists
109   - if (n.m_B != NULL) {
  123 + if (n.m_B != NULL && n.m_n_B > 0) {
110 124 m_B = new double[n.m_n_B];
111 125 if (m_B == NULL)
112 126 throw mem_alloc("GSparseNumeric::operator= (const GSparseNumeric&)",
... ... @@ -138,9 +152,23 @@ GSparseNumeric&amp; GSparseNumeric::operator= (const GSparseNumeric&amp; n)
138 152 * Input: A Sparse matrix *
139 153 * S Symbolic analysis of sparse matrix *
140 154 ***************************************************************************/
141   -void GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix& A,
  155 +void GSparseNumeric::cholesky_numeric_analysis(const GSparseMatrix& A,
142 156 const GSparseSymbolic& S)
143 157 {
  158 + // De-allocate memory that has indeed been previously allocated
  159 + if (m_L != NULL) delete m_L;
  160 + if (m_U != NULL) delete m_U;
  161 + if (m_pinv != NULL) delete [] m_pinv;
  162 + if (m_B != NULL) delete [] m_B;
  163 +
  164 + // Initialise members
  165 + m_L = NULL;
  166 + m_U = NULL;
  167 + m_pinv = NULL;
  168 + m_B = NULL;
  169 + m_n_pinv = 0;
  170 + m_n_B = 0;
  171 +
144 172 // Return if arrays in the symbolic analysis have not been allocated
145 173 if (!S.m_cp || !S.m_parent) return;
146 174  
... ... @@ -173,7 +201,7 @@ void GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix&amp; A,
173 201 int* parent = S.m_parent;
174 202  
175 203 // Assign C = A(p,p) where A and C are symmetric and the upper part stored
176   - GSparseMatrix C = (pinv) ? S.cs_symperm(&A, pinv, 1) : (A);
  204 + GSparseMatrix C = (pinv) ? cs_symperm(A, pinv) : (A);
177 205  
178 206 // Assign workspace pointer
179 207 int* c = wrk_int;
... ... @@ -186,7 +214,7 @@ void GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix&amp; A,
186 214 double* Cx = C.m_data;
187 215  
188 216 // Allocate L matrix
189   - m_L = new GSparseMatrix(n,n,cp[n]);
  217 + m_L = new GSparseMatrix(n, n, cp[n]);
190 218 if (m_L == NULL)
191 219 throw mem_alloc(
192 220 "GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix&, const GSparseSymbolic&)",
... ... @@ -204,7 +232,8 @@ void GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix&amp; A,
204 232 // Compute L(:,k) for L*L' = C
205 233 for (k = 0; k < n; k++) {
206 234  
207   - // Nonzero pattern of L(k,:)
  235 + // Nonzero pattern of L(k,:).
  236 + // Returns -1 if parent = NULL, s = NULL or c = NULL
208 237 top = cs_ereach(&C, k, parent, s, c); // find pattern of L(k,:)
209 238 x[k] = 0; // x (0:k) is now zero
210 239  
... ... @@ -222,14 +251,14 @@ void GSparseNumeric::cholesky_numeric_analysis(GSparseMatrix&amp; A,
222 251  
223 252 // Triangular solve: Solve L(0:k-1,0:k-1) * x = C(:,k)
224 253 for ( ; top < n; top++) {
225   - i = s[top]; // s [top..n-1] is pattern of L(k,:)
226   - lki = x[i]/Lx[Lp[i]]; // L(k,i) = x (i) / L(i,i)
227   - x[i] = 0; // clear x for k+1st iteration
  254 + i = s[top]; // s [top..n-1] is pattern of L(k,:)
  255 + lki = x[i]/Lx[Lp[i]]; // L(k,i) = x (i) / L(i,i)
  256 + x[i] = 0; // clear x for k+1st iteration
228 257 for (p = Lp[i]+1; p < c[i]; p++)
229 258 x[Li[p]] -= Lx[p] * lki;
230   - d -= lki * lki; // d = d - L(k,i)*L(k,i)
231   - p = c [i]++ ;
232   - Li[p] = k; // store L(k,i) in column i
  259 + d -= lki * lki; // d = d - L(k,i)*L(k,i)
  260 + p = c[i]++;
  261 + Li[p] = k; // store L(k,i) in column i
233 262 Lx[p] = lki;
234 263 }
235 264  
... ... @@ -319,7 +348,7 @@ int GSparseNumeric::cs_ereach(const GSparseMatrix* A, int k,
319 348 CS_MARK(w, k);
320 349  
321 350 // Return index top, where s[top..n-1] contains pattern of L(k,:)
322   - return (top);
  351 + return top;
323 352 }
324 353  
325 354  
... ...
src/GSparseNumeric.hpp
  View file @cb738b2
... ... @@ -51,7 +51,7 @@ public:
51 51 GSparseNumeric& operator= (const GSparseNumeric& n);
52 52  
53 53 // Functions
54   - void cholesky_numeric_analysis(GSparseMatrix& m, const GSparseSymbolic& s);
  54 + void cholesky_numeric_analysis(const GSparseMatrix& m, const GSparseSymbolic& s);
55 55  
56 56 private:
57 57 // Functions
... ...
src/GSparseSymbolic.cpp
  View file @cb738b2
... ... @@ -24,6 +24,7 @@ using namespace std;
24 24  
25 25 /* __ Definitions ________________________________________________________ */
26 26 //#define G_DEBUG_SPARSE_SYM_AMD_ORDERING // Debug AMD ordering
  27 +//#define G_DEBUG_SPARSE_CHOLESKY // Analyse Cholesky decomposition
27 28  
28 29  
29 30 /*==========================================================================
... ... @@ -63,11 +64,11 @@ GSparseSymbolic::GSparseSymbolic()
63 64 GSparseSymbolic::~GSparseSymbolic()
64 65 {
65 66 // De-allocate only if memory has indeed been allocated
66   - if (m_pinv != NULL) delete[] m_pinv;
67   - if (m_q != NULL) delete[] m_q;
68   - if (m_parent != NULL) delete[] m_parent;
69   - if (m_cp != NULL) delete[] m_cp;
70   - if (m_leftmost != NULL) delete[] m_leftmost;
  67 + if (m_pinv != NULL) delete [] m_pinv;
  68 + if (m_q != NULL) delete [] m_q;
  69 + if (m_parent != NULL) delete [] m_parent;
  70 + if (m_cp != NULL) delete [] m_cp;
  71 + if (m_leftmost != NULL) delete [] m_leftmost;
71 72  
72 73 // Return
73 74 return;
... ... @@ -88,6 +89,13 @@ GSparseSymbolic&amp; GSparseSymbolic::operator= (const GSparseSymbolic&amp; s)
88 89 // Execute only if object is not identical
89 90 if (this != &s) {
90 91  
  92 + // De-allocate only if memory has indeed been allocated
  93 + if (m_pinv != NULL) delete [] m_pinv;
  94 + if (m_q != NULL) delete [] m_q;
  95 + if (m_parent != NULL) delete [] m_parent;
  96 + if (m_cp != NULL) delete [] m_cp;
  97 + if (m_leftmost != NULL) delete [] m_leftmost;
  98 +
91 99 // Initialise private members for clean destruction
92 100 m_n_pinv = 0;
93 101 m_n_q = 0;
... ... @@ -109,7 +117,7 @@ GSparseSymbolic&amp; GSparseSymbolic::operator= (const GSparseSymbolic&amp; s)
109 117 m_unz = s.m_unz;
110 118  
111 119 // Copy m_pinv array if it exists
112   - if (s.m_pinv != NULL) {
  120 + if (s.m_pinv != NULL && s.m_n_pinv > 0) {
113 121 m_pinv = new int[s.m_n_pinv];
114 122 if (m_pinv == NULL)
115 123 throw mem_alloc("GSparseSymbolic::operator= (const GSparseSymbolic&)",
... ... @@ -120,7 +128,7 @@ GSparseSymbolic&amp; GSparseSymbolic::operator= (const GSparseSymbolic&amp; s)
120 128 }
121 129  
122 130 // Copy m_q array if it exists
123   - if (s.m_q != NULL) {
  131 + if (s.m_q != NULL && s.m_n_q > 0) {
124 132 m_q = new int[s.m_n_q];
125 133 if (m_q == NULL)
126 134 throw mem_alloc("GSparseSymbolic::operator= (const GSparseSymbolic&)",
... ... @@ -131,7 +139,7 @@ GSparseSymbolic&amp; GSparseSymbolic::operator= (const GSparseSymbolic&amp; s)
131 139 }
132 140  
133 141 // Copy m_parent array if it exists
134   - if (s.m_parent != NULL) {
  142 + if (s.m_parent != NULL && s.m_n_parent > 0) {
135 143 m_parent = new int[s.m_n_parent];
136 144 if (m_parent == NULL)
137 145 throw mem_alloc("GSparseSymbolic::operator= (const GSparseSymbolic&)",
... ... @@ -142,7 +150,7 @@ GSparseSymbolic&amp; GSparseSymbolic::operator= (const GSparseSymbolic&amp; s)
142 150 }
143 151  
144 152 // Copy m_cp array if it exists
145   - if (s.m_cp != NULL) {
  153 + if (s.m_cp != NULL && s.m_n_cp > 0) {
146 154 m_cp = new int[s.m_n_cp];
147 155 if (m_cp == NULL)
148 156 throw mem_alloc("GSparseSymbolic::operator= (const GSparseSymbolic&)",
... ... @@ -153,20 +161,20 @@ GSparseSymbolic&amp; GSparseSymbolic::operator= (const GSparseSymbolic&amp; s)
153 161 }
154 162  
155 163 // Copy m_leftmost array if it exists
156   - if (s.m_leftmost != NULL) {
  164 + if (s.m_leftmost != NULL && s.m_n_leftmost > 0) {
157 165 m_leftmost = new int[s.m_n_leftmost];
158 166 if (m_leftmost == NULL)
159 167 throw mem_alloc("GSparseSymbolic::operator= (const GSparseSymbolic&)",
160 168 s.m_n_leftmost);
161 169 for (int i = 0; i < s.m_n_leftmost; ++i)
162 170 m_leftmost[i] = s.m_leftmost[i];
163   - m_leftmost = s.m_leftmost;
  171 + m_n_leftmost = s.m_n_leftmost;
164 172 }
165 173  
166 174 } // endif: object was not identical
167 175  
168 176 // Return
169   - return (*this);
  177 + return *this;
170 178 }
171 179  
172 180  
... ... @@ -186,44 +194,120 @@ GSparseSymbolic&amp; GSparseSymbolic::operator= (const GSparseSymbolic&amp; s)
186 194 * 1: Cholesky decomposition *
187 195 * m Sparse matrix *
188 196 ***************************************************************************/
189   -void GSparseSymbolic::cholesky_symbolic_analysis(int order, GSparseMatrix& m)
  197 +void GSparseSymbolic::cholesky_symbolic_analysis(int order,
  198 + const GSparseMatrix& m)
190 199 {
  200 + // Debug
  201 + #if defined(G_DEBUG_SPARSE_CHOLESKY)
  202 + cout << endl << "GSparseSymbolic::cholesky_symbolic_analysis entered" << endl;
  203 + cout << " order ......................: " << order << endl;
  204 + cout << " number of rows .............: " << m.m_rows << endl;
  205 + cout << " number of columns ..........: " << m.m_cols << endl;
  206 + cout << " number of non-zero elements : " << m.m_elements << endl;
  207 + #endif
  208 +
  209 + // De-allocate memory that has indeed been previously allocated
  210 + if (m_pinv != NULL) delete [] m_pinv;
  211 + if (m_q != NULL) delete [] m_q;
  212 + if (m_parent != NULL) delete [] m_parent;
  213 + if (m_cp != NULL) delete [] m_cp;
  214 + if (m_leftmost != NULL) delete [] m_leftmost;
  215 +
  216 + // Initialise members
  217 + m_n_pinv = 0;
  218 + m_n_q = 0;
  219 + m_n_parent = 0;
  220 + m_n_cp = 0;
  221 + m_n_leftmost = 0;
  222 + m_pinv = NULL;
  223 + m_q = NULL;
  224 + m_parent = NULL;
  225 + m_cp = NULL;
  226 + m_leftmost = NULL;
  227 + m_m2 = 0;
  228 + m_lnz = 0.0;
  229 + m_unz = 0.0;
  230 +
191 231 // Check if order type is valid
192 232 if (order < 0 || order > 1)
193   - throw invalid_order("GSparseSymbolic::cholesky_symbolic_analysis(int, GSparseMatrix*)",
194   - order, 0, 1);
  233 + throw invalid_order(
  234 + "GSparseSymbolic::cholesky_symbolic_analysis(int, GSparseMatrix*)",
  235 + order, 0, 1);
195 236  
196   - // Fill pending element
197   - m.fill_pending();
198   -
199 237 // Assign input matrix attributes
200 238 int n = m.m_cols;
201 239  
202 240 // P = amd(A+A'), or natural
203 241 int* P = cs_amd(order, &m);
  242 + #if defined(G_DEBUG_SPARSE_CHOLESKY)
  243 + cout << " AMD ordering permutation (" << P << ")" << endl;
  244 + if (P != NULL) {
  245 + for (int i = 0; i < n; ++i)
  246 + cout << " " << P[i];
  247 + cout << endl;
  248 + }
  249 + #endif
204 250  
205   - // Find inverse permutation and store it in class member 'm_pinv'
  251 + // Find inverse permutation and store it in class member 'm_pinv'.
  252 + // Note that if P = NULL or n < 1 this function returns NULL.
206 253 m_pinv = cs_pinv(P, n);
207 254 m_n_pinv = n;
  255 + #if defined(G_DEBUG_SPARSE_CHOLESKY)
  256 + cout << " Inverse permutation (" << m_pinv << ")" << endl;
  257 + if (m_pinv != NULL) {
  258 + for (int i = 0; i < n; ++i)
  259 + cout << " " << m_pinv[i];
  260 + cout << endl;
  261 + }
  262 + #endif
208 263  
209 264 // Delete workspace
210   - delete [] P;
  265 + if (P != NULL) delete [] P;
211 266  
212 267 // C = spones(triu(A(P,P)))
213   - GSparseMatrix C = cs_symperm(&m, m_pinv, 0);
  268 + GSparseMatrix C = cs_symperm(m, m_pinv);
  269 + #if defined(G_DEBUG_SPARSE_CHOLESKY)
  270 + cout << " C = spones(triu(A(P,P))) " << C << endl;
  271 + #endif
214 272  
215 273 // Find etree of C and store it in class member 'm_parent'
216 274 m_parent = cs_etree(&C, 0);
217 275 m_n_parent = n;
  276 + #if defined(G_DEBUG_SPARSE_CHOLESKY)
  277 + cout << " Elimination tree (" << m_parent << ")" << endl;
  278 + if (m_parent != NULL) {
  279 + for (int i = 0; i < n; ++i)
  280 + cout << " " << m_parent[i];
  281 + cout << endl;
  282 + }
  283 + #endif
218 284  
219 285 // Post order the etree
  286 + // Note that if m_parent = NULL or n < 1 this function returns NULL.
220 287 int* post = cs_post(m_parent, n);
  288 + #if defined(G_DEBUG_SPARSE_CHOLESKY)
  289 + cout << " Post ordered elimination tree (" << post << ")" << endl;
  290 + if (post != NULL) {
  291 + for (int i = 0; i < n; ++i)
  292 + cout << " " << post[i];
  293 + cout << endl;
  294 + }
  295 + #endif
221 296  
222 297 // Find column counts of chol(C)
  298 + // Note that if m_parent = NULL or post = NULL this function returns NULL.
223 299 int* c = cs_counts(&C, m_parent, post, 0);
  300 + #if defined(G_DEBUG_SPARSE_CHOLESKY)
  301 + cout << " Column counts (" << c << ")" << endl;
  302 + if (c != NULL) {
  303 + for (int i = 0; i < n; ++i)
  304 + cout << " " << c[i];
  305 + cout << endl;
  306 + }
  307 + #endif
224 308  
225 309 // Delete workspace
226   - delete [] post;
  310 + if (post != NULL) delete [] post;
227 311  
228 312 // Allocate column pointers for Cholesky decomposition
229 313 m_n_cp = n+1;
... ... @@ -233,11 +317,52 @@ void GSparseSymbolic::cholesky_symbolic_analysis(int order, GSparseMatrix&amp; m)
233 317  
234 318 // Find column pointers for L
235 319 m_unz = m_lnz = cs_cumsum(m_cp, c, n);
  320 + #if defined(G_DEBUG_SPARSE_CHOLESKY)
  321 + cout << " Column pointers for L (" << m_cp << ")" << endl;
  322 + if (m_cp != NULL) {
  323 + for (int i = 0; i < m_n_cp; ++i)
  324 + cout << " " << m_cp[i];
  325 + cout << endl;
  326 + }
  327 + cout << " Number of non-zero elements in L: " << m_lnz << endl;
  328 + #endif
236 329  
237 330 // Delete workspace
238   - delete [] c;
  331 + if (c != NULL) delete [] c;
  332 +
  333 + // Delete symbolic analysis if it is not valid
  334 + if (m_lnz < 0) {
  335 +
  336 + // De-allocate memory that has indeed been previously allocated
  337 + if (m_pinv != NULL) delete [] m_pinv;
  338 + if (m_q != NULL) delete [] m_q;
  339 + if (m_parent != NULL) delete [] m_parent;
  340 + if (m_cp != NULL) delete [] m_cp;
  341 + if (m_leftmost != NULL) delete [] m_leftmost;
239 342  
240   -// return ((S->lnz >= 0) ? S : cs_sfree (S)) ;
  343 + // Initialise members
  344 + m_n_pinv = 0;
  345 + m_n_q = 0;
  346 + m_n_parent = 0;
  347 + m_n_cp = 0;
  348 + m_n_leftmost = 0;
  349 + m_pinv = NULL;
  350 + m_q = NULL;
  351 + m_parent = NULL;
  352 + m_cp = NULL;
  353 + m_leftmost = NULL;
  354 + m_m2 = 0;
  355 + m_lnz = 0.0;
  356 + m_unz = 0.0;
  357 + }
  358 +
  359 + // Debug
  360 + #if defined(G_DEBUG_SPARSE_CHOLESKY)
  361 + cout << "GSparseSymbolic::cholesky_symbolic_analysis finished" << endl;
  362 + #endif
  363 +
  364 + // Return
  365 + return;
241 366 }
242 367  
243 368  
... ... @@ -261,11 +386,11 @@ void GSparseSymbolic::cholesky_symbolic_analysis(int order, GSparseMatrix&amp; m)
261 386 * A Sparse matrix *
262 387 * Output: p[0..n] integer array of n+1 elements *
263 388 ***************************************************************************/
264   -int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
  389 +int* GSparseSymbolic::cs_amd(int order, const GSparseMatrix* A)
265 390 {
266 391 // Throw exception if order is not a decomposition
267 392 if (order < 0 || order > 3)
268   - throw invalid_order("GSparseSymbolic::cs_amd(int, GSparseMatrix*)",
  393 + throw invalid_order("GSparseSymbolic::cs_amd(int, const GSparseMatrix*)",
269 394 order, 0, 3);
270 395  
271 396 // Debug
... ... @@ -305,17 +430,20 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
305 430 GSparseMatrix C(m,n);
306 431  
307 432 // Get (logical) transpose of A: AT = A'
308   - GSparseMatrix AT = A->cs_transpose(A,0);
  433 +// NOTE: WE ONLY NEED THE LOGICAL TRANSPOSE HERE. HOWEVER, WE HAVE NO
  434 +// LOGICAL ADDITION SO FAR ...
  435 +// GSparseMatrix AT = cs_transpose(*A, 0);
  436 + GSparseMatrix AT = cs_transpose(*A, 1);
309 437  
310 438 // Find dense threshold
311 439 dense = (int)CS_MAX(16, 10 * sqrt((double)n));
312 440 dense = CS_MIN(n-2, dense);
313 441  
314 442 // Case A: Cholesky decomposition of a symmetric matrix: C = A + A'
315   - if (order == 1 && n == m) {
  443 + if (order == 1 && n == m)
  444 +// NOTE: HERE WE ONLY NEED LOGICAL ADDITION
316 445 // C = cs_add (A, AT, 0, 0);
317 446 C = *A + AT;
318   - }
319 447  
320 448 // Case B: LU decomposition: C = A' * A with no dense rows
321 449 else if (order == 2) {
... ... @@ -337,24 +465,35 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
337 465 ATp[m] = p2;
338 466  
339 467 // Get (logical) transpose of AT: A2 = AT'
340   - GSparseMatrix A2 = AT.cs_transpose(&AT,0);
  468 +// NOTE: WE ONLY NEED THE LOGICAL TRANSPOSE HERE. HOWEVER, WE HAVE NO
  469 +// LOGICAL MULTIPLICATION SO FAR ...
  470 +// GSparseMatrix A2 = cs_transpose(AT, 0);
  471 + GSparseMatrix A2 = cs_transpose(AT, 1);
341 472  
342 473 // C = A' * A with no dense rows
  474 +// NOTE: cs_multiply NOT YET IMPLEMENTED
343 475 // C = A2 ? cs_multiply(AT, A2) : NULL;
344 476 C = AT * A2;
345 477 }
346 478  
347 479 // Case C: QR decomposition: C = A' * A
348   - else {
  480 + else
  481 +// NOTE: cs_multiply NOT YET IMPLEMENTED
349 482 // C = cs_multiply(AT, A);
350 483 C = AT * *A;
351   - }
352   -
353   - // Stop in case of an error
354   -// if (!C) return (NULL) ;
355 484  
  485 + // Debug
  486 + #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
  487 + cout << " Matrix C " << C << endl;
  488 + #endif
  489 +
356 490 // Drop diagonal entries from C
357   - cs_fkeep(&C, &cs_diag, NULL);
  491 + int cnz_diag = cs_fkeep(&C, &cs_diag, NULL);
  492 +
  493 + // Debug
  494 + #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
  495 + cout << " Dropped diagonal entries from matrix C " << C << endl;
  496 + #endif
358 497  
359 498 // Assign C matrix attributes
360 499 int* Cp = C.m_colstart;
... ... @@ -372,10 +511,19 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
372 511 throw mem_alloc("GSparseSymbolic::cs_amd (int, GSparseMatrix*)", wrk_size);
373 512  
374 513 // Add elbow room to C
375   -// t = cnz + cnz/5 + 2*n;
376   -// if (!P || !W || !cs_sprealloc (C, t)) return (cs_idone (P, C, W, 0)) ;
377   - int t = cnz/5 + 2*n;
378   - C.alloc_elements(cnz, t);
  514 + int elbow_room = cnz/5 + 2*n; // Request additional # of elements
  515 + C.alloc_elements(cnz, elbow_room); // Appand elements
  516 + int nzmax = C.m_elements; // Save the maximum # of entries for garbage coll.
  517 + C.m_elements -= elbow_room; // Reverse element increase of alloc_elements
  518 +
  519 + // (Re-)assign C matrix attributes
  520 + Cp = C.m_colstart;
  521 + int* Ci = C.m_rowinx;
  522 +
  523 + // Debug
  524 + #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
  525 + cout << " Added elbow room to matrix C " << C << endl;
  526 + #endif
379 527  
380 528 // Set workspace pointers
381 529 int* len = wrk_int;
... ... @@ -386,7 +534,7 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
386 534 int* degree = wrk_int + 5*(n+1);
387 535 int* w = wrk_int + 6*(n+1);
388 536 int* hhead = wrk_int + 7*(n+1);
389   - int* last = P; // use P as workspace for last
  537 + int* last = P; // use P as workspace for last
390 538  
391 539  
392 540 // ======================================================================
... ... @@ -398,35 +546,38 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
398 546 cout << " Step 2: Initialize quotient graph" << endl;
399 547 #endif
400 548  
401   - // (Re-)assign C matrix attributes
402   - int nzmax = C.m_elements; // Maximum number of entries
403   - Cp = C.m_colstart;
404   - int* Ci = C.m_rowinx;
405   -
406 549 // Setup array that contains for each column the number of non-zero elements
407 550 for (k = 0 ; k < n ; k++)
408 551 len[k] = Cp[k+1] - Cp[k];
409 552 len[n] = 0;
  553 +
  554 + // Debug
  555 + #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
  556 + cout << " Non-zero elements per column: ";
  557 + for (k = 0 ; k <= n ; k++)
  558 + cout << " " << len[k];
  559 + cout << endl;
  560 + #endif
410 561  
411   - // Loop over all columns
  562 + // Loop over all nodes (i.e. matrix columns)
412 563 for (i = 0 ; i <= n ; i++) {
413   - head[i] = -1; // degree list i is empty
  564 + head[i] = -1; // degree list i is empty
414 565 last[i] = -1;
415 566 next[i] = -1;
416   - hhead[i] = -1; // hash list i is empty
417   - nv[i] = 1; // node i is just one node
418   - w[i] = 1; // node i is alive
419   - elen[i] = 0; // Ek of node i is empty
420   - degree[i] = len[i]; // degree of node i
  567 + hhead[i] = -1; // hash list i is empty
  568 + nv[i] = 1; // node i is just one node
  569 + w[i] = 1; // node i is alive
  570 + elen[i] = 0; // Ek of node i is empty
  571 + degree[i] = len[i]; // degree of node i
421 572 }
422 573  
423 574 // Clear w
424 575 mark = cs_wclear(0, 0, w, n);
425 576  
426 577 // End-point initialisations
427   - elen[n] = -2; // n is a dead element
428   - Cp[n] = -1; // n is a root of assembly tree
429   - w[n] = 0; // n is a dead element
  578 + elen[n] = -2; // n is a dead element
  579 + Cp[n] = -1; // n is a root of assembly tree
  580 + w[n] = 0; // n is a dead element
430 581  
431 582  
432 583 // ======================================================================
... ... @@ -438,7 +589,7 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
438 589 cout << " Step 3: Initialize degree lists" << endl;
439 590 #endif
440 591  
441   - // Loop over all columns
  592 + // Loop over all nodes
442 593 for (i = 0 ; i < n ; i++) {
443 594  
444 595 // Get degree of node
... ... @@ -446,25 +597,25 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
446 597  
447 598 // Case A: Node i is empty
448 599 if (d == 0) {
449   - elen[i] = -2; // element i is dead
  600 + elen[i] = -2; // element i is dead
450 601 nel++ ;
451   - Cp[i] = -1; // i is a root of assembly tree
452   - w[i] = 0; // i is a dead element
  602 + Cp[i] = -1; // i is a root of assembly tree
  603 + w[i] = 0; // i is a dead element
453 604 }
454 605  
455 606 // Case B: Node is dense
456 607 else if (d > dense) {
457   - nv[i] = 0; // absorb i into element n
458   - elen[i] = -1; // node i is dead
459   - nel++ ;
460   - Cp[i] = CS_FLIP(n);
461   - nv[n]++ ;
  608 + nv[i] = 0; // absorb i into element n
  609 + elen[i] = -1; // node i is dead
  610 + nel++;
  611 + Cp[i] = CS_FLIP(n);
  612 + nv[n]++;
462 613 }
463 614  
464 615 // Case C: Node is sparse
465 616 else {
466 617 if (head[d] != -1) last[head[d]] = i;
467   - next[i] = head[d]; // put node i in degree list d
  618 + next[i] = head[d]; // put node i in degree list d
468 619 head[d] = i;
469 620 }
470 621 } // endfor: looped over all columns
... ... @@ -485,65 +636,65 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
485 636 // Select node of minimum approximate degree
486 637 for (k = -1; mindeg < n && (k = head [mindeg]) == -1; mindeg++) ;
487 638 if (next[k] != -1) last [next[k]] = -1;
488   - head[mindeg] = next[k]; // remove k from degree list
489   - elenk = elen[k]; // elenk = |Ek|
490   - nvk = nv[k]; // # of nodes k represents
491   - nel += nvk; // nv[k] nodes of A eliminated
  639 + head[mindeg] = next[k]; // remove k from degree list
  640 + elenk = elen[k]; // elenk = |Ek|
  641 + nvk = nv[k]; // # of nodes k represents
  642 + nel += nvk; // nv[k] nodes of A eliminated
492 643  
493   - // Garbage collection
  644 + // Garbage collection. This frees the elbow rooms in matrix C
494 645 if (elenk > 0 && cnz + mindeg >= nzmax) {
495 646  
496 647 // Loop over all columns
497   - for (j = 0 ; j < n ; j++) {
  648 + for (j = 0; j < n; j++) {
498 649 // If j is a live node or element, then ...
499 650 if ((p = Cp[j]) >= 0) {
500   - Cp[j] = Ci[p]; // save first entry of object
501   - Ci[p] = CS_FLIP(j); // first entry is now CS_FLIP(j)
  651 + Cp[j] = Ci[p]; // save first entry of object
  652 + Ci[p] = CS_FLIP(j); // first entry is now CS_FLIP(j)
502 653 }
503 654 }
504 655  
505 656 // Scan all memory
506   - for (q = 0, p = 0; p < cnz ; ) {
  657 + for (q = 0, p = 0; p < cnz; ) {
507 658 // If we found object j, then ...
508 659 if ((j = CS_FLIP(Ci[p++])) >= 0) {
509   - Ci[q] = Cp[j]; // restore first entry of object
510   - Cp[j] = q++; // new pointer to object j
511   - for (k3 = 0; k3 < len [j]-1; k3++)
  660 + Ci[q] = Cp[j]; // restore first entry of object
  661 + Cp[j] = q++; // new pointer to object j
  662 + for (k3 = 0; k3 < len[j]-1; k3++)
512 663 Ci[q++] = Ci[p++];
513 664 }
514 665 }
515   - cnz = q; // Ci[cnz...nzmax-1] now free
  666 + cnz = q; // Ci[cnz...nzmax-1] now free
516 667  
517 668 } // endif: garbage collection done
518 669  
519 670 // Construct new element ---------------------------------------- */
520 671 dk = 0;
521   - nv[k] = -nvk; // flag k as in Lk
  672 + nv[k] = -nvk; // flag k as in Lk
522 673 p = Cp[k];
523   - pk1 = (elenk == 0) ? p : cnz; // do in place if elen[k] == 0
  674 + pk1 = (elenk == 0) ? p : cnz; // do in place if elen[k] == 0
524 675 pk2 = pk1;
525 676 for (k1 = 1; k1 <= elenk + 1; k1++) {
526 677  
527 678 // ...
528 679 if (k1 > elenk) {
529   - e = k; // search the nodes in k
530   - pj = p; // list of nodes starts at Ci[pj]
531   - ln = len [k] - elenk; // length of list of nodes in k
  680 + e = k; // search the nodes in k
  681 + pj = p; // list of nodes starts at Ci[pj]
  682 + ln = len[k] - elenk; // length of list of nodes in k
532 683 }
533 684 else {
534   - e = Ci[p++]; // search the nodes in e
  685 + e = Ci[p++]; // search the nodes in e
535 686 pj = Cp[e];
536   - ln = len[e]; // length of list of nodes in e
  687 + ln = len[e]; // length of list of nodes in e
537 688 }
538 689  
539 690 // ...
540 691 for (k2 = 1 ; k2 <= ln ; k2++) {
541 692 i = Ci[pj++];
542   - if ((nvi = nv [i]) <= 0) // Skip of node i is dead or seen
  693 + if ((nvi = nv [i]) <= 0) // Skip of node i is dead or seen
543 694 continue;
544   - dk += nvi; // degree[Lk] += size of node i
545   - nv[i] = -nvi; // negate nv[i] to denote i in Lk
546   - Ci[pk2++] = i; // place i in Lk
  695 + dk += nvi; // degree[Lk] += size of node i
  696 + nv[i] = -nvi; // negate nv[i] to denote i in Lk
  697 + Ci[pk2++] = i; // place i in Lk
547 698 if (next[i] != -1)
548 699 last[next[i]] = last[i];
549 700  
... ... @@ -556,16 +707,16 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
556 707  
557 708 // ...
558 709 if (e != k) {
559   - Cp[e] = CS_FLIP (k); // absorb e into k
560   - w[e] = 0; // e is now a dead element
  710 + Cp[e] = CS_FLIP(k); // absorb e into k
  711 + w[e] = 0; // e is now a dead element
561 712 }
562 713 } // endfor: k1
563 714  
564   - if (elenk != 0) cnz = pk2; // Ci [cnz...nzmax] is free
565   - degree[k] = dk; // external degree of k - |Lk\i|
566   - Cp[k] = pk1; // element k is in Ci[pk1..pk2-1]
  715 + if (elenk != 0) cnz = pk2; // Ci [cnz...nzmax] is free
  716 + degree[k] = dk; // external degree of k - |Lk\i|
  717 + Cp[k] = pk1; // element k is in Ci[pk1..pk2-1]
567 718 len[k] = pk2 - pk1;
568   - elen[k] = -2; // k is now an element
  719 + elen[k] = -2; // k is now an element
569 720  
570 721 // Clear w if necessary
571 722 mark = cs_wclear(mark, lemax, w, n);
... ... @@ -573,85 +724,86 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
573 724 // Scan 1: Find set differences |Le\Lk|
574 725 for (pk = pk1 ; pk < pk2 ; pk++) {
575 726 i = Ci [pk];
576   - if ((eln = elen [i]) <= 0) // skip if elen[i] empty
  727 + if ((eln = elen[i]) <= 0) // skip if elen[i] empty
577 728 continue;
578   - nvi = -nv[i]; // nv [i] was negated
  729 + nvi = -nv[i]; // nv [i] was negated
579 730 wnvi = mark - nvi;
580 731  
581 732 // Scan Ei
582 733 for (p = Cp[i]; p <= Cp[i] + eln - 1; p++) {
583   - e = Ci [p] ;
  734 + e = Ci[p];
584 735 if (w[e] >= mark)
585   - w[e] -= nvi; // decrement |Le\Lk|
586   - else if (w[e] != 0) // ensure e is a live element
587   - w[e] = degree[e] + wnvi; // 1st time e seen in scan 1
  736 + w[e] -= nvi; // decrement |Le\Lk|
  737 + else if (w[e] != 0) // ensure e is a live element
  738 + w[e] = degree[e] + wnvi; // 1st time e seen in scan 1
588 739 }
589 740 } // endfor: scan 1 finished
590 741  
591 742 // Scan 2: Degree update
592 743 for (pk = pk1; pk < pk2; pk++) {
593   - i = Ci [pk] ; /* consider node i in Lk */
594   - p1 = Cp [i] ;
595   - p2 = p1 + elen [i] - 1 ;
596   - pn = p1 ;
  744 + i = Ci[pk]; // consider node i in Lk
  745 + p1 = Cp[i];
  746 + p2 = p1 + elen [i] - 1;
  747 + pn = p1;
597 748  
598 749 // Scan Ei
599 750 for (h = 0, d = 0, p = p1 ; p <= p2 ; p++) {
600   - e = Ci [p] ;
  751 + e = Ci[p];
601 752 // If e is an unabsorbed element, then ...
602   - if (w [e] != 0) {
603   - dext = w [e] - mark ; /* dext = |Le\Lk| */
  753 + if (w[e] != 0) {
  754 + dext = w[e] - mark; // dext = |Le\Lk|
604 755 if (dext > 0) {
605   - d += dext ; /* sum up the set differences */
606   - Ci [pn++] = e ; /* keep e in Ei */
607   - h += e ; /* compute the hash of node i */
  756 + d += dext; // sum up the set differences
  757 + Ci[pn++] = e; // keep e in Ei
  758 + h += e; // compute the hash of node i
608 759 }
609 760 else {
610   - Cp [e] = CS_FLIP (k) ; /* aggressive absorb. e->k */
611   - w [e] = 0 ; /* e is a dead element */
  761 + Cp[e] = CS_FLIP(k); // aggressive absorb. e->k
  762 + w[e] = 0; // e is a dead element
612 763 }
613 764 }
614 765 } // endfor: scanned Ei
615 766  
616 767 // elen[i] = |Ei|
617   - elen [i] = pn - p1 + 1 ;
  768 + elen[i] = pn - p1 + 1;
618 769 p3 = pn ;
619   - p4 = p1 + len [i] ;
  770 + p4 = p1 + len[i];
620 771  
621 772 // Prune edges in Ai
622 773 for (p = p2 + 1 ; p < p4 ; p++) {
623   - j = Ci [p] ;
624   - if ((nvj = nv [j]) <= 0) continue ; /* node j dead or in Lk */
625   - d += nvj ; /* degree(i) += |j| */
626   - Ci [pn++] = j ; /* place j in node list of i */
627   - h += j ; /* compute hash for node i */
  774 + j = Ci[p];
  775 + if ((nvj = nv[j]) <= 0) // node j dead or in Lk
  776 + continue;
  777 + d += nvj; // degree(i) += |j|
  778 + Ci[pn++] = j; // place j in node list of i
  779 + h += j; // compute hash for node i
628 780 }
629 781  
630 782 // Check for mass elimination
631 783 if (d == 0) {
632   - Cp [i] = CS_FLIP (k) ; /* absorb i into k */
633   - nvi = -nv [i] ;
634   - dk -= nvi ; /* |Lk| -= |i| */
635   - nvk += nvi ; /* |k| += nv[i] */
636   - nel += nvi ;
637   - nv [i] = 0 ;
638   - elen [i] = -1 ; /* node i is dead */
  784 + Cp[i] = CS_FLIP(k); // absorb i into k
  785 + nvi = -nv[i];
  786 + dk -= nvi; // |Lk| -= |i|
  787 + nvk += nvi; // |k| += nv[i]
  788 + nel += nvi;
  789 + nv[i] = 0;
  790 + elen[i] = -1; // node i is dead
639 791 }
640 792 else {
641   - degree [i] = CS_MIN (degree [i], d) ; /* update degree(i) */
642   - Ci [pn] = Ci [p3] ; /* move first node to end */
643   - Ci [p3] = Ci [p1] ; /* move 1st el. to end of Ei */
644   - Ci [p1] = k ; /* add k as 1st element in of Ei */
645   - len [i] = pn - p1 + 1 ; /* new len of adj. list of node i */
646   - h %= n ; /* finalize hash of i */
647   - next [i] = hhead [h] ; /* place i in hash bucket */
648   - hhead [h] = i ;
649   - last [i] = h ; /* save hash of i in last[i] */
  793 + degree[i] = CS_MIN(degree[i], d); // update degree(i)
  794 + Ci[pn] = Ci[p3]; // move first node to end
  795 + Ci[p3] = Ci[p1]; // move 1st el. to end of Ei
  796 + Ci[p1] = k; // add k as 1st element in of Ei
  797 + len[i] = pn - p1 + 1; // new len of adj. list of node i
  798 + h %= n; // finalize hash of i
  799 + next[i] = hhead[h]; // place i in hash bucket
  800 + hhead[h] = i;
  801 + last[i] = h; // save hash of i in last[i]
650 802 }
651 803 } // endfor: scan 2 is done
652 804  
653 805 // finalize |Lk|
654   - degree[k] = dk ;
  806 + degree[k] = dk;
655 807 lemax = CS_MAX(lemax, dk);
656 808  
657 809 // Clear w
... ... @@ -659,33 +811,35 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
659 811  
660 812 // Supernode detection
661 813 for (pk = pk1 ; pk < pk2 ; pk++) {
662   - i = Ci [pk] ;
663   - if (nv [i] >= 0) continue ; /* skip if i is dead */
664   - h = last [i] ; /* scan hash bucket of node i */
665   - i = hhead [h] ;
666   - hhead [h] = -1 ; /* hash bucket will be empty */
  814 + i = Ci[pk];
  815 + if (nv[i] >= 0) continue; // skip if i is dead
  816 + h = last[i]; // scan hash bucket of node i
  817 + i = hhead[h];
  818 + hhead[h] = -1; // hash bucket will be empty
667 819  
668 820 // ...
669   - for ( ; i != -1 && next [i] != -1 ; i = next [i], mark++) {
670   - ln = len [i] ;
671   - eln = elen [i] ;
672   - for (p = Cp [i]+1 ; p <= Cp [i] + ln-1 ; p++) w [Ci [p]] = mark;
673   - jlast = i ;
  821 + for ( ; i != -1 && next[i] != -1 ; i = next[i], mark++) {
  822 + ln = len[i];
  823 + eln = elen[i];
  824 + for (p = Cp[i]+1 ; p <= Cp[i] + ln-1 ; p++)
  825 + w[Ci[p]] = mark;
  826 + jlast = i;
674 827  
675 828 // Compare i with all j
676   - for (j = next [i] ; j != -1 ; ) {
677   - ok = (len [j] == ln) && (elen [j] == eln) ;
678   - for (p = Cp [j] + 1 ; ok && p <= Cp [j] + ln - 1 ; p++) {
679   - if (w [Ci [p]] != mark) ok = 0 ; /* compare i and j*/
  829 + for (j = next[i]; j != -1; ) {
  830 + ok = (len[j] == ln) && (elen[j] == eln) ;
  831 + for (p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++) {
  832 + if (w[Ci[p]] != mark) // compare i and j
  833 + ok = 0;
680 834 }
681 835  
682 836 // i and j are identical
683 837 if (ok) {
684   - Cp[j] = CS_FLIP(i); /* absorb j into i */
  838 + Cp[j] = CS_FLIP(i); // absorb j into i
685 839 nv[i] += nv[j];
686 840 nv[j] = 0;
687   - elen[j] = -1; /* node j is dead */
688   - j = next[j]; /* delete j from hash bucket */
  841 + elen[j] = -1; // node j is dead
  842 + j = next[j]; // delete j from hash bucket
689 843 next[jlast] = j;
690 844 }
691 845 // j and i are different
... ... @@ -697,57 +851,75 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
697 851 } //
698 852 } // endfor: supernode detection
699 853  
700   - // Finalize new element
701   - for (p = pk1, pk = pk1 ; pk < pk2 ; pk++) /* finalize Lk */
702   - {
703   - i = Ci [pk] ;
704   - if ((nvi = -nv [i]) <= 0) continue ;/* skip if i is dead */
705   - nv [i] = nvi ; /* restore nv[i] */
706   - d = degree [i] + dk - nvi ; /* compute external degree(i) */
707   - d = CS_MIN (d, n - nel - nvi) ;
708   - if (head [d] != -1) last [head [d]] = i ;
709   - next [i] = head [d] ; /* put i back in degree list */
710   - last [i] = -1 ;
711   - head [d] = i ;
712   - mindeg = CS_MIN (mindeg, d) ; /* find new minimum degree */
713   - degree [i] = d ;
714   - Ci [p++] = i ; /* place i in Lk */
  854 + // Finalize new element Lk
  855 + for (p = pk1, pk = pk1; pk < pk2; pk++) {
  856 + i = Ci[pk];
  857 + if ((nvi = -nv[i]) <= 0) // skip if i is dead
  858 + continue;
  859 + nv[i] = nvi; // restore nv[i]
  860 + d = degree[i] + dk - nvi; // compute external degree(i)
  861 + d = CS_MIN(d, n - nel - nvi);
  862 + if (head[d] != -1) last[head[d]] = i;
  863 + next [i] = head [d]; // put i back in degree list
  864 + last [i] = -1 ;
  865 + head [d] = i ;
  866 + mindeg = CS_MIN (mindeg, d); // find new minimum degree
  867 + degree [i] = d ;
  868 + Ci [p++] = i ; // place i in Lk
715 869 }
716   - nv [k] = nvk ; /* # nodes absorbed into k */
717   - if ((len [k] = p-pk1) == 0) /* length of adj list of element k*/
718   - {
719   - Cp [k] = -1 ; /* k is a root of the tree */
720   - w [k] = 0 ; /* k is now a dead element */
  870 +
  871 + // Number of nodes absorbed into k
  872 + nv[k] = nvk;
  873 +
  874 + // ...
  875 + if ((len [k] = p-pk1) == 0) { // length of adj list of element k
  876 + Cp [k] = -1 ; // k is a root of the tree
  877 + w [k] = 0 ; // k is now a dead element
721 878 }
722   - if (elenk != 0) cnz = p ; /* free unused space in Lk */
723   - }
  879 + if (elenk != 0) cnz = p ; // free unused space in Lk
  880 +
  881 + } // endwhile
  882 +
  883 +
  884 + // ======================================================================
  885 + // Step 5: Postordering
  886 + // ======================================================================
  887 +
  888 + // Debug
  889 + #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
  890 + cout << " Step 5: Postordering" << endl;
  891 + #endif
724 892  
725   - // Postordering
726   - for (i = 0 ; i < n ; i++) Cp [i] = CS_FLIP (Cp [i]) ;/* fix assembly tree */
727   - for (j = 0 ; j <= n ; j++) head [j] = -1 ;
728   - for (j = n ; j >= 0 ; j--) /* place unordered nodes in lists */
729   - {
730   - if (nv [j] > 0) continue ; /* skip if j is an element */
731   - next [j] = head [Cp [j]] ; /* place j in list of its parent */
732   - head [Cp [j]] = j ;
733   - }
734   - for (e = n ; e >= 0 ; e--) /* place elements in lists */
735   - {
736   - if (nv [e] <= 0) continue ; /* skip unless e is an element */
737   - if (Cp [e] != -1)
738   - {
739   - next [e] = head [Cp [e]] ; /* place e in list of its parent */
740   - head [Cp [e]] = e ;
  893 + // Fix assembly tree
  894 + for (i = 0; i < n; i++)
  895 + Cp[i] = CS_FLIP(Cp[i]);
  896 + for (j = 0 ; j <= n ; j++)
  897 + head[j] = -1;
  898 +
  899 + // Place unordered nodes in lists
  900 + for (j = n ; j >= 0 ; j--) {
  901 + if (nv[j] > 0) continue; // skip if j is an element
  902 + next[j] = head [Cp[j]]; // place j in list of its parent
  903 + head[Cp[j]] = j;
  904 + }
  905 +
  906 + // Place elements in lists
  907 + for (e = n ; e >= 0 ; e--) {
  908 + if (nv[e] <= 0) continue; // skip unless e is an element
  909 + if (Cp[e] != -1) {
  910 + next[e] = head[Cp[e]]; // place e in list of its parent
  911 + head[Cp[e]] = e;
741 912 }
742 913 }
743 914  
  915 +
744 916 // ======================================================================
745   - // Step 5: Postorder the assembly tree
  917 + // Step 6: Postorder the assembly tree
746 918 // ======================================================================
747 919  
748 920 // Debug
749 921 #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
750   - cout << " Step 5: Postorder the assembly tree" << endl;
  922 + cout << " Step 6: Postorder the assembly tree" << endl;
751 923 #endif
752 924  
753 925 // Postorder the assembly tree
... ... @@ -769,7 +941,7 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
769 941 #endif
770 942  
771 943 // Return result
772   - return (P);
  944 + return P;
773 945 }
774 946  
775 947  
... ... @@ -784,7 +956,9 @@ int* GSparseSymbolic::cs_amd(int order, GSparseMatrix* A)
784 956 * ata 0: count LL'=A, 1: count LL'=A'A *
785 957 * Output: p[0..n] integer array of n+1 elements *
786 958 ***************************************************************************/
787   -int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent,
  959 +#define HEAD(k,j) (ata ? head [k] : j)
  960 +#define NEXT(J) (ata ? next [J] : -1)
  961 +int* GSparseSymbolic::cs_counts(const GSparseMatrix* A, const int* parent,
788 962 const int* post, int ata)
789 963 {
790 964 // Declare loop variables
... ... @@ -792,7 +966,7 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent,
792 966  
793 967 // Return NULL pointer if input arrays are NULL or sparse matrix is
794 968 // empty
795   - if (!parent || !post) return (NULL);
  969 + if (!parent || !post) return NULL;
796 970  
797 971 // Assign input matrix attributes
798 972 int m = A->m_rows;
... ... @@ -801,8 +975,9 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent,
801 975 // Allocate result
802 976 int* colcount = new int[n];
803 977 if (colcount == NULL)
804   - throw mem_alloc("GSparseSymbolic::cs_counts(GSparseMatrix*, const int*, const int*, int)",
805   - n);
  978 + throw mem_alloc(
  979 + "GSparseSymbolic::cs_counts(GSparseMatrix*, const int*, const int*, int)",
  980 + n);
806 981 int* delta = colcount;
807 982  
808 983 // Allocate workspace
... ... @@ -812,8 +987,8 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent,
812 987 throw mem_alloc("GSparseSymbolic::cs_counts(GSparseMatrix*, const int*, const int*, int)",
813 988 wrk_size);
814 989  
815   - // Get transpose of A: AT = A'
816   - GSparseMatrix AT = A->cs_transpose(A,0);
  990 + // Get (logical) transpose of A: AT = A'
  991 + GSparseMatrix AT = cs_transpose(*A, 0);
817 992  
818 993 // Set-up pointers to workspace
819 994 int* ancestor = wrk_int;
... ... @@ -836,7 +1011,9 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent,
836 1011 int* ATp = AT.m_colstart;
837 1012 int* ATi = AT.m_rowinx;
838 1013  
839   - // Initialise
  1014 + // Initialise. Note that init_ata wires the working array into the *head
  1015 + // and *next pointers, so both arrays are just aliases of the wrk_int
  1016 + // array.
840 1017 int* head = NULL;
841 1018 int* next = NULL;
842 1019 if (ata) init_ata(&AT, post, wrk_int, &head, &next);
... ... @@ -845,8 +1022,6 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent,
845 1022 for (i = 0; i < n; i++) ancestor[i] = i;
846 1023  
847 1024 // Loop over all columns of matrix A (or rows of matrix AT)
848   - #define HEAD(k,j) (ata ? head [k] : j)
849   - #define NEXT(J) (ata ? next [J] : -1)
850 1025 for (k = 0; k < n; k++) {
851 1026 j = post[k]; // j is the kth node in postordered etree
852 1027 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,
872 1047 delete [] wrk_int;
873 1048  
874 1049 // Return result
875   - return (colcount);
  1050 + return colcount;
876 1051 }
877 1052  
878 1053  
... ... @@ -885,7 +1060,7 @@ int* GSparseSymbolic::cs_counts(GSparseMatrix* A, const int* parent,
885 1060 * ata flag *
886 1061 * Output: etree(A) Evaluation tree of matrix A *
887 1062 ***************************************************************************/
888   -int* GSparseSymbolic::cs_etree(GSparseMatrix* A, int ata)
  1063 +int* GSparseSymbolic::cs_etree(const GSparseMatrix* A, int ata)
889 1064 {
890 1065 // Declare loop variables
891 1066 int i, k, p;
... ... @@ -910,7 +1085,10 @@ int* GSparseSymbolic::cs_etree(GSparseMatrix* A, int ata)
910 1085 int* prev = wrk_int + n;
911 1086  
912 1087 // If 'prev' is requested, initialise array with -1
913   - if (ata) for (i = 0; i < m ; ++i) prev[i] = -1;
  1088 + if (ata) {
  1089 + for (i = 0; i < m ; ++i)
  1090 + prev[i] = -1;
  1091 + }
914 1092  
915 1093 // Loop over all nodes (columns)
916 1094 for (k = 0; k < n; ++k) {
... ... @@ -940,7 +1118,7 @@ int* GSparseSymbolic::cs_etree(GSparseMatrix* A, int ata)
940 1118 delete [] wrk_int;
941 1119  
942 1120 // Return result
943   - return (parent);
  1121 + return parent;
944 1122  
945 1123 }
946 1124  
... ... @@ -956,7 +1134,8 @@ int* GSparseSymbolic::cs_etree(GSparseMatrix* A, int ata)
956 1134 * other Other function *
957 1135 * Output: result 0..n-1 if ok, -1 in case of error *
958 1136 ***************************************************************************/
959   -int GSparseSymbolic::cs_fkeep(GSparseMatrix* A, int(*fkeep)(int, int, double, void*),
  1137 +int GSparseSymbolic::cs_fkeep(GSparseMatrix* A,
  1138 + int(*fkeep)(int, int, double, void*),
960 1139 void* other)
961 1140 {
962 1141 // 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
984 1163 }
985 1164  
986 1165 // Finalise A
987   - Ap[n] = nz;
  1166 + Ap[n] = nz;
  1167 + A->m_elements = nz;
988 1168  
989 1169 // Remove extra space from A
990   -// cs_sprealloc (A, 0) ; /* remove extra space from A */
  1170 +// cs_sprealloc(A, 0);
  1171 + A->free_elements(nz, (A->m_elements-nz));
991 1172  
992 1173 // Return number of non-zero elements
993   - return (nz) ;
  1174 + return nz;
994 1175 }
995 1176  
996 1177  
... ... @@ -1061,7 +1242,7 @@ int GSparseSymbolic::cs_leaf(int i, int j, const int *first, int *maxfirst,
1061 1242 int* GSparseSymbolic::cs_pinv(int const *p, int n)
1062 1243 {
1063 1244 // Return NULL pointer if input pointer is NULL. This denotes identity
1064   - if (!p || n < 1) return (NULL);
  1245 + if (!p || n < 1) return NULL;
1065 1246  
1066 1247 // Allocate result array
1067 1248 int* pinv = new int[n];
... ... @@ -1073,7 +1254,7 @@ int* GSparseSymbolic::cs_pinv(int const *p, int n)
1073 1254 pinv[p[k]] = k;
1074 1255  
1075 1256 // Return result
1076   - return (pinv);
  1257 + return pinv;
1077 1258 }
1078 1259  
1079 1260  
... ... @@ -1137,90 +1318,7 @@ int* GSparseSymbolic::cs_post(const int* parent, int n)
1137 1318 delete [] wrk_int;
1138 1319  
1139 1320 // Return result
1140   - return (post);
1141   -
1142   -}
1143   -
1144   -
1145   -/***************************************************************************
1146   - * cs_symperm *
1147   - * ----------------------------------------------------------------------- *
1148   - * C = A(p,p) where A and C are symmetric the upper part stored. *
1149   - * pinv not p *
1150   - * ----------------------------------------------------------------------- *
1151   - * Input: A Sparse matrix *
1152   - * pinv pinv[0..n-1] *
1153   - * values =1 allocate memory for values (NOT USED) *
1154   - * Output: C Sparse matrix *
1155   - ***************************************************************************/
1156   -GSparseMatrix GSparseSymbolic::cs_symperm(GSparseMatrix* A, const int* pinv,
1157   - int values) const
1158   -{
1159   - // Declare loop variables
1160   - int i, j, p, q, i2, j2;
1161   -
1162   - // Assign matrix attributes
1163   - int n = A->m_cols;
1164   - int* Ap = A->m_colstart;
1165   - int* Ai = A->m_rowinx;
1166   - double* Ax = A->m_data;
1167   -
1168   - // Allocate result matrix
1169   - GSparseMatrix C(n, n, Ap[n]);
1170   -
1171   - // Allocate and initialise workspace
1172   - int wrk_size = n;
1173   - int* wrk_int = new int[wrk_size];
1174   - if (wrk_int == NULL)
1175   - throw mem_alloc("GSparseSymbolic::cs_symperm(GSparseMatrix*, const int*, int",
1176   - wrk_size);
1177   - for (i = 0; i < wrk_size; ++i)
1178   - wrk_int[i] = 0;
1179   -
1180   - // Assign result matriix attributes
1181   - int* Cp = C.m_colstart;
1182   - int* Ci = C.m_rowinx;
1183   - double* Cx = C.m_data;
1184   -
1185   - // Count entries in each column of C
1186   - for (j = 0; j < n; j++) {
1187   -
1188   - // Column j of A is column j2 of C
1189   - j2 = pinv ? pinv[j] : j;
1190   -
1191   - // Loop over entries in column j
1192   - for (p = Ap[j]; p < Ap[j+1]; p++) {
1193   - i = Ai [p];
1194   - if (i > j) continue; // skip lower triangular part of A
1195   - i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C
1196   - wrk_int[CS_MAX(i2, j2)]++; // column count of C
1197   - }
1198   - }
1199   -
1200   - // Compute column pointers of C
1201   - cs_cumsum(Cp, wrk_int, n);
1202   -
1203   - // Loop over all columns of A
1204   - for (j = 0 ; j < n ; j++) {
1205   -
1206   - // Column j of A is column j2 of C
1207   - j2 = pinv ? pinv[j] : j;
1208   -
1209   - // Loop over entries in column j
1210   - for (p = Ap[j]; p < Ap[j+1]; p++) {
1211   - i = Ai [p] ;
1212   - if (i > j) continue; // skip lower triangular part of A
1213   - i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C
1214   - Ci[q = wrk_int[CS_MAX(i2,j2)]++] = CS_MIN(i2,j2);
1215   - Cx[q] = Ax[p];
1216   - }
1217   - }
1218   -
1219   - // Free workspace
1220   - delete [] wrk_int;
1221   -
1222   - // Return result
1223   - return (C);
  1321 + return post;
1224 1322  
1225 1323 }
1226 1324  
... ... @@ -1266,7 +1364,7 @@ int GSparseSymbolic::cs_tdfs(int j, int k, int* head, const int* next,
1266 1364 }
1267 1365  
1268 1366 // Return result
1269   - return (k) ;
  1367 + return k;
1270 1368 }
1271 1369  
1272 1370  
... ... @@ -1290,7 +1388,7 @@ int GSparseSymbolic::cs_tdfs(int j, int k, int* head, const int* next,
1290 1388 * next pointer to next[0..n-1] linked list info *
1291 1389 * Output: - - *
1292 1390 ***************************************************************************/
1293   -void GSparseSymbolic::init_ata(GSparseMatrix* AT, const int* post,
  1391 +void GSparseSymbolic::init_ata(const GSparseMatrix* AT, const int* post,
1294 1392 int* wrk_int, int** head, int** next)
1295 1393 {
1296 1394 // Declare loop variables
... ... @@ -1321,6 +1419,9 @@ void GSparseSymbolic::init_ata(GSparseMatrix* AT, const int* post,
1321 1419 (*next)[i] = (*head)[k];
1322 1420 (*head)[k] = i;
1323 1421 }
  1422 +
  1423 + // Return
  1424 + return;
1324 1425 }
1325 1426  
1326 1427  
... ...
src/GSparseSymbolic.hpp
  View file @cb738b2
... ... @@ -61,7 +61,7 @@ public:
61 61 // Unary operators
62 62  
63 63 // Functions
64   - void cholesky_symbolic_analysis(int order, GSparseMatrix& m);
  64 + void cholesky_symbolic_analysis(int order, const GSparseMatrix& m);
65 65  
66 66 // Exception: Invalid ordering scheme
67 67 class invalid_order : public GException {
... ... @@ -71,18 +71,17 @@ public:
71 71  
72 72 private:
73 73 // Functions
74   - int* cs_amd (int order, GSparseMatrix* A);
75   - int* cs_counts(GSparseMatrix* A, const int* parent, const int* post, int ata);
76   - int* cs_etree(GSparseMatrix* A, int ata);
77   - int cs_fkeep(GSparseMatrix* A, int(*fkeep)(int, int, double, void*), void* other);
78   - int cs_leaf(int i, int j, const int* first, int* maxfirst, int* prevleaf, int* ancestor, int* jleaf);
79   - int* cs_pinv(int const* p, int n);
80   - int* cs_post(const int* parent, int n);
81   - GSparseMatrix cs_symperm(GSparseMatrix* A, const int* pinv, int values) const;
82   - int cs_tdfs(int j, int k, int* head, const int* next, int* post, int* stack);
83   - static void init_ata(GSparseMatrix* AT, const int* post, int* wrk_int, int** head, int** next);
84   - static int cs_diag(int i, int j, double aij, void* other);
85   - static int cs_wclear(int mark, int lemax, int* w, int n);
  74 + int* cs_amd (int order, const GSparseMatrix* A);
  75 + int* cs_counts(const GSparseMatrix* A, const int* parent, const int* post, int ata);
  76 + int* cs_etree(const GSparseMatrix* A, int ata);
  77 + int cs_fkeep(GSparseMatrix* A, int(*fkeep)(int, int, double, void*), void* other);
  78 + int cs_leaf(int i, int j, const int* first, int* maxfirst, int* prevleaf, int* ancestor, int* jleaf);
  79 + int* cs_pinv(int const* p, int n);
  80 + int* cs_post(const int* parent, int n);
  81 + int cs_tdfs(int j, int k, int* head, const int* next, int* post, int* stack);
  82 + static void init_ata(const GSparseMatrix* AT, const int* post, int* wrk_int, int** head, int** next);
  83 + static int cs_diag(int i, int j, double aij, void* other);
  84 + static int cs_wclear(int mark, int lemax, int* w, int n);
86 85  
87 86 // Data
88 87 int* m_pinv; // Inverse row permutation for QR, fill reduce permutation for Cholesky
... ...
src/GVector.cpp
  View file @cb738b2
... ... @@ -46,15 +46,32 @@ GVector::GVector(int num)
46 46  
47 47  
48 48 /***************************************************************************
  49 + * GVector copy constructor *
  50 + ***************************************************************************/
  51 +GVector::GVector(const GVector& v)
  52 +{
  53 + // Allocate vector array. Throw exception if allocation failed
  54 + m_data = new double[v.m_num];
  55 + if (m_data == NULL)
  56 + throw mem_alloc("GVector copy constructor", v.m_num);
  57 +
  58 + // Store vector size and copy elements
  59 + m_num = v.m_num;
  60 + for (int i = 0; i < m_num; ++i)
  61 + m_data[i] = v.m_data[i];
  62 +
  63 + // Return
  64 + return;
  65 +}
  66 +
  67 +
  68 +/***************************************************************************
49 69 * GVector destructor *
50 70 ***************************************************************************/
51 71 GVector::~GVector()
52 72 {
53 73 // Deallocate only if memory has indeed been allocated
54   - if (m_num > 0) {
55   - m_num = 0;
56   - delete[] m_data;
57   - }
  74 + if (m_data != NULL) delete[] m_data;
58 75  
59 76 // Return
60 77 return;
... ... @@ -67,6 +84,33 @@ GVector::~GVector()
67 84 = =
68 85 ==========================================================================*/
69 86  
  87 +/***************************************************************************
  88 + * GVector assignment operator *
  89 + ***************************************************************************/
  90 +GVector& GVector::operator= (const GVector& v)
  91 +{
  92 + // Assign only if not identical
  93 + if (this != &v) {
  94 +
  95 + // First delete old data if it exists
  96 + if (m_data != NULL) delete [] m_data;
  97 +
  98 + // Allocate vector array. Throw exception if allocation failed
  99 + m_data = new double[v.m_num];
  100 + if (m_data == NULL)
  101 + throw mem_alloc("GVector assignment operator", v.m_num);
  102 +
  103 + // Store vector size and copy elements
  104 + m_num = v.m_num;
  105 + for (int i = 0; i < m_num; ++i)
  106 + m_data[i] = v.m_data[i];
  107 + }
  108 +
  109 + // Return vector
  110 + return *this;
  111 +}
  112 +
  113 +
70 114 /*==========================================================================
71 115 = =
72 116 = GVector member functions =
... ...
src/GVector.hpp
  View file @cb738b2
... ... @@ -134,18 +134,6 @@ private:
134 134 * Inline members *
135 135 ***************************************************************************/
136 136  
137   -// Copy constructor
138   -inline
139   -GVector::GVector(const GVector& v)
140   -{
141   - m_data = new double[v.m_num];
142   - if (m_data == NULL)
143   - throw mem_alloc("GVector copy constructor", v.m_num);
144   - m_num = v.m_num;
145   - for (int i = 0; i < m_num; ++i)
146   - m_data[i] = v.m_data[i];
147   -}
148   -
149 137 // Vector element access operator
150 138 inline
151 139 double& GVector::operator() (int inx)
... ... @@ -168,22 +156,6 @@ const double&amp; GVector::operator() (int inx) const
168 156 return m_data[inx];
169 157 }
170 158  
171   -// Vector assignment operator
172   -inline
173   -GVector& GVector::operator= (const GVector& v)
174   -{
175   - if (this != &v) {
176   - if (m_data != NULL) delete [] m_data;
177   - m_data = new double[v.m_num];
178   - if (m_data == NULL)
179   - throw mem_alloc("GVector assignment operator", v.m_num);
180   - m_num = v.m_num;
181   - for (int i = 0; i < m_num; ++i)
182   - m_data[i] = v.m_data[i];
183   - }
184   - return *this;
185   -}
186   -
187 159 // Vector unary addition operator
188 160 inline
189 161 GVector& GVector::operator+= (const GVector& v)
... ...