Replace more instances of Eigen GEMV with Ceres BLAS.

With this ITERATIVE_SCHUR with JACOBI preconditioner went down from
280 seconds to 150 seconds on problem-744-543562-pre.txt.

Change-Id: I4f319c1108421e8d59f58654a4c0576ad65df609

CMakeLists.txt

@@ -654,9 +654,13 @@ IF ("${UNIX}")
       "${CMAKE_CXX_FLAGS} -Werror -Wall -Wextra -Wno-unknown-pragmas -Wno-sign-compare -Wno-unused-parameter")
 ENDIF ("${UNIX}")
 
-# We can be even stricter when using CLang
+# Use a larger inlining threshold for Clang, since it hobbles Eigen,
+# resulting in an unreasonably slow version of the blas routines. The
+# -Qunused-arguments is needed because CMake passes the inline
+# threshold to the linker and clang complains about it and dies.
 IF ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
-  SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-return-type-c-linkage")
+  SET(CMAKE_CXX_FLAGS
+      "${CMAKE_CXX_FLAGS} -Qunused-arguments -mllvm -inline-threshold=600 -Wno-return-type-c-linkage")
 ENDIF()
 
 ADD_SUBDIRECTORY(internal/ceres)

internal/ceres/block_sparse_matrix.cc

@@ -33,6 +33,7 @@
 #include <cstddef>
 #include <algorithm>
 #include <vector>
+#include "ceres/blas.h"
 #include "ceres/block_structure.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/matrix_proto.h"
@@ -117,16 +118,15 @@ void BlockSparseMatrix::RightMultiply(const double* x,  double* y) const {
   for (int i = 0; i < block_structure_->rows.size(); ++i) {
     int row_block_pos = block_structure_->rows[i].block.position;
     int row_block_size = block_structure_->rows[i].block.size;
-    VectorRef yref(y + row_block_pos, row_block_size);
     const vector<Cell>& cells = block_structure_->rows[i].cells;
     for (int j = 0; j < cells.size(); ++j) {
       int col_block_id = cells[j].block_id;
       int col_block_size = block_structure_->cols[col_block_id].size;
       int col_block_pos = block_structure_->cols[col_block_id].position;
-      ConstVectorRef xref(x + col_block_pos, col_block_size);
-      MatrixRef m(values_.get() + cells[j].position,
-                  row_block_size, col_block_size);
-      yref += m.lazyProduct(xref);
+      MatrixVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
+          values_.get() + cells[j].position, row_block_size, col_block_size,
+          x + col_block_pos,
+          y + row_block_pos);
     }
   }
 }
@@ -138,16 +138,16 @@ void BlockSparseMatrix::LeftMultiply(const double* x, double* y) const {
   for (int i = 0; i < block_structure_->rows.size(); ++i) {
     int row_block_pos = block_structure_->rows[i].block.position;
     int row_block_size = block_structure_->rows[i].block.size;
-    const ConstVectorRef xref(x + row_block_pos, row_block_size);
     const vector<Cell>& cells = block_structure_->rows[i].cells;
     for (int j = 0; j < cells.size(); ++j) {
       int col_block_id = cells[j].block_id;
       int col_block_size = block_structure_->cols[col_block_id].size;
       int col_block_pos = block_structure_->cols[col_block_id].position;
-      VectorRef yref(y + col_block_pos, col_block_size);
-      MatrixRef m(values_.get() + cells[j].position,
-                  row_block_size, col_block_size);
-      yref += m.transpose().lazyProduct(xref);
+      MatrixTransposeVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
+          values_.get() + cells[j].position, row_block_size, col_block_size,
+          x + row_block_pos,
+          y + col_block_pos);
+
     }
   }
 }

internal/ceres/partitioned_matrix_view.cc

@@ -35,6 +35,7 @@
 #include <algorithm>
 #include <cstring>
 #include <vector>
+#include "ceres/blas.h"
 #include "ceres/block_sparse_matrix.h"
 #include "ceres/block_structure.h"
 #include "ceres/internal/eigen.h"
@@ -103,13 +104,10 @@ void PartitionedMatrixView::RightMultiplyE(const double* x, double* y) const {
     const int col_block_id = cell.block_id;
     const int col_block_pos = bs->cols[col_block_id].position;
     const int col_block_size = bs->cols[col_block_id].size;
-
-    ConstVectorRef xref(x + col_block_pos, col_block_size);
-    VectorRef yref(y + row_block_pos, row_block_size);
-    ConstMatrixRef m(row_values + cell.position,
-                     row_block_size,
-                     col_block_size);
-    yref += m.lazyProduct(xref);
+    MatrixVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
+        row_values + cell.position, row_block_size, col_block_size,
+        x + col_block_pos,
+        y + row_block_pos);
   }
 }
 
@@ -124,20 +122,16 @@ void PartitionedMatrixView::RightMultiplyF(const double* x, double* y) const {
   for (int r = 0; r < bs->rows.size(); ++r) {
     const int row_block_pos = bs->rows[r].block.position;
     const int row_block_size = bs->rows[r].block.size;
-    VectorRef yref(y + row_block_pos, row_block_size);
     const vector<Cell>& cells = bs->rows[r].cells;
     for (int c = (r < num_row_blocks_e_) ? 1 : 0; c < cells.size(); ++c) {
       const double* row_values = matrix_.RowBlockValues(r);
       const int col_block_id = cells[c].block_id;
       const int col_block_pos = bs->cols[col_block_id].position;
       const int col_block_size = bs->cols[col_block_id].size;
-
-      ConstVectorRef xref(x + col_block_pos - num_cols_e(),
-                          col_block_size);
-      ConstMatrixRef m(row_values + cells[c].position,
-                       row_block_size,
-                       col_block_size);
-      yref += m.lazyProduct(xref);
+      MatrixVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
+          row_values + cells[c].position, row_block_size, col_block_size,
+          x + col_block_pos - num_cols_e(),
+          y + row_block_pos);
     }
   }
 }
@@ -155,13 +149,10 @@ void PartitionedMatrixView::LeftMultiplyE(const double* x, double* y) const {
     const int col_block_id = cell.block_id;
     const int col_block_pos = bs->cols[col_block_id].position;
     const int col_block_size = bs->cols[col_block_id].size;
-
-    ConstVectorRef xref(x + row_block_pos, row_block_size);
-    VectorRef yref(y + col_block_pos, col_block_size);
-    ConstMatrixRef m(row_values + cell.position,
-                     row_block_size,
-                     col_block_size);
-    yref += m.transpose().lazyProduct(xref);
+    MatrixTransposeVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
+        row_values + cell.position, row_block_size, col_block_size,
+        x + row_block_pos,
+        y + col_block_pos);
   }
 }
 
@@ -176,19 +167,16 @@ void PartitionedMatrixView::LeftMultiplyF(const double* x, double* y) const {
   for (int r = 0; r < bs->rows.size(); ++r) {
     const int row_block_pos = bs->rows[r].block.position;
     const int row_block_size = bs->rows[r].block.size;
-    ConstVectorRef xref(x + row_block_pos, row_block_size);
     const vector<Cell>& cells = bs->rows[r].cells;
     for (int c = (r < num_row_blocks_e_) ? 1 : 0; c < cells.size(); ++c) {
       const double* row_values = matrix_.RowBlockValues(r);
       const int col_block_id = cells[c].block_id;
       const int col_block_pos = bs->cols[col_block_id].position;
       const int col_block_size = bs->cols[col_block_id].size;
-
-      VectorRef yref(y + col_block_pos - num_cols_e(), col_block_size);
-      ConstMatrixRef m(row_values + cells[c].position,
-                       row_block_size,
-                       col_block_size);
-      yref += m.transpose().lazyProduct(xref);
+      MatrixTransposeVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
+        row_values + cells[c].position, row_block_size, col_block_size,
+        x + row_block_pos,
+        y + col_block_pos - num_cols_e());
     }
   }
 }
@@ -267,15 +255,15 @@ void PartitionedMatrixView::UpdateBlockDiagonalEtE(
     const int row_block_size = bs->rows[r].block.size;
     const int block_id = cell.block_id;
     const int col_block_size = bs->cols[block_id].size;
-    ConstMatrixRef m(row_values + cell.position,
-                           row_block_size,
-                           col_block_size);
-
     const int cell_position =
         block_diagonal_structure->rows[block_id].cells[0].position;
 
-    MatrixRef(block_diagonal->mutable_values() + cell_position,
-              col_block_size, col_block_size).noalias() += m.transpose() * m;
+    MatrixTransposeMatrixMultiply
+        <Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic, 1>(
+            row_values + cell.position, row_block_size, col_block_size,
+            row_values + cell.position, row_block_size, col_block_size,
+            block_diagonal->mutable_values() + cell_position,
+            0, 0, col_block_size, col_block_size);
   }
 }
 
@@ -298,15 +286,16 @@ void PartitionedMatrixView::UpdateBlockDiagonalFtF(
     for (int c = (r < num_row_blocks_e_) ? 1 : 0; c < cells.size(); ++c) {
       const int col_block_id = cells[c].block_id;
       const int col_block_size = bs->cols[col_block_id].size;
-      ConstMatrixRef m(row_values + cells[c].position,
-                       row_block_size,
-                       col_block_size);
       const int diagonal_block_id = col_block_id - num_col_blocks_e_;
       const int cell_position =
           block_diagonal_structure->rows[diagonal_block_id].cells[0].position;
 
-      MatrixRef(block_diagonal->mutable_values() + cell_position,
-                col_block_size, col_block_size).noalias() += m.transpose() * m;
+      MatrixTransposeMatrixMultiply
+          <Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic, 1>(
+              row_values + cells[c].position, row_block_size, col_block_size,
+              row_values + cells[c].position, row_block_size, col_block_size,
+              block_diagonal->mutable_values() + cell_position,
+              0, 0, col_block_size, col_block_size);
     }
   }
 }