00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_PART_H
00027 #define EIGEN_PART_H
00028
00046 template<typename MatrixType, unsigned int Mode>
00047 struct ei_traits<Part<MatrixType, Mode> > : ei_traits<MatrixType>
00048 {
00049 typedef typename ei_nested<MatrixType>::type MatrixTypeNested;
00050 typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
00051 enum {
00052 Flags = (_MatrixTypeNested::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode,
00053 CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ConditionalJumpCost
00054 };
00055 };
00056
00057 template<typename MatrixType, unsigned int Mode> class Part
00058 : public MatrixBase<Part<MatrixType, Mode> >
00059 {
00060 public:
00061
00062 EIGEN_GENERIC_PUBLIC_INTERFACE(Part)
00063
00064 inline Part(const MatrixType& matrix) : m_matrix(matrix)
00065 { ei_assert(ei_are_flags_consistent<Mode>::ret); }
00066
00068 template<typename Other> Part& operator+=(const Other& other);
00070 template<typename Other> Part& operator-=(const Other& other);
00072 Part& operator*=(const typename ei_traits<MatrixType>::Scalar& other);
00074 Part& operator/=(const typename ei_traits<MatrixType>::Scalar& other);
00075
00077 template<typename Other> void lazyAssign(const Other& other);
00079 template<typename Other> Part& operator=(const Other& other);
00080
00081 inline int rows() const { return m_matrix.rows(); }
00082 inline int cols() const { return m_matrix.cols(); }
00083 inline int stride() const { return m_matrix.stride(); }
00084
00085 inline Scalar coeff(int row, int col) const
00086 {
00087
00088
00089 if( ((Flags & LowerTriangularBit) && (col>row)) || ((Flags & UpperTriangularBit) && (row>col)) )
00090 return (Scalar)0;
00091 if(Flags & UnitDiagBit)
00092 return col==row ? (Scalar)1 : m_matrix.coeff(row, col);
00093 else if(Flags & ZeroDiagBit)
00094 return col==row ? (Scalar)0 : m_matrix.coeff(row, col);
00095 else
00096 return m_matrix.coeff(row, col);
00097 }
00098
00099 inline Scalar& coeffRef(int row, int col)
00100 {
00101 EIGEN_STATIC_ASSERT(!(Flags & UnitDiagBit), WRITING_TO_TRIANGULAR_PART_WITH_UNIT_DIAGONAL_IS_NOT_SUPPORTED)
00102 EIGEN_STATIC_ASSERT(!(Flags & SelfAdjointBit), COEFFICIENT_WRITE_ACCESS_TO_SELFADJOINT_NOT_SUPPORTED)
00103 ei_assert( (Mode==UpperTriangular && col>=row)
00104 || (Mode==LowerTriangular && col<=row)
00105 || (Mode==StrictlyUpperTriangular && col>row)
00106 || (Mode==StrictlyLowerTriangular && col<row));
00107 return m_matrix.const_cast_derived().coeffRef(row, col);
00108 }
00109
00111 const MatrixType& _expression() const { return m_matrix; }
00112
00114 const Block<Part, 1, ColsAtCompileTime> row(int i) { return Base::row(i); }
00115 const Block<Part, 1, ColsAtCompileTime> row(int i) const { return Base::row(i); }
00117 const Block<Part, RowsAtCompileTime, 1> col(int i) { return Base::col(i); }
00118 const Block<Part, RowsAtCompileTime, 1> col(int i) const { return Base::col(i); }
00119
00120 template<typename OtherDerived>
00121 void swap(const MatrixBase<OtherDerived>& other)
00122 {
00123 Part<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived());
00124 }
00125
00126 protected:
00127 const typename MatrixType::Nested m_matrix;
00128
00129 private:
00130 Part& operator=(const Part&);
00131 };
00132
00146 template<typename Derived>
00147 template<unsigned int Mode>
00148 const Part<Derived, Mode> MatrixBase<Derived>::part() const
00149 {
00150 return derived();
00151 }
00152
00153 template<typename MatrixType, unsigned int Mode>
00154 template<typename Other>
00155 inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator=(const Other& other)
00156 {
00157 if(Other::Flags & EvalBeforeAssigningBit)
00158 {
00159 typename MatrixBase<Other>::PlainMatrixType other_evaluated(other.rows(), other.cols());
00160 other_evaluated.template part<Mode>().lazyAssign(other);
00161 lazyAssign(other_evaluated);
00162 }
00163 else
00164 lazyAssign(other.derived());
00165 return *this;
00166 }
00167
00168 template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount>
00169 struct ei_part_assignment_impl
00170 {
00171 enum {
00172 col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
00173 row = (UnrollCount-1) % Derived1::RowsAtCompileTime
00174 };
00175
00176 inline static void run(Derived1 &dst, const Derived2 &src)
00177 {
00178 ei_part_assignment_impl<Derived1, Derived2, Mode, UnrollCount-1>::run(dst, src);
00179
00180 if(Mode == SelfAdjoint)
00181 {
00182 if(row == col)
00183 dst.coeffRef(row, col) = ei_real(src.coeff(row, col));
00184 else if(row < col)
00185 dst.coeffRef(col, row) = ei_conj(dst.coeffRef(row, col) = src.coeff(row, col));
00186 }
00187 else
00188 {
00189 ei_assert(Mode == UpperTriangular || Mode == LowerTriangular || Mode == StrictlyUpperTriangular || Mode == StrictlyLowerTriangular);
00190 if((Mode == UpperTriangular && row <= col)
00191 || (Mode == LowerTriangular && row >= col)
00192 || (Mode == StrictlyUpperTriangular && row < col)
00193 || (Mode == StrictlyLowerTriangular && row > col))
00194 dst.copyCoeff(row, col, src);
00195 }
00196 }
00197 };
00198
00199 template<typename Derived1, typename Derived2, unsigned int Mode>
00200 struct ei_part_assignment_impl<Derived1, Derived2, Mode, 1>
00201 {
00202 inline static void run(Derived1 &dst, const Derived2 &src)
00203 {
00204 if(!(Mode & ZeroDiagBit))
00205 dst.copyCoeff(0, 0, src);
00206 }
00207 };
00208
00209
00210 template<typename Derived1, typename Derived2, unsigned int Mode>
00211 struct ei_part_assignment_impl<Derived1, Derived2, Mode, 0>
00212 {
00213 inline static void run(Derived1 &, const Derived2 &) {}
00214 };
00215
00216 template<typename Derived1, typename Derived2>
00217 struct ei_part_assignment_impl<Derived1, Derived2, UpperTriangular, Dynamic>
00218 {
00219 inline static void run(Derived1 &dst, const Derived2 &src)
00220 {
00221 for(int j = 0; j < dst.cols(); ++j)
00222 for(int i = 0; i <= j; ++i)
00223 dst.copyCoeff(i, j, src);
00224 }
00225 };
00226
00227 template<typename Derived1, typename Derived2>
00228 struct ei_part_assignment_impl<Derived1, Derived2, LowerTriangular, Dynamic>
00229 {
00230 inline static void run(Derived1 &dst, const Derived2 &src)
00231 {
00232 for(int j = 0; j < dst.cols(); ++j)
00233 for(int i = j; i < dst.rows(); ++i)
00234 dst.copyCoeff(i, j, src);
00235 }
00236 };
00237
00238 template<typename Derived1, typename Derived2>
00239 struct ei_part_assignment_impl<Derived1, Derived2, StrictlyUpperTriangular, Dynamic>
00240 {
00241 inline static void run(Derived1 &dst, const Derived2 &src)
00242 {
00243 for(int j = 0; j < dst.cols(); ++j)
00244 for(int i = 0; i < j; ++i)
00245 dst.copyCoeff(i, j, src);
00246 }
00247 };
00248 template<typename Derived1, typename Derived2>
00249 struct ei_part_assignment_impl<Derived1, Derived2, StrictlyLowerTriangular, Dynamic>
00250 {
00251 inline static void run(Derived1 &dst, const Derived2 &src)
00252 {
00253 for(int j = 0; j < dst.cols(); ++j)
00254 for(int i = j+1; i < dst.rows(); ++i)
00255 dst.copyCoeff(i, j, src);
00256 }
00257 };
00258 template<typename Derived1, typename Derived2>
00259 struct ei_part_assignment_impl<Derived1, Derived2, SelfAdjoint, Dynamic>
00260 {
00261 inline static void run(Derived1 &dst, const Derived2 &src)
00262 {
00263 for(int j = 0; j < dst.cols(); ++j)
00264 {
00265 for(int i = 0; i < j; ++i)
00266 dst.coeffRef(j, i) = ei_conj(dst.coeffRef(i, j) = src.coeff(i, j));
00267 dst.coeffRef(j, j) = ei_real(src.coeff(j, j));
00268 }
00269 }
00270 };
00271
00272 template<typename MatrixType, unsigned int Mode>
00273 template<typename Other>
00274 void Part<MatrixType, Mode>::lazyAssign(const Other& other)
00275 {
00276 const bool unroll = MatrixType::SizeAtCompileTime * Other::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT;
00277 ei_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
00278
00279 ei_part_assignment_impl
00280 <MatrixType, Other, Mode,
00281 unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic
00282 >::run(m_matrix.const_cast_derived(), other.derived());
00283 }
00284
00298 template<typename Derived>
00299 template<unsigned int Mode>
00300 inline Part<Derived, Mode> MatrixBase<Derived>::part()
00301 {
00302 return Part<Derived, Mode>(derived());
00303 }
00304
00310 template<typename Derived>
00311 bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const
00312 {
00313 if(cols() != rows()) return false;
00314 RealScalar maxAbsOnUpperTriangularPart = static_cast<RealScalar>(-1);
00315 for(int j = 0; j < cols(); ++j)
00316 for(int i = 0; i <= j; ++i)
00317 {
00318 RealScalar absValue = ei_abs(coeff(i,j));
00319 if(absValue > maxAbsOnUpperTriangularPart) maxAbsOnUpperTriangularPart = absValue;
00320 }
00321 for(int j = 0; j < cols()-1; ++j)
00322 for(int i = j+1; i < rows(); ++i)
00323 if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperTriangularPart, prec)) return false;
00324 return true;
00325 }
00326
00332 template<typename Derived>
00333 bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const
00334 {
00335 if(cols() != rows()) return false;
00336 RealScalar maxAbsOnLowerTriangularPart = static_cast<RealScalar>(-1);
00337 for(int j = 0; j < cols(); ++j)
00338 for(int i = j; i < rows(); ++i)
00339 {
00340 RealScalar absValue = ei_abs(coeff(i,j));
00341 if(absValue > maxAbsOnLowerTriangularPart) maxAbsOnLowerTriangularPart = absValue;
00342 }
00343 for(int j = 1; j < cols(); ++j)
00344 for(int i = 0; i < j; ++i)
00345 if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerTriangularPart, prec)) return false;
00346 return true;
00347 }
00348
00349 template<typename MatrixType, unsigned int Mode>
00350 template<typename Other>
00351 inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator+=(const Other& other)
00352 {
00353 return *this = m_matrix + other;
00354 }
00355
00356 template<typename MatrixType, unsigned int Mode>
00357 template<typename Other>
00358 inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator-=(const Other& other)
00359 {
00360 return *this = m_matrix - other;
00361 }
00362
00363 template<typename MatrixType, unsigned int Mode>
00364 inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator*=
00365 (const typename ei_traits<MatrixType>::Scalar& other)
00366 {
00367 return *this = m_matrix * other;
00368 }
00369
00370 template<typename MatrixType, unsigned int Mode>
00371 inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator/=
00372 (const typename ei_traits<MatrixType>::Scalar& other)
00373 {
00374 return *this = m_matrix / other;
00375 }
00376
00377 #endif // EIGEN_PART_H