Memory leaks and initialization (#73)

* initial work on sparse matrix representation

* store/restore functionality

* addLastRow functionality

* getRow and getColumn

* column-merging functionality

* added an interface class

* introducing also sparse vectors

* added addLastColumn functionality

* another unittest

* get sparse columns/matrices in dense form

* WIP on storing the constraint matrix inside the tableau in sparse form

* more WIP, fixed a few bugs, still have a couple of failing tests

* fixed some test issues

* initialization

* resize _a along with the rest

* sparse lu factors

* store also the transposed versions of F and V

* starting work on the sparse GE

* some work on changing the values within an existing sparse
representation, needed for sparse factorization.
still WIP

* refactoring and new functionality for CSRMatrix: any kind of
insertions and deletions

* support for empty initialization and counting elements

* sparse GE is now working. minor bug fixes elsewhere

* compute Ft and Vt as part of the G-elimination process

* tests

* basis oracles can return also sparse columns, not just dense

* sparse LU factorization class

* switch to using the sparse factorization in the engine/tableau

* bug fix in mergeColumns, and nicer printing

* bug fix

* bug fix

* bug fix: merging columns does not delete the actual column, just
leaves it empty

* configuration changes

* optimization: since the sparse columns of A are needed all the time,
just compute them once-and-for-all

* a more efficient implementation of sparse vectors

* comments and unit tests

* cleanup

* keep _A in dense column-major format, too, instead of repeatedly
invoking toDense() to gets its columns

* bad deletes

* bug fix in test

* bug fixes: relu constraint propagation, and the handling of merged
variables in the preprocessor

* new test

* compute Ft incrementally, use it whenever sparse columns are required

* did the todo

* valgrind fixes

src/basis_factorization/IBasisFactorization.h

@@ -108,6 +108,11 @@ class IBasisFactorization
      */
     virtual void invertBasis( double *result ) = 0;
 
+    /*
+      For debugging
+    */
+    virtual void dump() const {};
+
 private:
     /*
       A flag that controls whether factorization is enabled or

src/basis_factorization/SparseGaussianEliminator.cpp

@@ -82,6 +82,7 @@ void SparseGaussianEliminator::initializeFactorization( const SparseMatrix *A, S
         so we just leave it empty for now.
     */
     A->storeIntoOther( _sparseLUFactors->_V );
+    _sparseLUFactors->_V->transposeIntoOther( _sparseLUFactors->_Vt );
     _sparseLUFactors->_F->initializeToEmpty( _m, _m );
     _sparseLUFactors->_P.resetToIdentity();
     _sparseLUFactors->_Q.resetToIdentity();
@@ -148,12 +149,9 @@ void SparseGaussianEliminator::run( const SparseMatrix *A, SparseLUFactors *spar
         eliminate();
     }
 
-    // Execute the changes in F
+    // Execute the changes in F, compute its transpose
     _sparseLUFactors->_F->executeChanges();
-
-    // Compute the transposed F, V
     _sparseLUFactors->_F->transposeIntoOther( _sparseLUFactors->_Ft );
-    _sparseLUFactors->_V->transposeIntoOther( _sparseLUFactors->_Vt );
 
     // DEBUG({
     //         // Check that the factorization is correct
@@ -229,7 +227,7 @@ void SparseGaussianEliminator::choosePivot()
             _vPivotColumn = _sparseLUFactors->_Q._rowOrdering[i];
 
             // Get the singleton element
-            _sparseLUFactors->_V->getColumn( _vPivotColumn, &sparseColumn );
+            _sparseLUFactors->_Vt->getRow( _vPivotColumn, &sparseColumn );
 
             // There may be some elements in higher rows - we need just the one
             // in the active submatrix.
@@ -276,7 +274,7 @@ void SparseGaussianEliminator::choosePivot()
     for ( unsigned uColumn = _eliminationStep; uColumn < _m; ++uColumn )
     {
         unsigned vColumn = _sparseLUFactors->_Q._rowOrdering[uColumn];
-        _sparseLUFactors->_V->getColumn( vColumn, &sparseColumn );
+        _sparseLUFactors->_Vt->getRow( vColumn, &sparseColumn );
 
         double maxInColumn = 0;
         for ( unsigned entry = 0; entry < sparseColumn.getNnz(); ++entry )
@@ -352,7 +350,7 @@ void SparseGaussianEliminator::eliminate()
       Eliminate all entries below the pivot element U[k,k]
       We know that V[_vPivotRow, _vPivotColumn] = U[k,k].
     */
-    _sparseLUFactors->_V->getColumn( _vPivotColumn, &sparseColumn );
+    _sparseLUFactors->_Vt->getRow( _vPivotColumn, &sparseColumn );
 
     // Get the pivot row in dense format, due to repeated access
     _sparseLUFactors->_V->getRowDense( _vPivotRow, _work );
@@ -389,6 +387,7 @@ void SparseGaussianEliminator::eliminate()
         --_numUColumnElements[_eliminationStep];
         --_numURowElements[uRow];
         _sparseLUFactors->_V->commitChange( vRow, _vPivotColumn, 0.0 );
+        _sparseLUFactors->_Vt->commitChange( _vPivotColumn, vRow, 0.0 );
         _sparseLUFactors->_V->getRow( vRow, &sparseRow );
 
         Set<unsigned> columnsAlreadyHandled;
@@ -418,6 +417,7 @@ void SparseGaussianEliminator::eliminate()
             }
 
             _sparseLUFactors->_V->commitChange( vRow, vColumnIndex, newValue );
+            _sparseLUFactors->_Vt->commitChange( vColumnIndex, vRow, newValue );
         }
 
         // Next, handle entries that were zero in the eliminated row
@@ -433,7 +433,9 @@ void SparseGaussianEliminator::eliminate()
 
             ++_numUColumnElements[uColumnIndex];
             ++_numURowElements[uRow];
-            _sparseLUFactors->_V->commitChange( vRow, vColumnIndex, rowMultiplier * _work[vColumnIndex] );
+            double newVal = rowMultiplier * _work[vColumnIndex];
+            _sparseLUFactors->_V->commitChange( vRow, vColumnIndex, newVal );
+            _sparseLUFactors->_Vt->commitChange( vColumnIndex, vRow, newVal );
         }
 
         /*
@@ -445,6 +447,7 @@ void SparseGaussianEliminator::eliminate()
 
     // Execute the changes in V
     _sparseLUFactors->_V->executeChanges();
+    _sparseLUFactors->_Vt->executeChanges();
 }
 
 void SparseGaussianEliminator::log( const String &message )

src/basis_factorization/SparseLUFactorization.cpp

@@ -52,6 +52,12 @@ void SparseLUFactorization::freeIfNeeded()
 		_B = NULL;
 	}
 
+    if ( _z )
+    {
+        delete[] _z;
+        _z = NULL;
+    }
+
     List<EtaMatrix *>::iterator it;
     for ( it = _etas.begin(); it != _etas.end(); ++it )
         delete *it;

src/basis_factorization/tests/Test_SparseGaussianEliminator.h

@@ -301,11 +301,6 @@ class SparseGaussianEliminatorTestSuite : public CxxTest::TestSuite
             TS_ASSERT_THROWS_NOTHING( delete ge );
         }
     }
-
-    void test_todo()
-    {
-        TS_TRACE( "Currently we don't update Vt, Ft during the factorization, and sometimes we call getColumn() which is inefficient. Consider changing." );
-    }
 };
 
 //

src/engine/Tableau.cpp

@@ -923,6 +923,13 @@ void Tableau::pickLeavingVariable( double *changeColumn )
                 {
                     _changeRatio = ratio;
                     _leavingVariable = i;
+
+                    if ( FloatUtils::isZero( _changeRatio ) )
+                    {
+                        // All ratios are non-positive, so if we already hit 0 we are done.
+                        _changeRatio = 0.0;
+                        break;
+                    }
                 }
             }
         }
@@ -950,6 +957,13 @@ void Tableau::pickLeavingVariable( double *changeColumn )
                 {
                     _changeRatio = ratio;
                     _leavingVariable = i;
+
+                    if ( FloatUtils::isZero( _changeRatio ) )
+                    {
+                        // All ratios are non-negative, so if we already hit 0 we are done.
+                        _changeRatio = 0.0;
+                        break;
+                    }
                 }
             }
         }
@@ -1483,7 +1497,11 @@ void Tableau::addRow()
         throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "Tableau::newDenseA" );
 
     for ( unsigned column = 0; column < _n; ++column )
+    {
         memcpy( newDenseA + ( column * newM ), _denseA + ( column * _m ), sizeof(double) * _m );
+        newDenseA[column*newM + newM - 1] = 0.0;
+    }
+
     delete[] _denseA;
     _denseA = newDenseA;