From d58bb54e7fe600232a4ecd882445b33446172b86 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Sun, 18 Jan 2009 18:00:19 +0000
Subject: [PATCH] add a smart realloc algorithm when filling a sparse matrix

---
 Eigen/src/Sparse/CompressedStorage.h |  3 ++-
 Eigen/src/Sparse/SparseMatrix.h      | 21 ++++++++++++++++-----
 Eigen/src/Sparse/SparseVector.h      |  4 ++--
 3 files changed, 20 insertions(+), 8 deletions(-)

diff --git a/Eigen/src/Sparse/CompressedStorage.h b/Eigen/src/Sparse/CompressedStorage.h
index f6f1ae687..da85aad0b 100644
--- a/Eigen/src/Sparse/CompressedStorage.h
+++ b/Eigen/src/Sparse/CompressedStorage.h
@@ -82,7 +82,7 @@ class CompressedStorage
         reallocate(m_size);
     }
 
-    void resize(int size, int reserveSizeFactor = 0)
+    void resize(int size, float reserveSizeFactor = 0)
     {
       if (m_allocatedSize<size)
         reallocate(size + reserveSizeFactor*size);
@@ -98,6 +98,7 @@ class CompressedStorage
     }
 
     int size() const { return m_size; }
+    int allocatedSize() const { return m_allocatedSize; }
     void clear() { m_size = 0; }
 
     Scalar& value(int i) { return m_values[i]; }
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index d3b504839..d90a7c69f 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -200,7 +200,21 @@ class SparseMatrix
       int startId = m_outerIndex[outer];
       int id = m_outerIndex[outer+1]-1;
       ++m_outerIndex[outer+1];
-      m_data.resize(id+2);
+      
+      float reallocRatio = 1;
+      if (m_data.allocatedSize()<id+2)
+      {
+        // we need to reallocate the data, to reduce multiple reallocations
+        // we use a smart resize algorithm based on the current filling ratio
+        // we use float to avoid overflows
+        float nnzEstimate = float(m_outerIndex[outer])*float(m_outerSize)/float(outer);
+        reallocRatio = (nnzEstimate-float(m_data.size()))/float(m_data.size());
+        // let's bounds the realloc ratio to
+        //   1) reduce multiple minor realloc when the matrix is almost filled
+        //   2) avoid to allocate too much memory when the matrix is almost empty
+        reallocRatio = std::min(std::max(reallocRatio,1.5f),8.f);
+      }
+      m_data.resize(id+2,reallocRatio);
 
       while ( (id >= startId) && (m_data.index(id) > inner) )
       {
@@ -209,10 +223,7 @@ class SparseMatrix
         --id;
       }
       m_data.index(id+1) = inner;
-      //return (m_data.value(id+1) = 0);
-      m_data.value(id+1) = 0;
-//       std::cerr << m_outerIndex[outer] << " " << m_outerIndex[outer+1] << "\n";
-      return m_data.value(id+1);
+      return (m_data.value(id+1) = 0);
     }
 
 //     inline void
diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h
index 901332e08..479b678e1 100644
--- a/Eigen/src/Sparse/SparseVector.h
+++ b/Eigen/src/Sparse/SparseVector.h
@@ -151,7 +151,7 @@ class SparseVector
     {
       int startId = 0;
       int id = m_data.size() - 1;
-      m_data.resize(id+2);
+      m_data.resize(id+2,1);
 
       while ( (id >= startId) && (m_data.index(id) > i) )
       {
@@ -172,7 +172,7 @@ class SparseVector
 
     void resizeNonZeros(int size) { m_data.resize(size); }
 
-    inline SparseVector() : m_size(0) { resize(0, 0); }
+    inline SparseVector() : m_size(0) { resize(0); }
 
     inline SparseVector(int size) : m_size(0) { resize(size); }