diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/gnp/triple_range_distance_sq.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/gnp/triple_range_distance_sq.h
index c43c6e73a6..4885ed4dbb 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/gnp/triple_range_distance_sq.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/gnp/triple_range_distance_sq.h
@@ -41,16 +41,16 @@ class TripleRangeDistanceSq {
             num_tuples_[1] =
               num_tuples_[2] =
                 core::math::BinomialCoefficient<double>(
-                  table_in.get_node_count(nodes_[0]) - 1, 2);
+                  nodes_[0]->count() - 1, 2);
         }
 
         // node_0 = node_1, node_1 \not = node_2
         else {
           num_tuples_[0] = num_tuples_[1] =
-                             (table_in.get_node_count(nodes_[0]) - 1) *
-                             table_in.get_node_count(nodes_[2]);
+                             (nodes_[0]->count() - 1) *
+                             nodes_[2]->count();
           num_tuples_[2] = core::math::BinomialCoefficient<double>(
-                             table_in.get_node_count(nodes_[0]), 2);
+                             nodes_[0]->count(), 2);
         }
       }
       else {
@@ -58,20 +58,17 @@ class TripleRangeDistanceSq {
         // node_0 \not = node_1, node_1 = node_2
         if(nodes_[1] == nodes_[2]) {
           num_tuples_[1] = num_tuples_[2] =
-                             (table_in.get_node_count(nodes_[1]) - 1) *
-                             table_in.get_node_count(nodes_[0]);
+                             (nodes_[1]->count() - 1) *
+                             nodes_[0]->count();
           num_tuples_[0] = core::math::BinomialCoefficient<double>(
-                             table_in.get_node_count(nodes_[1]), 2);
+                             nodes_[1]->count(), 2);
         }
 
         // node_0 \not = node_1, node_1 \not = node_2
         else {
-          num_tuples_[0] = table_in.get_node_count(nodes_[1]) *
-                           table_in.get_node_count(nodes_[2]);
-          num_tuples_[1] = table_in.get_node_count(nodes_[0]) *
-                           table_in.get_node_count(nodes_[2]);
-          num_tuples_[2] = table_in.get_node_count(nodes_[0]) *
-                           table_in.get_node_count(nodes_[1]);
+          num_tuples_[0] = nodes_[1]->count() * nodes_[2]->count();
+          num_tuples_[1] = nodes_[0]->count() * nodes_[2]->count();
+          num_tuples_[2] = nodes_[0]->count() * nodes_[1]->count();
         }
       }
     }
@@ -138,8 +135,7 @@ class TripleRangeDistanceSq {
       int node_index_in) {
 
       nodes_[node_index_in] = new_node_in;
-      const typename TreeType::BoundType &new_node_bound =
-        table_in.get_node_bound(new_node_in);
+      const typename TreeType::BoundType &new_node_bound = new_node_in->bound();
 
       for(int existing_node_index = node_index_in + 1;
           existing_node_index < 3; existing_node_index++) {
@@ -147,8 +143,7 @@ class TripleRangeDistanceSq {
         // Change for the first existing node.
         core::math::Range existing_range_distance_sq =
           new_node_bound.RangeDistanceSq(
-            metric_in,
-            table_in.get_node_bound(nodes_[existing_node_index]));
+            metric_in, nodes_[existing_node_index]->bound());
         set_range_distance_sq(
           node_index_in, existing_node_index, existing_range_distance_sq);
       }
@@ -164,8 +159,7 @@ class TripleRangeDistanceSq {
       int node_index_in) {
 
       nodes_[node_index_in] = new_node_in;
-      const typename TreeType::BoundType &new_node_bound =
-        table_in.get_node_bound(new_node_in);
+      const typename TreeType::BoundType &new_node_bound = new_node_in->bound();
 
       for(int existing_node_index = 0; existing_node_index < node_index_in;
           existing_node_index++) {
@@ -173,8 +167,7 @@ class TripleRangeDistanceSq {
         // Change for the first existing node.
         core::math::Range existing_range_distance_sq =
           new_node_bound.RangeDistanceSq(
-            metric_in,
-            table_in.get_node_bound(nodes_[existing_node_index]));
+            metric_in, nodes_[existing_node_index]->bound());
         set_range_distance_sq(
           node_index_in, existing_node_index, existing_range_distance_sq);
       }
@@ -191,11 +184,9 @@ class TripleRangeDistanceSq {
         nodes_[j] = nodes_in[j];
       }
       for(unsigned int j = 0; j < nodes_.size(); j++) {
-        const typename TreeType::BoundType &outer_bound =
-          table.get_node_bound(nodes_[j]);
+        const typename TreeType::BoundType &outer_bound = nodes_[j]->bound();
         for(unsigned int i = j + 1; i < nodes_.size(); i++) {
-          const typename TreeType::BoundType &inner_bound =
-            table.get_node_bound(nodes_[i]);
+          const typename TreeType::BoundType &inner_bound = nodes_[i]->bound();
           core::math::Range range_distance_sq =
             outer_bound.RangeDistanceSq(metric_in, inner_bound);
           set_range_distance_sq(i, j, range_distance_sq);
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/math/linear_algebra.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/math/linear_algebra.h
index 3d084752c7..78ede86137 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/math/linear_algebra.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/math/linear_algebra.h
@@ -135,9 +135,9 @@ template<typename VectorType>
 static void AddTo(
   const VectorType &vec_in, VectorType *vec_out) {
 
-  for(unsigned int i = 0; i < vec_in.n_elem; i++) {
-    (*vec_out)[i] += vec_in[i];
-  }
+  arma::vec vec_in_alias(vec_in.ptr(), vec_in.length());
+  arma::vec vec_out_alias(vec_out->ptr(), vec_in.length(), false);
+  vec_out_alias = vec_out_alias + vec_in_alias;
 }
 
 /** @brief Computes $c = c + \alpha * a b^T$.
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h
index e3be7ff57c..fdc243d882 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h
@@ -74,8 +74,7 @@ void TrustRegion::ComputeSteihaugDirection_(
     double alpha = arma::dot(r, r) / quadratic_form.at(0, 0);
 
     // z_{j + 1} = z_j + \alpha_j d_j
-    core::math::ScaleOverwrite(alpha, d, &z_next);
-    core::math::AddTo(z, &z_next);
+    z_next = z + alpha * d;
 
     // If the z_{j + 1} violates the trust region bound,
     if(arma::norm(z_next, 2) >= radius) {
@@ -104,8 +103,7 @@ void TrustRegion::ComputeSteihaugDirection_(
     double beta_next = arma::dot(r_next, r_next) / arma::dot(r, r);
 
     // d_{j + 1} = -r_{j + 1} + \beta_{j + 1} d_j
-    core::math::ScaleOverwrite(beta_next, d, &d_next);
-    core::math::SubFrom(r_next, &d_next);
+    d_next = beta_next * d - r_next;
 
     // Set the variables for the next iteration.
     core::math::CopyValues(r_next, &r);
@@ -141,7 +139,7 @@ void TrustRegion::ComputeCauchyPoint_(
             1.0, core::math::Pow<3, 1>(gradient_norm) / (
               radius * quadratic_form.at(0, 0)));
   }
-  core::math::ScaleInit(- tau * radius / gradient_norm , gradient, p);
+  (*p) = (- tau * radius / gradient_norm)  *  gradient;
 }
 
 void TrustRegion::ComputeDoglegDirection_(
@@ -200,7 +198,7 @@ void TrustRegion::ComputeDoglegDirection_(
       // If the norm of p_u is beyond the trust radius, then the
       // solution lies on the boundary along p_u.
       if(p_u_norm >= radius) {
-        core::math::ScaleInit(radius / p_u_norm, p_u, p);
+        (*p) = (radius / p_u_norm) * p_u;
       }
 
       // Otherwise the quadratic equation composed of the gradient and
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/table/dense_point.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/table/dense_point.h
index b4a6661bc7..59dee7fe50 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/table/dense_point.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/table/dense_point.h
@@ -167,46 +167,6 @@ class DensePoint {
       is_alias_ = true;
     }
 
-    void Add(
-      double scale_factor, const core::table::DensePoint &point_in) {
-      for(int i = 0; i < point_in.length(); i++) {
-        ptr_.get()[i] += scale_factor * point_in[i];
-      }
-    }
-
-    void SubOverwrite(
-      const core::table::DensePoint &subtracted,
-      const core::table::DensePoint &subtract_from) {
-      for(int i = 0; i < n_rows_; i++) {
-        ptr_.get()[i] = subtract_from[i] - subtracted[i];
-      }
-    }
-
-    void ScaleOverwrite(
-      double scale_in, const core::table::DensePoint &point_in) {
-      for(int i = 0; i < n_rows_; i++) {
-        ptr_.get()[i] = scale_in * point_in[i];
-      }
-    }
-
-    void operator+=(const core::table::DensePoint &point_in) {
-      for(int i = 0; i < point_in.length(); i++) {
-        ptr_.get()[i] += point_in[i];
-      }
-    }
-
-    void operator/=(double scale_factor) {
-      for(int i = 0; i < n_rows_; i++) {
-        ptr_.get()[i] /= scale_factor;
-      }
-    }
-
-    void operator*=(double scale_factor) {
-      for(int i = 0; i < n_rows_; i++) {
-        ptr_.get()[i] *= scale_factor;
-      }
-    }
-
     void Print() const {
       printf("Vector of length: %d\n", n_rows_);
       for(int i = 0; i < n_rows_; i++) {
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/table/distributed_table.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/table/distributed_table.h
index 777f76d641..41d57b044f 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/table/distributed_table.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/table/distributed_table.h
@@ -163,10 +163,12 @@ class DistributedTable: public boost::noncopyable {
         tmp_point.SetZero();
         for(int i = 0; i < num_samples; i++) {
           tmp_point.SetZero();
-          tmp_point += top_leaf_nodes[
-                         core::math::RandInt(top_leaf_nodes.size())]->bound().center();
+          core::math::AddTo(
+            top_leaf_nodes[
+              core::math::RandInt(top_leaf_nodes.size())]->bound().center(),
+            &tmp_point);
         }
-        tmp_point /= static_cast<double>(num_samples);
+        core::math::Scale(1.0 / static_cast<double>(num_samples), &tmp_point);
         top_leaf_nodes.push_back(new TreeType());
         top_leaf_nodes[ top_leaf_nodes.size() - 1 ]->bound().center().Copy(
           tmp_point);
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/distributed_local_kmeans.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/distributed_local_kmeans.h
index 3c8326dcd2..699bf0de92 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/distributed_local_kmeans.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/distributed_local_kmeans.h
@@ -8,6 +8,7 @@
 
 #include <boost/mpi.hpp>
 #include <boost/serialization/string.hpp>
+#include "core/math/linear_algebra.h"
 #include "core/table/dense_point.h"
 
 namespace core {
@@ -52,16 +53,17 @@ class DistributedLocalKMeans {
           double factor =
             static_cast<double>(num_points_) /
             static_cast<double>(num_points_ + centroid_in.num_points());
-          centroid_ *= factor;
-          centroid_.Add(1.0 - factor, centroid_in.centroid());
+          core::math::Scale(factor, &centroid_);
+          core::math::AddExpert(
+            1.0 - factor, centroid_in.centroid(), &centroid_);
           num_points_ = num_points_ + centroid_in.num_points();
         }
 
         void Add(const core::table::DensePoint &point_in) {
           double factor = static_cast<double>(num_points_) /
                           static_cast<double>(num_points_ + 1);
-          centroid_ *= factor;
-          centroid_.Add(1.0 - factor, point_in);
+          core::math::Scale(factor, &centroid_);
+          core::math::AddExpert(1.0 - factor, point_in, &centroid_);
           num_points_++;
         }
 
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/gen_metric_tree.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/gen_metric_tree.h
index e5e41817e4..d1fb53c219 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/gen_metric_tree.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/gen_metric_tree.h
@@ -9,6 +9,7 @@
 #include <vector>
 #include "ball_bound.h"
 #include "general_spacetree.h"
+#include "core/math/linear_algebra.h"
 #include "core/metric_kernels/abstract_metric.h"
 #include "core/table/dense_matrix.h"
 #include "core/table/memory_mapped_file.h"
@@ -70,9 +71,9 @@ class GenMetricTree {
       core::table::DensePoint col_point;
       for(int i = begin; i < end; i++) {
         matrix.MakeColumnVector(i, &col_point);
-        bounds->center() += col_point;
+        core::math::AddTo(col_point, &(bounds->center()));
       }
-      bounds->center() /= ((double) count);
+      core::math::Scale(1.0 / static_cast<double>(count), &(bounds->center()));
 
       double furthest_distance;
       FurthestColumnIndex_(
@@ -88,9 +89,11 @@ class GenMetricTree {
 
       // Compute the weighted sum of the two pivots
       node->bound().center().CopyValues(left->bound().center());
-      node->bound().center() *= left->count();
-      node->bound().center().Add(right->count(), right->bound().center());
-      node->bound().center() /= ((double) node->count());
+      core::math::Scale(left->count(), &(node->bound().center()));
+      core::math::AddExpert(
+        right->count(), right->bound().center(), & (node->bound().center()));
+      core::math::Scale(
+        1.0 / static_cast<double>(node->count()), & (node->bound().center()));
 
       double left_max_dist, right_max_dist;
       FurthestColumnIndex_(