diff --git a/fastlib/trunk/contrib/dongryel/CMakeLists.txt b/fastlib/trunk/contrib/dongryel/CMakeLists.txt
index 3b0894a9ed..dc0ed53255 100644
--- a/fastlib/trunk/contrib/dongryel/CMakeLists.txt
+++ b/fastlib/trunk/contrib/dongryel/CMakeLists.txt
@@ -4,6 +4,7 @@ cmake_minimum_required(VERSION 2.8)
 set(DIRS
 #   compression
 #   fast_multipole_method  ## awaiting resolution of #10 (libint dependency)
+   gp_regression
 #   kde  ## does not compile
    linear_algebra
    linear_regression  
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/CMakeLists.txt b/fastlib/trunk/contrib/dongryel/gp_regression/CMakeLists.txt
new file mode 100644
index 0000000000..de25c6b0e1
--- /dev/null
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/CMakeLists.txt
@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.8)
+
+# test executable
+add_executable(bilinear_form_test
+               EXCLUDE_FROM_ALL
+               bilinear_form_test.cc
+              )
+# link dependencies of test executable
+target_link_libraries(bilinear_form_test
+                      fastlib
+                     )
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/bilinear_form_estimator.h b/fastlib/trunk/contrib/dongryel/gp_regression/bilinear_form_estimator.h
new file mode 100644
index 0000000000..e95414a91c
--- /dev/null
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/bilinear_form_estimator.h
@@ -0,0 +1,128 @@
+/** @author Dongryeol Lee
+ *
+ *  @file bilinear_form_estimator.h
+ */
+
+#ifndef FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_BILINEAR_FORM_ESTIMATOR_H
+#define FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_BILINEAR_FORM_ESTIMATOR_H
+
+#include <vector>
+#include "fastlib/la/matrix.h"
+#include "fastlib/la/uselapack.h"
+#include "linear_operator.h"
+
+namespace fl {
+namespace ml {
+
+class SquareRootTransformation {
+  public:
+    static double Transform(double val_in) {
+      return sqrt(val_in);
+    }
+};
+
+class IdentityTransformation {
+  public:
+    static double Transform(double val_in) {
+      return val_in;
+    }
+};
+
+class InverseTransformation {
+  public:
+    static double Transform(double val_in) {
+      return 1.0 / val_in;
+    }
+};
+
+class LogTransformation {
+  public:
+    static double Transform(double val_in) {
+      return log(val_in);
+    }
+};
+
+template<typename TransformationType>
+class BilinearFormEstimator {
+
+  private:
+
+    class TridiagonalLinearOperator {
+      private:
+        const std::vector<double> *diagonal_entries_;
+        const std::vector<double> *offdiagonal_entries_;
+
+      public:
+
+        int n_rows() const;
+
+        int n_cols() const;
+
+        double get(int row, int col) const;
+
+        const std::vector<double> *diagonal_entries() const;
+
+        TridiagonalLinearOperator(
+          const std::vector<double> &diagonal_entries_in,
+          const std::vector<double> &offdiagonal_entries_in);
+
+        int Apply(const Epetra_MultiVector &vecs,
+                  Epetra_MultiVector &prods) const;
+
+        void PrintDebug(const char *name = "", FILE *stream = stderr) const;
+    };
+
+  private:
+
+#ifdef EPETRA_MPI
+    Epetra_MpiComm comm_;
+#else
+    Epetra_SerialComm comm_;
+#endif
+
+    Anasazi::LinearOperator *linear_operator_;
+
+    const Epetra_Map *map_;
+
+  private:
+    void AddExpert_(double scalar,
+                    const Epetra_MultiVector &source,
+                    Epetra_MultiVector *destination) const;
+
+    void Scale_(double scalar,
+                const Epetra_MultiVector &source,
+                Epetra_MultiVector *destination) const;
+
+    double Dot_(const Epetra_MultiVector &first_vec,
+                const Epetra_MultiVector &second_vec) const;
+
+    template<typename LinearOperatorType>
+    double ComputeQuadraticForm_(
+      int num_iterations,
+      const GenVector<double> &starting_vector,
+      const LinearOperatorType &linear_operator_in,
+      const Epetra_Map &map_in,
+      int level,
+      bool *broke_down);
+
+  public:
+
+    BilinearFormEstimator();
+
+    Anasazi::LinearOperator *linear_operator();
+
+    void Init(Anasazi::LinearOperator *linear_operator_in);
+
+    double Compute(
+      const GenVector<double> &left_argument,
+      const GenVector<double> &right_argument,
+      bool naive_compute);
+
+    double Compute(const GenVector<double> &argument);
+
+    double NaiveCompute(const GenVector<double> &argument);
+};
+};
+};
+
+#endif
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/bilinear_form_estimator_dev.h b/fastlib/trunk/contrib/dongryel/gp_regression/bilinear_form_estimator_dev.h
new file mode 100644
index 0000000000..e2b4169fc6
--- /dev/null
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/bilinear_form_estimator_dev.h
@@ -0,0 +1,356 @@
+/** @author Dongryeol Lee
+ *
+ *  @file bilinear_form_estimator_dev.h
+ */
+
+#ifndef FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_BILINEAR_FORM_ESTIMATOR_DEV_H
+#define FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_BILINEAR_FORM_ESTIMATOR_DEV_H
+
+#include <vector>
+#include "fastlib/la/matrix.h"
+#include "fastlib/la/la.h"
+#include "fastlib/la/uselapack.h"
+#include "bilinear_form_estimator.h"
+
+namespace fl {
+namespace ml {
+
+template<typename TransformationType>
+void BilinearFormEstimator<TransformationType>::
+TridiagonalLinearOperator::PrintDebug(
+  const char *name, FILE *stream) const {
+  fprintf(stream, "----- MATRIX ------: %s\n", name);
+  for (int r = 0; r < this->n_rows(); r++) {
+    for (int c = 0; c < this->n_cols(); c++) {
+      fprintf(stream, "%+3.3f ", this->get(r, c));
+    }
+    fprintf(stream, "\n");
+  }
+}
+
+template<typename TransformationType>
+int BilinearFormEstimator<TransformationType>::
+TridiagonalLinearOperator::n_rows() const {
+  return diagonal_entries_->size();
+}
+
+template<typename TransformationType>
+int BilinearFormEstimator<TransformationType>::
+TridiagonalLinearOperator::n_cols() const {
+  return diagonal_entries_->size();
+}
+
+template<typename TransformationType>
+double BilinearFormEstimator<TransformationType>::
+TridiagonalLinearOperator::get(int row, int col) const {
+  if (row == col) {
+    return (*diagonal_entries_)[row];
+  }
+  else if (row == col + 1 || col == row + 1) {
+    return (*offdiagonal_entries_)[ std::min(row, col)];
+  }
+  else {
+    return 0;
+  }
+}
+
+template<typename TransformationType>
+const std::vector<double> *BilinearFormEstimator<TransformationType>::
+TridiagonalLinearOperator::diagonal_entries() const {
+
+  return diagonal_entries_;
+}
+
+template<typename TransformationType>
+BilinearFormEstimator<TransformationType>::
+TridiagonalLinearOperator::TridiagonalLinearOperator(
+  const std::vector<double> &diagonal_entries_in,
+  const std::vector<double> &offdiagonal_entries_in) {
+
+  diagonal_entries_ = &diagonal_entries_in;
+  offdiagonal_entries_ = &offdiagonal_entries_in;
+}
+
+template<typename TransformationType>
+int BilinearFormEstimator<TransformationType>::
+TridiagonalLinearOperator::Apply(
+  const Epetra_MultiVector &vecs,
+  Epetra_MultiVector &prods) const {
+
+  prods.PutScalar(0);
+
+  for (int j = 0; j < diagonal_entries_->size(); j++) {
+
+    for (int k = 0; k < vecs.NumVectors(); k++) {
+
+      // Apply the diagonal entry.
+      prods.Pointers()[k][j] += ((*diagonal_entries_)[j]) *
+                                vecs.Pointers()[k][j];
+
+      // Apply the lower diagonal entry.
+      if (j > 0) {
+        prods.Pointers()[k][j - 1] += ((*offdiagonal_entries_)[j - 1]) *
+                                      vecs.Pointers()[k][j - 1];
+      }
+
+      // Apply the upper diagonal entry.
+      if (j < diagonal_entries_->size() - 1) {
+        prods.Pointers()[k][j + 1] += ((*offdiagonal_entries_)[j]) *
+                                      vecs.Pointers()[k][j + 1];
+      }
+    }
+  }
+  return 0;
+}
+
+#ifdef EPETRA_MPI
+template<typename TransformationType>
+BilinearFormEstimator<TransformationType>::BilinearFormEstimator(): comm_(MPI_COMM_WORLD) {
+
+}
+#else
+template<typename TransformationType>
+BilinearFormEstimator<TransformationType>::BilinearFormEstimator() {
+
+}
+#endif
+
+template<typename TransformationType>
+Anasazi::LinearOperator *BilinearFormEstimator<TransformationType>::
+linear_operator() {
+  return linear_operator_;
+}
+
+template<typename TransformationType>
+void BilinearFormEstimator<TransformationType>::Init(
+  Anasazi::LinearOperator *linear_operator_in) {
+
+  linear_operator_ = linear_operator_in;
+  map_ = &(linear_operator_->OperatorDomainMap());
+}
+
+template<typename TransformationType>
+double BilinearFormEstimator<TransformationType>::Dot_(
+  const Epetra_MultiVector &first_vec,
+  const Epetra_MultiVector &second_vec) const {
+
+  double dot_product = 0;
+  for (int i = 0; i < first_vec.GlobalLength(); i++) {
+    dot_product += first_vec.Pointers()[0][i] * second_vec.Pointers()[0][i];
+  }
+  return dot_product;
+}
+
+template<typename TransformationType>
+void BilinearFormEstimator<TransformationType>::AddExpert_(double scalar,
+    const Epetra_MultiVector &source,
+    Epetra_MultiVector *destination) const {
+  for (int i = 0; i < source.GlobalLength(); i++) {
+    destination->Pointers()[0][i] =
+      destination->Pointers()[0][i] + scalar * source.Pointers()[0][i];
+  }
+}
+
+template<typename TransformationType>
+void BilinearFormEstimator<TransformationType>::Scale_(double scalar,
+    const Epetra_MultiVector &source,
+    Epetra_MultiVector *destination) const {
+  for (int i = 0; i < source.GlobalLength(); i++) {
+    destination->Pointers()[0][i] = scalar * source.Pointers()[0][i];
+  }
+}
+
+template<typename TransformationType>
+template<typename LinearOperatorType>
+double BilinearFormEstimator<TransformationType>::ComputeQuadraticForm_(
+  int num_iterations,
+  const GenVector<double> &starting_vector,
+  const LinearOperatorType &linear_operator_in,
+  const Epetra_Map &map_in,
+  int level,
+  bool *break_down) {
+
+  // The threshold for determining the convergence.
+  const double convergence_threshold = 1e-7;
+
+  // If it is a 1 by 1 matrix, then apply the transformation.
+  if (linear_operator_in.n_rows() == 1 && linear_operator_in.n_cols() == 1) {
+    return TransformationType::Transform(linear_operator_in.get(0, 0));
+  }
+
+  // The diagonal entries and the offdiagonal entries with the wrapper
+  // class around it.
+  std::vector<double> diagonal_entries;
+  std::vector<double> offdiagonal_entries;
+  TridiagonalLinearOperator tridiagonal_linear_operator(
+    diagonal_entries, offdiagonal_entries);
+
+  // The basis vector in the previous iteration.
+  Epetra_MultiVector previous_vector(map_in, 1);
+  previous_vector.PutScalar(0);
+
+  // The basis vector in the current iteration.
+  Epetra_MultiVector current_vector(map_in, 1);
+  for (int i = 0; i < starting_vector.length(); i++) {
+    current_vector.Pointers()[0][i] = starting_vector[i];
+  }
+
+  // A temporary vector used for matrix-vector multiplication.
+  Epetra_MultiVector residual_vector(map_in, 1);
+  residual_vector.PutScalar(0);
+
+  // A temporary vector used for denoting the unit vector of varying
+  // dimension.
+  GenVector<double> unit_vector;
+  unit_vector.Init(starting_vector.length());
+  unit_vector.SetZero();
+  unit_vector[0] = 1;
+
+  // The old bilinear estimate.
+  double old_bilinear_estimate = std::numeric_limits<double>::max();
+
+  for (int j = 0; j < num_iterations; j++) {
+
+    linear_operator_in.Apply(current_vector, residual_vector);
+    if (j > 0) {
+      AddExpert_(- (offdiagonal_entries[j - 1]),
+                 previous_vector, &residual_vector);
+    }
+
+    // The dot product with the residual and the current basis
+    // vector. Compute the off-diagonal and the diagonal entries
+    // in this iteration.
+    double alpha_j = Dot_(residual_vector, current_vector);
+    diagonal_entries.push_back(alpha_j);
+    AddExpert_(- alpha_j, current_vector, &residual_vector);
+    double beta_j_plus_one =
+      sqrt(Dot_(residual_vector, residual_vector));
+
+    // Add in the offdiagonal entry.
+    if (fabs(beta_j_plus_one) <= convergence_threshold) {
+      *break_down = true;
+      break;
+    }
+    offdiagonal_entries.push_back(beta_j_plus_one);
+
+    // Take the current tridiagonal decomposition and estimate the
+    // bilinear form. It is essential that Epetra_Map is constructed
+    // here right before the recursive call.
+    GenVector<double> unit_vector_alias;
+    unit_vector_alias.Alias(unit_vector.ptr(),
+                            diagonal_entries.size());
+    Epetra_Map tridiagonal_linear_operator_map(
+      unit_vector_alias.length(), 0, comm_);
+
+    bool subcase_break_down = false;
+    double bilinear_estimate =
+      ComputeQuadraticForm_(unit_vector_alias.length(), unit_vector_alias,
+                            tridiagonal_linear_operator,
+                            tridiagonal_linear_operator_map, level + 1,
+                            &subcase_break_down);
+
+    // Check whether it converged.
+    if (subcase_break_down) {
+      *break_down = true;
+      break;
+    }
+    if (fabs(old_bilinear_estimate - bilinear_estimate) <=
+        convergence_threshold && false) {
+      *break_down = false;
+      break;
+    }
+    else {
+      old_bilinear_estimate = bilinear_estimate;
+    }
+
+    // Update the previous vector and the next vector.
+    for (int i = 0; i < current_vector.GlobalLength(); i++) {
+      previous_vector.Pointers()[0][i] = current_vector.Pointers()[0][i];
+    }
+    Scale_(1.0 / beta_j_plus_one, residual_vector, &current_vector);
+  }
+
+  return old_bilinear_estimate;
+}
+
+template<typename TransformationType>
+double BilinearFormEstimator<TransformationType>::Compute(
+  const GenVector<double> &left_argument,
+  const GenVector<double> &right_argument,
+  bool naive_compute) {
+
+  // Use the formula: $u^T f(A) v = 0.25 * (y^T f(A) y - z^T f(A) z )
+  // $ where $y = u + v$ and $z = u - v$.
+  GenVector<double> sum, difference;
+  la::AddInit(
+    left_argument, right_argument, &sum);
+  la::SubInit(
+    right_argument, left_argument, &difference);
+  double estimate =
+    (naive_compute) ?
+    0.25 * (NaiveCompute(sum) - NaiveCompute(difference)) :
+    0.25 * (Compute(sum) - Compute(difference));
+  return estimate;
+}
+
+template<typename TransformationType>
+double BilinearFormEstimator<TransformationType>::NaiveCompute(
+  const GenVector<double> &argument) {
+
+  // The naive estimate to return.
+  double naive_bilinear_estimate = 0;
+
+  // Compute the kernel matrix, and its eigendecomposition.
+  GenMatrix<double> kernel_matrix;
+  GenVector<double> eigenvalues;
+  GenMatrix<double> eigenvectors;
+  GenMatrix<double> eigenvectors_transposed;
+  kernel_matrix.Init(linear_operator_->n_rows(), linear_operator_->n_cols());
+  for (int j = 0; j < linear_operator_->n_cols(); j++) {
+    for (int i = 0; i < linear_operator_->n_rows(); i++) {
+      kernel_matrix.set(i, j, linear_operator_->get(i, j));
+    }
+  }
+  la::SVDInit(
+    kernel_matrix, &eigenvalues, &eigenvectors, &eigenvectors_transposed);
+  eigenvalues.PrintDebug();
+
+  // Project the argument to the eigenspace.
+  GenVector<double> projected_vector;
+  la::MulInit(
+    eigenvectors_transposed, argument, &projected_vector);
+  for (int i = 0; i < eigenvalues.length(); i++) {
+    naive_bilinear_estimate += math::Sqr(projected_vector[i]) *
+                               TransformationType::Transform(eigenvalues[i]);
+  }
+  return naive_bilinear_estimate;
+}
+
+template<typename TransformationType>
+double BilinearFormEstimator<TransformationType>::Compute(
+  const GenVector<double> &argument) {
+
+  // Pass in the normalized unit vector to the quadratic form
+  // computation and correct it afterwards.
+  GenVector<double> normalized_argument;
+  double length = la::LengthEuclidean(argument);
+  if (length > 0) {
+    la::ScaleInit(
+      1.0 / length, argument, &normalized_argument);
+  }
+  else {
+    normalized_argument.Copy(argument);
+  }
+  bool break_down = false;
+  return math::Sqr(length) * ComputeQuadraticForm_(
+           normalized_argument.length(),
+           normalized_argument,
+           *linear_operator_,
+           *map_,
+           0,
+           &break_down);
+}
+};
+};
+
+#endif
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/bilinear_form_test.cc b/fastlib/trunk/contrib/dongryel/gp_regression/bilinear_form_test.cc
new file mode 100644
index 0000000000..8505141ea9
--- /dev/null
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/bilinear_form_test.cc
@@ -0,0 +1,84 @@
+/** @author Dongryeol Lee
+ *
+ *  @file bilinear_form.test.cc
+ */
+
+#undef BOOST_ALL_DYN_LINK
+#include "fastlib/fastlib.h"
+#include "boost/program_options.hpp"
+#include "boost/test/included/unit_test.hpp"
+#include "boost/mpl/map.hpp"
+#include "boost/mpl/if.hpp"
+#include "bilinear_form_estimator_dev.h"
+#include "log_determinant_dev.h"
+
+#ifdef EPETRA_MPI
+#include "trilinos/Epetra_MpiComm.h"
+#else
+#include "trilinos/Epetra_SerialComm.h"
+#endif
+
+namespace fl {
+namespace ml {
+namespace bilinear_form_test {
+class BilinearFormTestSuite : public boost::unit_test_framework::test_suite {
+
+  public:
+
+    class BilinearFormTest {
+      public:
+
+        BilinearFormTest() {
+        }
+
+        void RunTests() {
+
+	  fprintf(stderr, "Running the tests:\n");
+
+          // Call MPI Finalize.
+          MPI_Finalize();
+        }
+    };
+
+  public:
+
+    BilinearFormTestSuite()
+        : boost::unit_test_framework::test_suite("Bilinear form test suite") {
+
+      // Create an instance of test.
+      boost::shared_ptr<BilinearFormTest> instance(new BilinearFormTest());
+
+      // Create the test cases.
+      boost::unit_test_framework::test_case* bilinear_form_test_case
+      = BOOST_CLASS_TEST_CASE(
+          &BilinearFormTest::RunTests, instance);
+      // add the test cases to the test suite
+      add(bilinear_form_test_case);
+    }
+};
+};
+};
+};
+
+boost::unit_test_framework::test_suite*
+init_unit_test_suite(int argc, char** argv) {
+
+  // Initialize MPI.
+#ifdef EPETRA_MPI
+  MPI_Init(&argc, &argv);
+#endif
+
+  // create the top test suite
+  boost::unit_test_framework::test_suite* top_test_suite
+  = BOOST_TEST_SUITE("Bilinear form tests");
+
+  if (argc != 2) {
+    NOTIFY("Wrong number of arguments for tree test. Expected test input files directory. Returning NULL.");
+    return NULL;
+  }
+
+  // add test suites to the top test suite
+  std::string input_files_directory = argv[1];
+  top_test_suite->add(new fl::ml::bilinear_form_test::BilinearFormTestSuite());
+  return top_test_suite;
+}
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/kernel_linear_operator.h b/fastlib/trunk/contrib/dongryel/gp_regression/kernel_linear_operator.h
new file mode 100644
index 0000000000..793aa4e777
--- /dev/null
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/kernel_linear_operator.h
@@ -0,0 +1,166 @@
+#ifndef FASTLIB_CONTRIB_DONGRYEL_TRILINOS_WRAPPERS_KERNEL_LINEAR_OPERATOR_H
+#define FASTLIB_CONTRIB_DONGRYEL_TRILINOS_WRAPPERS_KERNEL_LINEAR_OPERATOR_H
+
+#include "fastlib/math/fl_math.h"
+#include "fastlib/trilinos_wrappers/linear_operator.h"
+
+namespace Anasazi {
+
+template<bool dotproduct_selfcase_special>
+class DotProductTrait {
+  public:
+    template<typename KernelType, typename PointType>
+    DotProductTrait(const KernelType *kernel,
+                    const PointType &first_point,
+                    const PointType &second_point,
+                    bool point_indices_are_same,
+                    double *dotproduct);
+};
+
+template<>
+class DotProductTrait<true> {
+  public:
+    template<typename KernelType, typename PointType>
+    DotProductTrait(const KernelType *kernel,
+                    const PointType &first_point,
+                    const PointType &second_point,
+                    bool point_indices_are_same,
+                    double *dotproduct) {
+      *dotproduct = kernel->Dot(first_point, second_point,
+                                point_indices_are_same);
+    }
+};
+
+template<>
+class DotProductTrait<false> {
+  public:
+    template<typename KernelType, typename PointType>
+    DotProductTrait(const KernelType *kernel,
+                    const PointType &first_point,
+                    const PointType &second_point,
+                    bool point_indices_are_same,
+                    double *dotproduct) {
+      *dotproduct = kernel->Dot(first_point, second_point);
+    }
+};
+
+template < typename TableType, typename KernelType, bool do_centering,
+bool dotproduct_selfcase_special = false >
+class KernelLinearOperator: public LinearOperator {
+
+  private:
+
+    TableType *table_;
+
+    const KernelType *kernel_;
+
+    fl::data::MonolithicPoint<double> average_row_;
+
+    double average_;
+
+  public:
+
+    KernelLinearOperator(TableType &table_in,
+                         const KernelType &kernel_in,
+#ifdef EPETRA_MPI
+                         const Epetra_MpiComm &comm_in,
+#else
+                         const Epetra_SerialComm &comm_in,
+#endif
+                         const Epetra_Map &map_in) {
+
+      table_ = &table_in;
+      kernel_ = &kernel_in;
+      comm_ = &comm_in;
+      map_ = &map_in;
+
+      if (do_centering) {
+        average_row_.Init(table_in.n_entries());
+      }
+      average_ = 0;
+
+      if (do_centering) {
+
+        // Precompute the average. This is a naive way of computing it.
+        for (int i = 0; i < table_in.n_entries(); i++) {
+          double average_for_i_th_point = 0;
+          for (int j = 0; j < table_in.n_entries(); j++) {
+            average_for_i_th_point += kernel_value(i, j);
+          }
+          average_for_i_th_point /= ((double) table_in.n_entries());
+          average_row_[i] = average_for_i_th_point;
+        }
+        for (int i = 0; i < table_in.n_entries(); i++) {
+          average_ += average_row_[i];
+        }
+        average_ /= ((double) table_in.n_entries());
+      }
+    }
+
+    double centered_kernel_value(int row, int col) const {
+      typename TableType::Dataset_t::Point_t row_point;
+      typename TableType::Dataset_t::Point_t col_point;
+      table_->get(row, &row_point);
+      table_->get(col, &col_point);
+      double dotproduct = 0;
+      DotProductTrait<dotproduct_selfcase_special>(
+        kernel_, row_point, col_point, row == col, &dotproduct);
+      return dotproduct - average_row_[row] - average_row_[col] + average_;
+    }
+
+    double kernel_value(int row, int col) const {
+      typename TableType::Dataset_t::Point_t row_point;
+      typename TableType::Dataset_t::Point_t col_point;
+      table_->get(row, &row_point);
+      table_->get(col, &col_point);
+      double dotproduct = 0;
+      DotProductTrait<dotproduct_selfcase_special>(
+        kernel_, row_point, col_point, row == col, &dotproduct);
+      return dotproduct;
+    }
+
+    int Apply(
+      const Epetra_MultiVector &vecs,
+      Epetra_MultiVector &prods) const {
+
+      prods.PutScalar(0);
+
+      for (int j = 0; j < table_->n_entries(); j++) {
+        for (int i = 0; i < table_->n_entries(); i++) {
+          double pair_kernel_value = (do_centering) ?
+                                     centered_kernel_value(i, j) : kernel_value(i, j);
+          for (int k = 0; k < vecs.NumVectors(); k++) {
+            prods.Pointers()[k][i] += pair_kernel_value * vecs.Pointers()[k][j];
+          }
+        }
+      }
+      return 0;
+    }
+
+    int n_rows() const {
+      return table_->n_entries();
+    }
+
+    int n_cols() const {
+      return table_->n_entries();
+    }
+
+    double get(int row, int col) const {
+      return (do_centering) ? centered_kernel_value(row, col) : kernel_value(row, col);
+    }
+};
+
+template<typename TableType, typename KernelType, bool do_centering>
+class OperatorTraits < double, Epetra_MultiVector,
+      KernelLinearOperator<TableType, KernelType, do_centering> > {
+  public:
+
+    static void Apply(const Epetra_Operator& Op,
+                      const Epetra_MultiVector& x,
+                      Epetra_MultiVector& y) {
+      Op.Apply(x, y);
+    }
+};
+};
+
+#endif
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/linear_operator.h b/fastlib/trunk/contrib/dongryel/gp_regression/linear_operator.h
new file mode 100644
index 0000000000..b4c792062c
--- /dev/null
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/linear_operator.h
@@ -0,0 +1,129 @@
+
+#ifndef FASTLIB_CONTRIB_DONGRYEL_FASTLIB_TRILINOS_WRAPPERS_LINEAR_OPERATOR_H
+#define FASTLIB_CONTRIB_DONGRYEL_FASTLIB_TRILINOS_WRAPPERS_LINEAR_OPERATOR_H
+
+#undef F77_FUNC
+#undef LI
+#include "trilinos/AnasaziEpetraAdapter.hpp"
+#include "trilinos/AnasaziBasicEigenproblem.hpp"
+#include "trilinos/AnasaziBlockKrylovSchurSolMgr.hpp"
+#include "trilinos/AnasaziBasicSort.hpp"
+#include "trilinos/AztecOO.h"
+#include "trilinos/Epetra_BlockMap.h"
+#include "trilinos/Epetra_CrsMatrix.h"
+#include "trilinos/Epetra_DataAccess.h"
+#include "trilinos/Epetra_LinearProblem.h"
+#include "trilinos/Epetra_Map.h"
+#include "trilinos/Epetra_MultiVector.h"
+#include "trilinos/Epetra_Operator.h"
+
+#ifdef EPETRA_MPI
+#include "trilinos/Epetra_MpiComm.h"
+#else
+#include "trilinos/Epetra_SerialComm.h"
+#endif
+
+#include "trilinos/Epetra_Vector.h"
+#undef F77_FUNC
+#include "fastlib/la/matrix.h"
+#include <vector>
+
+namespace Anasazi {
+
+class LinearOperator: public virtual Epetra_Operator {
+
+  protected:
+
+#ifdef EPETRA_MPI
+    const Epetra_MpiComm *comm_;
+#else
+    const Epetra_SerialComm *comm_;
+#endif
+
+    const Epetra_Map *map_;
+
+  public:
+
+    virtual ~LinearOperator() {
+    }
+
+    LinearOperator() {
+      comm_ = NULL;
+      map_ = NULL;
+    }
+
+#ifdef EPETRA_MPI
+    LinearOperator(const Epetra_MpiComm &comm_in,
+                   const Epetra_Map &map_in) {
+      comm_ = &comm_in;
+      map_ = &map_in;
+    }
+#else
+    LinearOperator(const Epetra_SerialComm &comm_in,
+                   const Epetra_Map &map_in) {
+      comm_ = &comm_in;
+      map_ = &map_in;
+    }
+#endif
+
+    virtual int Apply(const Epetra_MultiVector &vec,
+                      Epetra_MultiVector &prod) const = 0;
+
+    int SetUseTranspose(bool use_transpose) {
+      return -1;
+    }
+
+    int ApplyInverse(const Epetra_MultiVector &X,
+                     Epetra_MultiVector &Y) const {
+      return -1;
+    }
+
+    double NormInf() const {
+      return -1;
+    }
+
+    const char *Label() const {
+      return "Generic linear operator";
+    }
+
+    bool UseTranspose() const {
+      return false;
+    }
+
+    bool HasNormInf() const {
+      return false;
+    }
+
+    const Epetra_Comm &Comm() const {
+      return *comm_;
+    }
+
+    const Epetra_Map &OperatorDomainMap() const {
+      const Epetra_Map &map_reference = *map_;
+      return map_reference;
+    }
+
+    const Epetra_Map &OperatorRangeMap() const {
+      const Epetra_Map &map_reference = *map_;
+      return map_reference;
+    }
+
+    void PrintDebug(const char *name = "", FILE *stream = stderr) const {
+      fprintf(stream, "----- MATRIX ------: %s\n", name);
+      for (int r = 0; r < this->n_rows(); r++) {
+        for (int c = 0; c < this->n_cols(); c++) {
+          fprintf(stream, "%+3.3f ", this->get(r, c));
+        }
+        fprintf(stream, "\n");
+      }
+    }
+
+    virtual int n_rows() const = 0;
+
+    virtual int n_cols() const = 0;
+
+    virtual double get(int row, int col) const = 0;
+};
+};
+
+#endif
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/log_determinant.h b/fastlib/trunk/contrib/dongryel/gp_regression/log_determinant.h
new file mode 100644
index 0000000000..7755a1f6fb
--- /dev/null
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/log_determinant.h
@@ -0,0 +1,39 @@
+/** @author Dongryeol Lee
+ *
+ *  @file log_determinant.h
+ */
+
+#ifndef FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_LOG_DETERMINANT_H
+#define FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_LOG_DETERMINANT_H
+
+#include "bilinear_form_estimator.h"
+#include "fastlib/la/matrix.h"
+
+namespace fl {
+namespace ml {
+class LogDeterminant {
+
+  private:
+
+    fl::ml::BilinearFormEstimator<fl::ml::LogTransformation> bilinear_log_form_;
+
+  private:
+
+    void RandomVector_(GenVector<double> &v);
+
+  public:
+
+    LogDeterminant();
+
+    void Init(Anasazi::LinearOperator *linear_operator_in);
+
+    double MonteCarloCompute();
+
+    double Compute();
+
+    double NaiveCompute();
+};
+};
+};
+
+#endif
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/log_determinant_dev.h b/fastlib/trunk/contrib/dongryel/gp_regression/log_determinant_dev.h
new file mode 100644
index 0000000000..db386efeec
--- /dev/null
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/log_determinant_dev.h
@@ -0,0 +1,85 @@
+/** @author Dongryeol Lee
+ *
+ *  @file log_determinant_dev.h
+ */
+
+#ifndef FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_LOG_DETERMINANT_DEV_H
+#define FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_LOG_DETERMINANT_DEV_H
+
+#include "fastlib/la/matrix.h"
+#include "bilinear_form_estimator_dev.h"
+#include "log_determinant.h"
+
+namespace fl {
+namespace ml {
+LogDeterminant::LogDeterminant() {
+}
+
+void LogDeterminant::Init(Anasazi::LinearOperator *linear_operator_in) {
+  bilinear_log_form_.Init(linear_operator_in);
+}
+
+void LogDeterminant::RandomVector_(GenVector<double> &v) {
+  for (int i = 0; i < v.length(); i++) {
+    v[i] = (math::Random() >= 0.5) ? 1 : -1;
+  }
+}
+
+double LogDeterminant::MonteCarloCompute() {
+
+  // A random vector for the samples.
+  GenVector<double> random_vector;
+  random_vector.Init(bilinear_log_form_.linear_operator()->n_rows());
+
+  double log_determinant = 0;
+  const int num_samples = 100;
+
+  for (int i = 0; i < num_samples; i++) {
+    RandomVector_(random_vector);
+    log_determinant += bilinear_log_form_.Compute(random_vector);
+  }
+  log_determinant /= ((double) num_samples);
+  return log_determinant;
+}
+
+double LogDeterminant::NaiveCompute() {
+
+  double log_determinant = 0;
+  GenMatrix<double> kernel_matrix;
+  kernel_matrix.Init(bilinear_log_form_.linear_operator()->n_rows(),
+                     bilinear_log_form_.linear_operator()->n_cols());
+  for (int j = 0; j < bilinear_log_form_.linear_operator()->n_cols(); j++) {
+    for (int i = 0; i < bilinear_log_form_.linear_operator()->n_rows(); i++) {
+      kernel_matrix.set(i, j, bilinear_log_form_.linear_operator()->get(i, j));
+    }
+  }
+  GenVector<double> eigenvalues;
+  la::SVDInit(kernel_matrix, &eigenvalues);
+
+  for (int i = 0; i < eigenvalues.length(); i++) {
+    log_determinant += log(eigenvalues[i]);
+  }
+  return log_determinant;
+}
+
+double LogDeterminant::Compute() {
+
+  // Do a naive for-loop over each row and apply $e_i^T log(A) e_i$.
+  GenVector<double> i_th_unit_vector;
+  i_th_unit_vector.Init(bilinear_log_form_.linear_operator()->n_rows());
+  i_th_unit_vector.SetZero();
+  double log_determinant = 0;
+
+  for (int i = 0; i < bilinear_log_form_.linear_operator()->n_rows(); i++) {
+    i_th_unit_vector[i] = 1;
+    if (i > 0) {
+      i_th_unit_vector[i - 1] = 0;
+    }
+    log_determinant += bilinear_log_form_.Compute(i_th_unit_vector);
+  }
+  return log_determinant;
+}
+};
+};
+
+#endif