← Back to branch summary

~ubuntu-branches/ubuntu/utopic/blender/utopic-proposed

« back to all changes in this revision

Viewing changes to extern/libmv/third_party/ceres/internal/ceres/dense_qr_solver.cc

  • Committer: Package Import Robot
  • Author(s): Matthias Klose
  • Date: 2014-02-19 11:24:23 UTC
  • mfrom: (14.2.23 sid)
  • Revision ID: package-import@ubuntu.com-20140219112423-rkmaz2m7ha06d4tk
Tags: 2.69-3ubuntu1
* Merge with Debian; remaining changes:
  - Configure without OpenImageIO on armhf, as it is not available on
    Ubuntu.

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/libmv/third_party/ceres/internal/ceres/dense_qr_solver.cc
30
30
 
31
31
#include "ceres/dense_qr_solver.h"
32
32
 
 
33
 
33
34
#include <cstddef>
34
 
 
35
35
#include "Eigen/Dense"
36
36
#include "ceres/dense_sparse_matrix.h"
37
37
#include "ceres/internal/eigen.h"
38
38
#include "ceres/internal/scoped_ptr.h"
 
39
#include "ceres/lapack.h"
39
40
#include "ceres/linear_solver.h"
40
41
#include "ceres/types.h"
41
42
#include "ceres/wall_time.h"
44
45
namespace internal {
45
46
 
46
47
DenseQRSolver::DenseQRSolver(const LinearSolver::Options& options)
47
 
    : options_(options) {}
 
48
    : options_(options) {
 
49
  work_.resize(1);
 
50
}
48
51
 
49
52
LinearSolver::Summary DenseQRSolver::SolveImpl(
50
53
    DenseSparseMatrix* A,
51
54
    const double* b,
52
55
    const LinearSolver::PerSolveOptions& per_solve_options,
53
56
    double* x) {
 
57
  if (options_.dense_linear_algebra_library_type == EIGEN) {
 
58
    return SolveUsingEigen(A, b, per_solve_options, x);
 
59
  } else {
 
60
    return SolveUsingLAPACK(A, b, per_solve_options, x);
 
61
  }
 
62
}
 
63
LinearSolver::Summary DenseQRSolver::SolveUsingLAPACK(
 
64
    DenseSparseMatrix* A,
 
65
    const double* b,
 
66
    const LinearSolver::PerSolveOptions& per_solve_options,
 
67
    double* x) {
 
68
  EventLogger event_logger("DenseQRSolver::Solve");
 
69
 
 
70
  const int num_rows = A->num_rows();
 
71
  const int num_cols = A->num_cols();
 
72
 
 
73
  if (per_solve_options.D != NULL) {
 
74
    // Temporarily append a diagonal block to the A matrix, but undo
 
75
    // it before returning the matrix to the user.
 
76
    A->AppendDiagonal(per_solve_options.D);
 
77
  }
 
78
 
 
79
  // TODO(sameeragarwal): Since we are copying anyways, the diagonal
 
80
  // can be appended to the matrix instead of doing it on A.
 
81
  lhs_ =  A->matrix();
 
82
 
 
83
  if (per_solve_options.D != NULL) {
 
84
    // Undo the modifications to the matrix A.
 
85
    A->RemoveDiagonal();
 
86
  }
 
87
 
 
88
  // rhs = [b;0] to account for the additional rows in the lhs.
 
89
  if (rhs_.rows() != lhs_.rows()) {
 
90
    rhs_.resize(lhs_.rows());
 
91
  }
 
92
  rhs_.setZero();
 
93
  rhs_.head(num_rows) = ConstVectorRef(b, num_rows);
 
94
 
 
95
  if (work_.rows() == 1) {
 
96
    const int work_size =
 
97
        LAPACK::EstimateWorkSizeForQR(lhs_.rows(), lhs_.cols());
 
98
    VLOG(3) << "Working memory for Dense QR factorization: "
 
99
            << work_size * sizeof(double);
 
100
    work_.resize(work_size);
 
101
  }
 
102
 
 
103
  const int info = LAPACK::SolveUsingQR(lhs_.rows(),
 
104
                                        lhs_.cols(),
 
105
                                        lhs_.data(),
 
106
                                        work_.rows(),
 
107
                                        work_.data(),
 
108
                                        rhs_.data());
 
109
  event_logger.AddEvent("Solve");
 
110
 
 
111
  LinearSolver::Summary summary;
 
112
  summary.num_iterations = 1;
 
113
  if (info == 0) {
 
114
    VectorRef(x, num_cols) = rhs_.head(num_cols);
 
115
    summary.termination_type = TOLERANCE;
 
116
  } else {
 
117
    summary.termination_type = FAILURE;
 
118
  }
 
119
 
 
120
  event_logger.AddEvent("TearDown");
 
121
  return summary;
 
122
}
 
123
 
 
124
LinearSolver::Summary DenseQRSolver::SolveUsingEigen(
 
125
    DenseSparseMatrix* A,
 
126
    const double* b,
 
127
    const LinearSolver::PerSolveOptions& per_solve_options,
 
128
    double* x) {
54
129
  EventLogger event_logger("DenseQRSolver::Solve");
55
130
 
56
131
  const int num_rows = A->num_rows();
73
148
  event_logger.AddEvent("Setup");
74
149
 
75
150
  // Solve the system.
76
 
  VectorRef(x, num_cols) = A->matrix().colPivHouseholderQr().solve(rhs_);
 
151
  VectorRef(x, num_cols) = A->matrix().householderQr().solve(rhs_);
77
152
  event_logger.AddEvent("Solve");
78
153
 
79
154
  if (per_solve_options.D != NULL) {