dense_qr.h 8.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199
  1. // Ceres Solver - A fast non-linear least squares minimizer
  2. // Copyright 2023 Google Inc. All rights reserved.
  3. // http://ceres-solver.org/
  4. //
  5. // Redistribution and use in source and binary forms, with or without
  6. // modification, are permitted provided that the following conditions are met:
  7. //
  8. // * Redistributions of source code must retain the above copyright notice,
  9. // this list of conditions and the following disclaimer.
  10. // * Redistributions in binary form must reproduce the above copyright notice,
  11. // this list of conditions and the following disclaimer in the documentation
  12. // and/or other materials provided with the distribution.
  13. // * Neither the name of Google Inc. nor the names of its contributors may be
  14. // used to endorse or promote products derived from this software without
  15. // specific prior written permission.
  16. //
  17. // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18. // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19. // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20. // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  21. // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22. // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23. // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24. // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25. // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26. // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27. // POSSIBILITY OF SUCH DAMAGE.
  28. //
  29. // Author: sameeragarwal@google.com (Sameer Agarwal)
  30. #ifndef CERES_INTERNAL_DENSE_QR_H_
  31. #define CERES_INTERNAL_DENSE_QR_H_
  32. // This include must come before any #ifndef check on Ceres compile options.
  33. // clang-format off
  34. #include "ceres/internal/config.h"
  35. // clang-format on
  36. #include <memory>
  37. #include <vector>
  38. #include "Eigen/Dense"
  39. #include "ceres/context_impl.h"
  40. #include "ceres/internal/disable_warnings.h"
  41. #include "ceres/internal/eigen.h"
  42. #include "ceres/internal/export.h"
  43. #include "ceres/linear_solver.h"
  44. #include "glog/logging.h"
  45. #ifndef CERES_NO_CUDA
  46. #include "ceres/context_impl.h"
  47. #include "ceres/cuda_buffer.h"
  48. #include "cublas_v2.h"
  49. #include "cuda_runtime.h"
  50. #include "cusolverDn.h"
  51. #endif // CERES_NO_CUDA
  52. namespace ceres::internal {
  53. // An interface that abstracts away the internal details of various dense linear
  54. // algebra libraries and offers a simple API for solving dense linear systems
  55. // using a QR factorization.
  56. class CERES_NO_EXPORT DenseQR {
  57. public:
  58. static std::unique_ptr<DenseQR> Create(const LinearSolver::Options& options);
  59. virtual ~DenseQR();
  60. // Computes the QR factorization of the given matrix.
  61. //
  62. // The input matrix lhs is assumed to be a column-major num_rows x num_cols
  63. // matrix.
  64. //
  65. // The input matrix lhs may be modified by the implementation to store the
  66. // factorization, irrespective of whether the factorization succeeds or not.
  67. // As a result it is the user's responsibility to ensure that lhs is valid
  68. // when Solve is called.
  69. virtual LinearSolverTerminationType Factorize(int num_rows,
  70. int num_cols,
  71. double* lhs,
  72. std::string* message) = 0;
  73. // Computes the solution to the equation
  74. //
  75. // lhs * solution = rhs
  76. //
  77. // Calling Solve without calling Factorize is undefined behaviour. It is the
  78. // user's responsibility to ensure that the input matrix lhs passed to
  79. // Factorize has not been freed/modified when Solve is called.
  80. virtual LinearSolverTerminationType Solve(const double* rhs,
  81. double* solution,
  82. std::string* message) = 0;
  83. // Convenience method which combines a call to Factorize and Solve. Solve is
  84. // only called if Factorize returns LinearSolverTerminationType::SUCCESS.
  85. //
  86. // The input matrix lhs may be modified by the implementation to store the
  87. // factorization, irrespective of whether the method succeeds or not. It is
  88. // the user's responsibility to ensure that lhs is valid if and when Solve is
  89. // called again after this call.
  90. LinearSolverTerminationType FactorAndSolve(int num_rows,
  91. int num_cols,
  92. double* lhs,
  93. const double* rhs,
  94. double* solution,
  95. std::string* message);
  96. };
  97. class CERES_NO_EXPORT EigenDenseQR final : public DenseQR {
  98. public:
  99. LinearSolverTerminationType Factorize(int num_rows,
  100. int num_cols,
  101. double* lhs,
  102. std::string* message) override;
  103. LinearSolverTerminationType Solve(const double* rhs,
  104. double* solution,
  105. std::string* message) override;
  106. private:
  107. using QRType = Eigen::HouseholderQR<Eigen::Ref<ColMajorMatrix>>;
  108. std::unique_ptr<QRType> qr_;
  109. };
  110. #ifndef CERES_NO_LAPACK
  111. class CERES_NO_EXPORT LAPACKDenseQR final : public DenseQR {
  112. public:
  113. LinearSolverTerminationType Factorize(int num_rows,
  114. int num_cols,
  115. double* lhs,
  116. std::string* message) override;
  117. LinearSolverTerminationType Solve(const double* rhs,
  118. double* solution,
  119. std::string* message) override;
  120. private:
  121. double* lhs_ = nullptr;
  122. int num_rows_;
  123. int num_cols_;
  124. LinearSolverTerminationType termination_type_ =
  125. LinearSolverTerminationType::FATAL_ERROR;
  126. Vector work_;
  127. Vector tau_;
  128. Vector q_transpose_rhs_;
  129. };
  130. #endif // CERES_NO_LAPACK
  131. #ifndef CERES_NO_CUDA
  132. // Implementation of DenseQR using the 32-bit cuSolverDn interface. A
  133. // requirement for using this solver is that the lhs must not be rank deficient.
  134. // This is because cuSolverDn does not implement the singularity-checking
  135. // wrapper trtrs, hence this solver directly uses trsv from CUBLAS for the
  136. // backsubstitution.
  137. class CERES_NO_EXPORT CUDADenseQR final : public DenseQR {
  138. public:
  139. static std::unique_ptr<CUDADenseQR> Create(
  140. const LinearSolver::Options& options);
  141. CUDADenseQR(const CUDADenseQR&) = delete;
  142. CUDADenseQR& operator=(const CUDADenseQR&) = delete;
  143. LinearSolverTerminationType Factorize(int num_rows,
  144. int num_cols,
  145. double* lhs,
  146. std::string* message) override;
  147. LinearSolverTerminationType Solve(const double* rhs,
  148. double* solution,
  149. std::string* message) override;
  150. private:
  151. explicit CUDADenseQR(ContextImpl* context);
  152. ContextImpl* context_ = nullptr;
  153. // Number of rowns in the A matrix, to be cached between calls to *Factorize
  154. // and *Solve.
  155. size_t num_rows_ = 0;
  156. // Number of columns in the A matrix, to be cached between calls to *Factorize
  157. // and *Solve.
  158. size_t num_cols_ = 0;
  159. // GPU memory allocated for the A matrix (lhs matrix).
  160. CudaBuffer<double> lhs_;
  161. // GPU memory allocated for the B matrix (rhs vector).
  162. CudaBuffer<double> rhs_;
  163. // GPU memory allocated for the TAU matrix (scaling of householder vectors).
  164. CudaBuffer<double> tau_;
  165. // Scratch space for cuSOLVER on the GPU.
  166. CudaBuffer<double> device_workspace_;
  167. // Required for error handling with cuSOLVER.
  168. CudaBuffer<int> error_;
  169. // Cache the result of Factorize to ensure that when Solve is called, the
  170. // factiorization of lhs is valid.
  171. LinearSolverTerminationType factorize_result_ =
  172. LinearSolverTerminationType::FATAL_ERROR;
  173. };
  174. #endif // CERES_NO_CUDA
  175. } // namespace ceres::internal
  176. #include "ceres/internal/reenable_warnings.h"
  177. #endif // CERES_INTERNAL_DENSE_QR_H_