src/Core/SolverBase.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SOLVERBASE_H
0011 #define EIGEN_SOLVERBASE_H
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016
0017 template<typename Derived>
0018 struct solve_assertion {
0019     template<bool Transpose_, typename Rhs>
0020     static void run(const Derived& solver, const Rhs& b) { solver.template _check_solve_assertion<Transpose_>(b); }
0021 };
0022
0023 template<typename Derived>
0024 struct solve_assertion<Transpose<Derived> >
0025 {
0026     typedef Transpose<Derived> type;
0027
0028     template<bool Transpose_, typename Rhs>
0029     static void run(const type& transpose, const Rhs& b)
0030     {
0031         internal::solve_assertion<typename internal::remove_all<Derived>::type>::template run<true>(transpose.nestedExpression(), b);
0032     }
0033 };
0034
0035 template<typename Scalar, typename Derived>
0036 struct solve_assertion<CwiseUnaryOp<Eigen::internal::scalar_conjugate_op<Scalar>, const Transpose<Derived> > >
0037 {
0038     typedef CwiseUnaryOp<Eigen::internal::scalar_conjugate_op<Scalar>, const Transpose<Derived> > type;
0039
0040     template<bool Transpose_, typename Rhs>
0041     static void run(const type& adjoint, const Rhs& b)
0042     {
0043         internal::solve_assertion<typename internal::remove_all<Transpose<Derived> >::type>::template run<true>(adjoint.nestedExpression(), b);
0044     }
0045 };
0046 } // end namespace internal
0047
0048 /** \class SolverBase
0049   * \brief A base class for matrix decomposition and solvers
0050   *
0051   * \tparam Derived the actual type of the decomposition/solver.
0052   *
0053   * Any matrix decomposition inheriting this base class provide the following API:
0054   *
0055   * \code
0056   * MatrixType A, b, x;
0057   * DecompositionType dec(A);
0058   * x = dec.solve(b);             // solve A   * x = b
0059   * x = dec.transpose().solve(b); // solve A^T * x = b
0060   * x = dec.adjoint().solve(b);   // solve A'  * x = b
0061   * \endcode
0062   *
0063   * \warning Currently, any other usage of transpose() and adjoint() are not supported and will produce compilation errors.
0064   *
0065   * \sa class PartialPivLU, class FullPivLU, class HouseholderQR, class ColPivHouseholderQR, class FullPivHouseholderQR, class CompleteOrthogonalDecomposition, class LLT, class LDLT, class SVDBase
0066   */
0067 template<typename Derived>
0068 class SolverBase : public EigenBase<Derived>
0069 {
0070   public:
0071
0072     typedef EigenBase<Derived> Base;
0073     typedef typename internal::traits<Derived>::Scalar Scalar;
0074     typedef Scalar CoeffReturnType;
0075
0076     template<typename Derived_>
0077     friend struct internal::solve_assertion;
0078
0079     enum {
0080       RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
0081       ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
0082       SizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::RowsAtCompileTime,
0083                                                           internal::traits<Derived>::ColsAtCompileTime>::ret),
0084       MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
0085       MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime,
0086       MaxSizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::MaxRowsAtCompileTime,
0087                                                              internal::traits<Derived>::MaxColsAtCompileTime>::ret),
0088       IsVectorAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime == 1
0089                            || internal::traits<Derived>::MaxColsAtCompileTime == 1,
0090       NumDimensions = int(MaxSizeAtCompileTime) == 1 ? 0 : bool(IsVectorAtCompileTime) ? 1 : 2
0091     };
0092
0093     /** Default constructor */
0094     SolverBase()
0095     {}
0096
0097     ~SolverBase()
0098     {}
0099
0100     using Base::derived;
0101
0102     /** \returns an expression of the solution x of \f$ A x = b \f$ using the current decomposition of A.
0103       */
0104     template<typename Rhs>
0105     inline const Solve<Derived, Rhs>
0106     solve(const MatrixBase<Rhs>& b) const
0107     {
0108       internal::solve_assertion<typename internal::remove_all<Derived>::type>::template run<false>(derived(), b);
0109       return Solve<Derived, Rhs>(derived(), b.derived());
0110     }
0111
0112     /** \internal the return type of transpose() */
0113     typedef typename internal::add_const<Transpose<const Derived> >::type ConstTransposeReturnType;
0114     /** \returns an expression of the transposed of the factored matrix.
0115       *
0116       * A typical usage is to solve for the transposed problem A^T x = b:
0117       * \code x = dec.transpose().solve(b); \endcode
0118       *
0119       * \sa adjoint(), solve()
0120       */
0121     inline ConstTransposeReturnType transpose() const
0122     {
0123       return ConstTransposeReturnType(derived());
0124     }
0125
0126     /** \internal the return type of adjoint() */
0127     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
0128                         CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, ConstTransposeReturnType>,
0129                         ConstTransposeReturnType
0130                      >::type AdjointReturnType;
0131     /** \returns an expression of the adjoint of the factored matrix
0132       *
0133       * A typical usage is to solve for the adjoint problem A' x = b:
0134       * \code x = dec.adjoint().solve(b); \endcode
0135       *
0136       * For real scalar types, this function is equivalent to transpose().
0137       *
0138       * \sa transpose(), solve()
0139       */
0140     inline AdjointReturnType adjoint() const
0141     {
0142       return AdjointReturnType(derived().transpose());
0143     }
0144
0145   protected:
0146
0147     template<bool Transpose_, typename Rhs>
0148     void _check_solve_assertion(const Rhs& b) const {
0149         EIGEN_ONLY_USED_FOR_DEBUG(b);
0150         eigen_assert(derived().m_isInitialized && "Solver is not initialized.");
0151         eigen_assert((Transpose_?derived().cols():derived().rows())==b.rows() && "SolverBase::solve(): invalid number of rows of the right hand side matrix b");
0152     }
0153 };
0154
0155 namespace internal {
0156
0157 template<typename Derived>
0158 struct generic_xpr_base<Derived, MatrixXpr, SolverStorage>
0159 {
0160   typedef SolverBase<Derived> type;
0161
0162 };
0163
0164 } // end namespace internal
0165
0166 } // end namespace Eigen
0167
0168 #endif // EIGEN_SOLVERBASE_H