← Back to branch summary

~ubuntu-branches/ubuntu/trusty/blender/trusty-proposed

« back to all changes in this revision

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

  • Committer: Package Import Robot
  • Author(s): Matteo F. Vescovi
  • Date: 2013-08-14 10:43:49 UTC
  • mfrom: (14.2.19 sid)
  • Revision ID: package-import@ubuntu.com-20130814104349-t1d5mtwkphp12dyj
Tags: 2.68a-3
* Upload to unstable
* debian/: python3.3 Depends simplified
  - debian/control: python3.3 Depends dropped
    for blender-data package
  - 0001-blender_thumbnailer.patch refreshed
* debian/control: libavcodec b-dep versioning dropped

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/libmv/third_party/ceres/internal/ceres/schur_complement_solver.cc
38
38
#endif  // CERES_NO_CXSPARSE
39
39
 
40
40
#include "Eigen/Dense"
41
 
#include "glog/logging.h"
42
41
#include "ceres/block_random_access_dense_matrix.h"
43
42
#include "ceres/block_random_access_matrix.h"
44
43
#include "ceres/block_random_access_sparse_matrix.h"
45
44
#include "ceres/block_sparse_matrix.h"
46
45
#include "ceres/block_structure.h"
47
46
#include "ceres/detect_structure.h"
 
47
#include "ceres/internal/eigen.h"
 
48
#include "ceres/internal/port.h"
 
49
#include "ceres/internal/scoped_ptr.h"
48
50
#include "ceres/linear_solver.h"
49
51
#include "ceres/schur_complement_solver.h"
50
52
#include "ceres/suitesparse.h"
51
53
#include "ceres/triplet_sparse_matrix.h"
52
 
#include "ceres/internal/eigen.h"
53
 
#include "ceres/internal/port.h"
54
 
#include "ceres/internal/scoped_ptr.h"
55
54
#include "ceres/types.h"
56
 
 
 
55
#include "ceres/wall_time.h"
57
56
 
58
57
namespace ceres {
59
58
namespace internal {
63
62
    const double* b,
64
63
    const LinearSolver::PerSolveOptions& per_solve_options,
65
64
    double* x) {
66
 
  const time_t start_time = time(NULL);
 
65
  EventLogger event_logger("SchurComplementSolver::Solve");
 
66
 
67
67
  if (eliminator_.get() == NULL) {
68
68
    InitStorage(A->block_structure());
69
69
    DetectStructure(*A->block_structure(),
70
 
                    options_.num_eliminate_blocks,
 
70
                    options_.elimination_groups[0],
71
71
                    &options_.row_block_size,
72
72
                    &options_.e_block_size,
73
73
                    &options_.f_block_size);
74
74
    eliminator_.reset(CHECK_NOTNULL(SchurEliminatorBase::Create(options_)));
75
 
    eliminator_->Init(options_.num_eliminate_blocks, A->block_structure());
 
75
    eliminator_->Init(options_.elimination_groups[0], A->block_structure());
76
76
  };
77
 
  const time_t init_time = time(NULL);
78
77
  fill(x, x + A->num_cols(), 0.0);
 
78
  event_logger.AddEvent("Setup");
79
79
 
80
80
  LinearSolver::Summary summary;
81
81
  summary.num_iterations = 1;
82
82
  summary.termination_type = FAILURE;
83
83
  eliminator_->Eliminate(A, b, per_solve_options.D, lhs_.get(), rhs_.get());
84
 
  const time_t eliminate_time = time(NULL);
 
84
  event_logger.AddEvent("Eliminate");
85
85
 
86
86
  double* reduced_solution = x + A->num_cols() - lhs_->num_cols();
87
87
  const bool status = SolveReducedLinearSystem(reduced_solution);
88
 
  const time_t solve_time = time(NULL);
 
88
  event_logger.AddEvent("ReducedSolve");
89
89
 
90
90
  if (!status) {
91
91
    return summary;
92
92
  }
93
93
 
94
94
  eliminator_->BackSubstitute(A, b, per_solve_options.D, reduced_solution, x);
95
 
  const time_t backsubstitute_time = time(NULL);
96
95
  summary.termination_type = TOLERANCE;
97
96
 
98
 
  VLOG(2) << "time (sec) total: " << (backsubstitute_time - start_time)
99
 
          << " init: " << (init_time - start_time)
100
 
          << " eliminate: " << (eliminate_time - init_time)
101
 
          << " solve: " << (solve_time - eliminate_time)
102
 
          << " backsubstitute: " << (backsubstitute_time - solve_time);
 
97
  event_logger.AddEvent("BackSubstitute");
103
98
  return summary;
104
99
}
105
100
 
107
102
// complement.
108
103
void DenseSchurComplementSolver::InitStorage(
109
104
    const CompressedRowBlockStructure* bs) {
110
 
  const int num_eliminate_blocks = options().num_eliminate_blocks;
 
105
  const int num_eliminate_blocks = options().elimination_groups[0];
111
106
  const int num_col_blocks = bs->cols.size();
112
107
 
113
108
  vector<int> blocks(num_col_blocks - num_eliminate_blocks, 0);
146
141
  return true;
147
142
}
148
143
 
 
144
#if !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARE)
149
145
 
150
146
SparseSchurComplementSolver::SparseSchurComplementSolver(
151
147
    const LinearSolver::Options& options)
179
175
// initialize a BlockRandomAccessSparseMatrix object.
180
176
void SparseSchurComplementSolver::InitStorage(
181
177
    const CompressedRowBlockStructure* bs) {
182
 
  const int num_eliminate_blocks = options().num_eliminate_blocks;
 
178
  const int num_eliminate_blocks = options().elimination_groups[0];
183
179
  const int num_col_blocks = bs->cols.size();
184
180
  const int num_row_blocks = bs->rows.size();
185
181
 
268
264
// CHOLMOD's sparse cholesky factorization routines.
269
265
bool SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
270
266
    double* solution) {
271
 
  const time_t start_time = time(NULL);
272
 
 
273
267
  TripletSparseMatrix* tsm =
274
268
      const_cast<TripletSparseMatrix*>(
275
269
          down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix());
286
280
  // The matrix is symmetric, and the upper triangular part of the
287
281
  // matrix contains the values.
288
282
  cholmod_lhs->stype = 1;
289
 
  const time_t lhs_time = time(NULL);
290
283
 
291
284
  cholmod_dense*  cholmod_rhs =
292
285
      ss_.CreateDenseVector(const_cast<double*>(rhs()), num_rows, num_rows);
293
 
  const time_t rhs_time = time(NULL);
294
286
 
295
287
  // Symbolic factorization is computed if we don't already have one handy.
296
288
  if (factor_ == NULL) {
297
 
    if (options().use_block_amd) {
298
 
      factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs, blocks_, blocks_);
299
 
    } else {
300
 
      factor_ = ss_.AnalyzeCholesky(cholmod_lhs);
301
 
    }
302
 
 
303
 
    if (VLOG_IS_ON(2)) {
304
 
      cholmod_print_common("Symbolic Analysis", ss_.mutable_cc());
305
 
    }
 
289
    factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs, blocks_, blocks_);
306
290
  }
307
291
 
308
 
  CHECK_NOTNULL(factor_);
309
 
 
310
 
  const time_t symbolic_time = time(NULL);
311
292
  cholmod_dense* cholmod_solution =
312
293
      ss_.SolveCholesky(cholmod_lhs, factor_, cholmod_rhs);
313
294
 
314
 
  const time_t solve_time = time(NULL);
315
 
 
316
295
  ss_.Free(cholmod_lhs);
317
296
  cholmod_lhs = NULL;
318
297
  ss_.Free(cholmod_rhs);
319
298
  cholmod_rhs = NULL;
320
299
 
321
300
  if (cholmod_solution == NULL) {
322
 
    LOG(ERROR) << "CHOLMOD solve failed.";
 
301
    LOG(WARNING) << "CHOLMOD solve failed.";
323
302
    return false;
324
303
  }
325
304
 
326
305
  VectorRef(solution, num_rows)
327
306
      = VectorRef(static_cast<double*>(cholmod_solution->x), num_rows);
328
307
  ss_.Free(cholmod_solution);
329
 
  const time_t final_time = time(NULL);
330
 
  VLOG(2) << "time: " << (final_time - start_time)
331
 
          << " lhs : " << (lhs_time - start_time)
332
 
          << " rhs:  " << (rhs_time - lhs_time)
333
 
          << " analyze: " <<  (symbolic_time - rhs_time)
334
 
          << " factor_and_solve: " << (solve_time - symbolic_time)
335
 
          << " cleanup: " << (final_time - solve_time);
336
308
  return true;
337
309
}
338
310
#else
384
356
}
385
357
#endif  // CERES_NO_CXPARSE
386
358
 
 
359
#endif  // !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARE)
387
360
}  // namespace internal
388
361
}  // namespace ceres