commit e2a6cdc0816af9d0c77933f5017f137da3d52a35
author Keir Mierle <mierle@gmail.com> Mon May 07 06:39:56 2012 -0700
committer Keir Mierle <mierle@gmail.com> Mon May 07 06:39:56 2012 -0700
tree 45b85c8fe076dfb158d94aecb9af00d7edd1e9de
parent 1b95dc580aa5d89be021c0915e26df83f18013bb

Address some of the comments on CGNR patch

- Rename BlockDiagonalPreconditioner to BlockJacobiPreconditioner
- Include the diagonal in the block jacobi preconditioner.
- Better flag help for eta.
- Enable test for CGNR
- Rename CONJUGATE_GRADIENTS to CGNR.
- etc.

internal/ceres/cgnr_linear_operator.h

diff --git a/internal/ceres/cgnr_linear_operator.h b/internal/ceres/cgnr_linear_operator.h
index 94767fb..f32d8d9 100644
--- a/internal/ceres/cgnr_linear_operator.h
+++ b/internal/ceres/cgnr_linear_operator.h
@@ -61,12 +61,11 @@
 //
 // Then multiply out to get:
 //
-//   = xA^TAx - 2xA^T b + b^Tb + xD^TDx
 //   = xA^TAx - 2b^T Ax + b^Tb + xD^TDx
 //
 // Take the derivative:
 //
-//   0 = 2A^TAx - 2b^T A + b^Tb + 2 D^TDx
+//   0 = 2A^TAx - 2A^T b + 2 D^TDx
 //   0 = A^TAx - A^T b + D^TDx
 //   0 = (A^TA + D^TD)x - A^T b
 //
@@ -80,26 +79,25 @@
 // Note: This class is not thread safe, since it uses some temporary storage.
 class CgnrLinearOperator : public LinearOperator {
  public:
-  CgnrLinearOperator(LinearOperator* A, double *D)
-      : A_(A), D_(D), z_(new double[A->num_rows()]) {
+  CgnrLinearOperator(const LinearOperator& A, const double *D)
+      : A_(A), D_(D), z_(new double[A.num_rows()]) {
   }
   virtual ~CgnrLinearOperator() {}
 
   virtual void RightMultiply(const double* x, double* y) const {
-    std::fill(z_.get(), z_.get() + A_->num_rows(), 0.0);
+    std::fill(z_.get(), z_.get() + A_.num_rows(), 0.0);
 
     // z = Ax
-    A_->RightMultiply(x, z_.get());
+    A_.RightMultiply(x, z_.get());
 
     // y = y + Atz
-    A_->LeftMultiply(z_.get(), y);
+    A_.LeftMultiply(z_.get(), y);
 
     // y = y + DtDx
     if (D_ != NULL) {
-      int n = A_->num_cols();
-      for (int i = 0; i < n; ++i) {
-        y[i] += D_[i] * D_[i] * x[i];
-      }
+      int n = A_.num_cols();
+      VectorRef(y, n).array() += ConstVectorRef(D_, n).array().square() *
+                                 ConstVectorRef(x, n).array();
     }
   }
 
@@ -107,12 +105,12 @@
     RightMultiply(x, y);
   }
 
-  virtual int num_rows() const { return A_->num_cols(); }
-  virtual int num_cols() const { return A_->num_cols(); }
+  virtual int num_rows() const { return A_.num_cols(); }
+  virtual int num_cols() const { return A_.num_cols(); }
 
  private:
-  LinearOperator* A_;
-  double* D_;
+  const LinearOperator& A_;
+  const double* D_;
   scoped_array<double> z_;
 };