* replace postfix ++ by prefix ++ wherever that makes sense in Eigen/

* fix some "unused variable" warnings in the tests; there remains a libstdc++ "deprecated"
warning which I haven't looked much into
This commit is contained in:
Benoit Jacob 2008-12-17 14:30:01 +00:00
parent 2110cca431
commit 89f468671d
28 changed files with 164 additions and 162 deletions

View File

@ -99,8 +99,8 @@ inline bool MatrixBase<Derived>::all(void) const
>::run(derived()); >::run(derived());
else else
{ {
for(int j = 0; j < cols(); j++) for(int j = 0; j < cols(); ++j)
for(int i = 0; i < rows(); i++) for(int i = 0; i < rows(); ++i)
if (!coeff(i, j)) return false; if (!coeff(i, j)) return false;
return true; return true;
} }
@ -123,8 +123,8 @@ inline bool MatrixBase<Derived>::any(void) const
>::run(derived()); >::run(derived());
else else
{ {
for(int j = 0; j < cols(); j++) for(int j = 0; j < cols(); ++j)
for(int i = 0; i < rows(); i++) for(int i = 0; i < rows(); ++i)
if (coeff(i, j)) return true; if (coeff(i, j)) return true;
return false; return false;
} }

View File

@ -214,8 +214,8 @@ struct ei_assign_impl<Derived1, Derived2, NoVectorization, NoUnrolling>
{ {
const int innerSize = dst.innerSize(); const int innerSize = dst.innerSize();
const int outerSize = dst.outerSize(); const int outerSize = dst.outerSize();
for(int j = 0; j < outerSize; j++) for(int j = 0; j < outerSize; ++j)
for(int i = 0; i < innerSize; i++) for(int i = 0; i < innerSize; ++i)
{ {
if(int(Derived1::Flags)&RowMajorBit) if(int(Derived1::Flags)&RowMajorBit)
dst.copyCoeff(j, i, src); dst.copyCoeff(j, i, src);
@ -243,7 +243,7 @@ struct ei_assign_impl<Derived1, Derived2, NoVectorization, InnerUnrolling>
const bool rowMajor = int(Derived1::Flags)&RowMajorBit; const bool rowMajor = int(Derived1::Flags)&RowMajorBit;
const int innerSize = rowMajor ? Derived1::ColsAtCompileTime : Derived1::RowsAtCompileTime; const int innerSize = rowMajor ? Derived1::ColsAtCompileTime : Derived1::RowsAtCompileTime;
const int outerSize = dst.outerSize(); const int outerSize = dst.outerSize();
for(int j = 0; j < outerSize; j++) for(int j = 0; j < outerSize; ++j)
ei_assign_novec_InnerUnrolling<Derived1, Derived2, 0, innerSize> ei_assign_novec_InnerUnrolling<Derived1, Derived2, 0, innerSize>
::run(dst, src, j); ::run(dst, src, j);
} }
@ -261,7 +261,7 @@ struct ei_assign_impl<Derived1, Derived2, InnerVectorization, NoUnrolling>
const int innerSize = dst.innerSize(); const int innerSize = dst.innerSize();
const int outerSize = dst.outerSize(); const int outerSize = dst.outerSize();
const int packetSize = ei_packet_traits<typename Derived1::Scalar>::size; const int packetSize = ei_packet_traits<typename Derived1::Scalar>::size;
for(int j = 0; j < outerSize; j++) for(int j = 0; j < outerSize; ++j)
for(int i = 0; i < innerSize; i+=packetSize) for(int i = 0; i < innerSize; i+=packetSize)
{ {
if(int(Derived1::Flags)&RowMajorBit) if(int(Derived1::Flags)&RowMajorBit)
@ -290,7 +290,7 @@ struct ei_assign_impl<Derived1, Derived2, InnerVectorization, InnerUnrolling>
const bool rowMajor = int(Derived1::Flags)&RowMajorBit; const bool rowMajor = int(Derived1::Flags)&RowMajorBit;
const int innerSize = rowMajor ? Derived1::ColsAtCompileTime : Derived1::RowsAtCompileTime; const int innerSize = rowMajor ? Derived1::ColsAtCompileTime : Derived1::RowsAtCompileTime;
const int outerSize = dst.outerSize(); const int outerSize = dst.outerSize();
for(int j = 0; j < outerSize; j++) for(int j = 0; j < outerSize; ++j)
ei_assign_innervec_InnerUnrolling<Derived1, Derived2, 0, innerSize> ei_assign_innervec_InnerUnrolling<Derived1, Derived2, 0, innerSize>
::run(dst, src, j); ::run(dst, src, j);
} }
@ -311,7 +311,7 @@ struct ei_assign_impl<Derived1, Derived2, LinearVectorization, NoUnrolling>
: ei_alignmentOffset(&dst.coeffRef(0), size); : ei_alignmentOffset(&dst.coeffRef(0), size);
const int alignedEnd = alignedStart + ((size-alignedStart)/packetSize)*packetSize; const int alignedEnd = alignedStart + ((size-alignedStart)/packetSize)*packetSize;
for(int index = 0; index < alignedStart; index++) for(int index = 0; index < alignedStart; ++index)
dst.copyCoeff(index, src); dst.copyCoeff(index, src);
for(int index = alignedStart; index < alignedEnd; index += packetSize) for(int index = alignedStart; index < alignedEnd; index += packetSize)
@ -319,7 +319,7 @@ struct ei_assign_impl<Derived1, Derived2, LinearVectorization, NoUnrolling>
dst.template copyPacket<Derived2, Aligned, ei_assign_traits<Derived1,Derived2>::SrcAlignment>(index, src); dst.template copyPacket<Derived2, Aligned, ei_assign_traits<Derived1,Derived2>::SrcAlignment>(index, src);
} }
for(int index = alignedEnd; index < size; index++) for(int index = alignedEnd; index < size; ++index)
dst.copyCoeff(index, src); dst.copyCoeff(index, src);
} }
}; };
@ -355,12 +355,12 @@ struct ei_assign_impl<Derived1, Derived2, SliceVectorization, NoUnrolling>
int alignedStart = ei_assign_traits<Derived1,Derived2>::DstIsAligned ? 0 int alignedStart = ei_assign_traits<Derived1,Derived2>::DstIsAligned ? 0
: ei_alignmentOffset(&dst.coeffRef(0), innerSize); : ei_alignmentOffset(&dst.coeffRef(0), innerSize);
for(int i = 0; i < outerSize; i++) for(int i = 0; i < outerSize; ++i)
{ {
const int alignedEnd = alignedStart + ((innerSize-alignedStart) & ~packetAlignedMask); const int alignedEnd = alignedStart + ((innerSize-alignedStart) & ~packetAlignedMask);
// do the non-vectorizable part of the assignment // do the non-vectorizable part of the assignment
for (int index = 0; index<alignedStart ; index++) for (int index = 0; index<alignedStart ; ++index)
{ {
if(Derived1::Flags&RowMajorBit) if(Derived1::Flags&RowMajorBit)
dst.copyCoeff(i, index, src); dst.copyCoeff(i, index, src);
@ -378,7 +378,7 @@ struct ei_assign_impl<Derived1, Derived2, SliceVectorization, NoUnrolling>
} }
// do the non-vectorizable part of the assignment // do the non-vectorizable part of the assignment
for (int index = alignedEnd; index<innerSize ; index++) for (int index = alignedEnd; index<innerSize ; ++index)
{ {
if(Derived1::Flags&RowMajorBit) if(Derived1::Flags&RowMajorBit)
dst.copyCoeff(i, index, src); dst.copyCoeff(i, index, src);

View File

@ -433,7 +433,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
{ {
/* explicit vectorization */ /* explicit vectorization */
// process initial unaligned coeffs // process initial unaligned coeffs
for (int j=0; j<alignedStart; j++) for (int j=0; j<alignedStart; ++j)
res[j] += ei_pfirst(ptmp0)*lhs0[j] + ei_pfirst(ptmp1)*lhs1[j] + ei_pfirst(ptmp2)*lhs2[j] + ei_pfirst(ptmp3)*lhs3[j]; res[j] += ei_pfirst(ptmp0)*lhs0[j] + ei_pfirst(ptmp1)*lhs1[j] + ei_pfirst(ptmp2)*lhs2[j] + ei_pfirst(ptmp3)*lhs3[j];
if (alignedSize>alignedStart) if (alignedSize>alignedStart)
@ -493,7 +493,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
} // end explicit vectorization } // end explicit vectorization
/* process remaining coeffs (or all if there is no explicit vectorization) */ /* process remaining coeffs (or all if there is no explicit vectorization) */
for (int j=alignedSize; j<size; j++) for (int j=alignedSize; j<size; ++j)
res[j] += ei_pfirst(ptmp0)*lhs0[j] + ei_pfirst(ptmp1)*lhs1[j] + ei_pfirst(ptmp2)*lhs2[j] + ei_pfirst(ptmp3)*lhs3[j]; res[j] += ei_pfirst(ptmp0)*lhs0[j] + ei_pfirst(ptmp1)*lhs1[j] + ei_pfirst(ptmp2)*lhs2[j] + ei_pfirst(ptmp3)*lhs3[j];
} }
@ -502,7 +502,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
int start = columnBound; int start = columnBound;
do do
{ {
for (int i=start; i<end; i++) for (int i=start; i<end; ++i)
{ {
Packet ptmp0 = ei_pset1(rhs[i]); Packet ptmp0 = ei_pset1(rhs[i]);
const Scalar* lhs0 = lhs + i*lhsStride; const Scalar* lhs0 = lhs + i*lhsStride;
@ -511,7 +511,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
{ {
/* explicit vectorization */ /* explicit vectorization */
// process first unaligned result's coeffs // process first unaligned result's coeffs
for (int j=0; j<alignedStart; j++) for (int j=0; j<alignedStart; ++j)
res[j] += ei_pfirst(ptmp0) * lhs0[j]; res[j] += ei_pfirst(ptmp0) * lhs0[j];
// process aligned result's coeffs // process aligned result's coeffs
@ -524,7 +524,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
} }
// process remaining scalars (or all if no explicit vectorization) // process remaining scalars (or all if no explicit vectorization)
for (int j=alignedSize; j<size; j++) for (int j=alignedSize; j<size; ++j)
res[j] += ei_pfirst(ptmp0) * lhs0[j]; res[j] += ei_pfirst(ptmp0) * lhs0[j];
} }
if (skipColumns) if (skipColumns)
@ -624,7 +624,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
// process initial unaligned coeffs // process initial unaligned coeffs
// FIXME this loop get vectorized by the compiler ! // FIXME this loop get vectorized by the compiler !
for (int j=0; j<alignedStart; j++) for (int j=0; j<alignedStart; ++j)
{ {
Scalar b = rhs[j]; Scalar b = rhs[j];
tmp0 += b*lhs0[j]; tmp1 += b*lhs1[j]; tmp2 += b*lhs2[j]; tmp3 += b*lhs3[j]; tmp0 += b*lhs0[j]; tmp1 += b*lhs1[j]; tmp2 += b*lhs2[j]; tmp3 += b*lhs3[j];
@ -695,7 +695,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
// process remaining coeffs (or all if no explicit vectorization) // process remaining coeffs (or all if no explicit vectorization)
// FIXME this loop get vectorized by the compiler ! // FIXME this loop get vectorized by the compiler !
for (int j=alignedSize; j<size; j++) for (int j=alignedSize; j<size; ++j)
{ {
Scalar b = rhs[j]; Scalar b = rhs[j];
tmp0 += b*lhs0[j]; tmp1 += b*lhs1[j]; tmp2 += b*lhs2[j]; tmp3 += b*lhs3[j]; tmp0 += b*lhs0[j]; tmp1 += b*lhs1[j]; tmp2 += b*lhs2[j]; tmp3 += b*lhs3[j];
@ -708,14 +708,14 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
int start = rowBound; int start = rowBound;
do do
{ {
for (int i=start; i<end; i++) for (int i=start; i<end; ++i)
{ {
Scalar tmp0 = Scalar(0); Scalar tmp0 = Scalar(0);
Packet ptmp0 = ei_pset1(tmp0); Packet ptmp0 = ei_pset1(tmp0);
const Scalar* lhs0 = lhs + i*lhsStride; const Scalar* lhs0 = lhs + i*lhsStride;
// process first unaligned result's coeffs // process first unaligned result's coeffs
// FIXME this loop get vectorized by the compiler ! // FIXME this loop get vectorized by the compiler !
for (int j=0; j<alignedStart; j++) for (int j=0; j<alignedStart; ++j)
tmp0 += rhs[j] * lhs0[j]; tmp0 += rhs[j] * lhs0[j];
if (alignedSize>alignedStart) if (alignedSize>alignedStart)
@ -732,7 +732,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
// process remaining scalars // process remaining scalars
// FIXME this loop get vectorized by the compiler ! // FIXME this loop get vectorized by the compiler !
for (int j=alignedSize; j<size; j++) for (int j=alignedSize; j<size; ++j)
tmp0 += rhs[j] * lhs0[j]; tmp0 += rhs[j] * lhs0[j];
res[i] += tmp0; res[i] += tmp0;
} }

View File

@ -67,7 +67,7 @@ struct CommaInitializer
ei_assert(m_col<m_matrix.cols() ei_assert(m_col<m_matrix.cols()
&& "Too many coefficients passed to comma initializer (operator<<)"); && "Too many coefficients passed to comma initializer (operator<<)");
ei_assert(m_currentBlockRows==1); ei_assert(m_currentBlockRows==1);
m_matrix.coeffRef(m_row, m_col++) = s; m_matrix.coeffRef(m_row, ++m_col) = s;
return *this; return *this;
} }

View File

@ -236,8 +236,8 @@ template<typename Derived>
bool MatrixBase<Derived>::isApproxToConstant bool MatrixBase<Derived>::isApproxToConstant
(const Scalar& value, RealScalar prec) const (const Scalar& value, RealScalar prec) const
{ {
for(int j = 0; j < cols(); j++) for(int j = 0; j < cols(); ++j)
for(int i = 0; i < rows(); i++) for(int i = 0; i < rows(); ++i)
if(!ei_isApprox(coeff(i, j), value, prec)) if(!ei_isApprox(coeff(i, j), value, prec))
return false; return false;
return true; return true;
@ -330,8 +330,8 @@ template<typename Derived>
bool MatrixBase<Derived>::isZero bool MatrixBase<Derived>::isZero
(RealScalar prec) const (RealScalar prec) const
{ {
for(int j = 0; j < cols(); j++) for(int j = 0; j < cols(); ++j)
for(int i = 0; i < rows(); i++) for(int i = 0; i < rows(); ++i)
if(!ei_isMuchSmallerThan(coeff(i, j), static_cast<Scalar>(1), prec)) if(!ei_isMuchSmallerThan(coeff(i, j), static_cast<Scalar>(1), prec))
return false; return false;
return true; return true;
@ -499,9 +499,9 @@ template<typename Derived>
bool MatrixBase<Derived>::isIdentity bool MatrixBase<Derived>::isIdentity
(RealScalar prec) const (RealScalar prec) const
{ {
for(int j = 0; j < cols(); j++) for(int j = 0; j < cols(); ++j)
{ {
for(int i = 0; i < rows(); i++) for(int i = 0; i < rows(); ++i)
{ {
if(i == j) if(i == j)
{ {
@ -534,7 +534,7 @@ struct ei_setIdentity_impl<Derived, true>
{ {
m.setZero(); m.setZero();
const int size = std::min(m.rows(), m.cols()); const int size = std::min(m.rows(), m.cols());
for(int i = 0; i < size; i++) m.coeffRef(i,i) = typename Derived::Scalar(1); for(int i = 0; i < size; ++i) m.coeffRef(i,i) = typename Derived::Scalar(1);
return m; return m;
} }
}; };

View File

@ -123,13 +123,13 @@ bool MatrixBase<Derived>::isDiagonal
{ {
if(cols() != rows()) return false; if(cols() != rows()) return false;
RealScalar maxAbsOnDiagonal = static_cast<RealScalar>(-1); RealScalar maxAbsOnDiagonal = static_cast<RealScalar>(-1);
for(int j = 0; j < cols(); j++) for(int j = 0; j < cols(); ++j)
{ {
RealScalar absOnDiagonal = ei_abs(coeff(j,j)); RealScalar absOnDiagonal = ei_abs(coeff(j,j));
if(absOnDiagonal > maxAbsOnDiagonal) maxAbsOnDiagonal = absOnDiagonal; if(absOnDiagonal > maxAbsOnDiagonal) maxAbsOnDiagonal = absOnDiagonal;
} }
for(int j = 0; j < cols(); j++) for(int j = 0; j < cols(); ++j)
for(int i = 0; i < j; i++) for(int i = 0; i < j; ++i)
{ {
if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnDiagonal, prec)) return false; if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnDiagonal, prec)) return false;
if(!ei_isMuchSmallerThan(coeff(j, i), maxAbsOnDiagonal, prec)) return false; if(!ei_isMuchSmallerThan(coeff(j, i), maxAbsOnDiagonal, prec)) return false;

View File

@ -156,7 +156,7 @@ struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling>
ei_assert(v1.size()>0 && "you are using a non initialized vector"); ei_assert(v1.size()>0 && "you are using a non initialized vector");
Scalar res; Scalar res;
res = v1.coeff(0) * ei_conj(v2.coeff(0)); res = v1.coeff(0) * ei_conj(v2.coeff(0));
for(int i = 1; i < v1.size(); i++) for(int i = 1; i < v1.size(); ++i)
res += v1.coeff(i) * ei_conj(v2.coeff(i)); res += v1.coeff(i) * ei_conj(v2.coeff(i));
return res; return res;
} }
@ -211,7 +211,7 @@ struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling>
} }
// do the remainder of the vector // do the remainder of the vector
for(int index = alignedSize; index < size; index++) for(int index = alignedSize; index < size; ++index)
{ {
res += v1.coeff(index) * v2.coeff(index); res += v1.coeff(index) * v2.coeff(index);
} }
@ -370,11 +370,11 @@ template<typename Derived>
bool MatrixBase<Derived>::isUnitary(RealScalar prec) const bool MatrixBase<Derived>::isUnitary(RealScalar prec) const
{ {
typename Derived::Nested nested(derived()); typename Derived::Nested nested(derived());
for(int i = 0; i < cols(); i++) for(int i = 0; i < cols(); ++i)
{ {
if(!ei_isApprox(nested.col(i).squaredNorm(), static_cast<Scalar>(1), prec)) if(!ei_isApprox(nested.col(i).squaredNorm(), static_cast<Scalar>(1), prec))
return false; return false;
for(int j = 0; j < i; j++) for(int j = 0; j < i; ++j)
if(!ei_isMuchSmallerThan(nested.col(i).dot(nested.col(j)), static_cast<Scalar>(1), prec)) if(!ei_isMuchSmallerThan(nested.col(i).dot(nested.col(j)), static_cast<Scalar>(1), prec))
return false; return false;
} }

View File

@ -202,7 +202,7 @@ struct ei_fuzzy_selector<Derived,OtherDerived,false>
ei_assert(self.rows() == other.rows() && self.cols() == other.cols()); ei_assert(self.rows() == other.rows() && self.cols() == other.cols());
typename Derived::Nested nested(self); typename Derived::Nested nested(self);
typename OtherDerived::Nested otherNested(other); typename OtherDerived::Nested otherNested(other);
for(int i = 0; i < self.cols(); i++) for(int i = 0; i < self.cols(); ++i)
if((nested.col(i) - otherNested.col(i)).squaredNorm() if((nested.col(i) - otherNested.col(i)).squaredNorm()
> std::min(nested.col(i).squaredNorm(), otherNested.col(i).squaredNorm()) * prec * prec) > std::min(nested.col(i).squaredNorm(), otherNested.col(i).squaredNorm()) * prec * prec)
return false; return false;
@ -211,7 +211,7 @@ struct ei_fuzzy_selector<Derived,OtherDerived,false>
static bool isMuchSmallerThan(const Derived& self, const RealScalar& other, RealScalar prec) static bool isMuchSmallerThan(const Derived& self, const RealScalar& other, RealScalar prec)
{ {
typename Derived::Nested nested(self); typename Derived::Nested nested(self);
for(int i = 0; i < self.cols(); i++) for(int i = 0; i < self.cols(); ++i)
if(nested.col(i).squaredNorm() > ei_abs2(other * prec)) if(nested.col(i).squaredNorm() > ei_abs2(other * prec))
return false; return false;
return true; return true;
@ -222,7 +222,7 @@ struct ei_fuzzy_selector<Derived,OtherDerived,false>
ei_assert(self.rows() == other.rows() && self.cols() == other.cols()); ei_assert(self.rows() == other.rows() && self.cols() == other.cols());
typename Derived::Nested nested(self); typename Derived::Nested nested(self);
typename OtherDerived::Nested otherNested(other); typename OtherDerived::Nested otherNested(other);
for(int i = 0; i < self.cols(); i++) for(int i = 0; i < self.cols(); ++i)
if(nested.col(i).squaredNorm() > otherNested.col(i).squaredNorm() * prec * prec) if(nested.col(i).squaredNorm() > otherNested.col(i).squaredNorm() * prec * prec)
return false; return false;
return true; return true;

View File

@ -129,8 +129,8 @@ std::ostream & ei_print_matrix(std::ostream & s, const MatrixBase<Derived> & _m,
if (fmt.flags & AlignCols) if (fmt.flags & AlignCols)
{ {
// compute the largest width // compute the largest width
for(int j = 1; j < m.cols(); j++) for(int j = 1; j < m.cols(); ++j)
for(int i = 0; i < m.rows(); i++) for(int i = 0; i < m.rows(); ++i)
{ {
std::stringstream sstr; std::stringstream sstr;
sstr.precision(fmt.precision); sstr.precision(fmt.precision);
@ -140,14 +140,14 @@ std::ostream & ei_print_matrix(std::ostream & s, const MatrixBase<Derived> & _m,
} }
s.precision(fmt.precision); s.precision(fmt.precision);
s << fmt.matPrefix; s << fmt.matPrefix;
for(int i = 0; i < m.rows(); i++) for(int i = 0; i < m.rows(); ++i)
{ {
if (i) if (i)
s << fmt.rowSpacer; s << fmt.rowSpacer;
s << fmt.rowPrefix; s << fmt.rowPrefix;
if(width) s.width(width); if(width) s.width(width);
s << m.coeff(i, 0); s << m.coeff(i, 0);
for(int j = 1; j < m.cols(); j++) for(int j = 1; j < m.cols(); ++j)
{ {
s << fmt.coeffSeparator; s << fmt.coeffSeparator;
if (width) s.width(width); if (width) s.width(width);

View File

@ -219,8 +219,8 @@ struct ei_part_assignment_impl<Derived1, Derived2, Upper, Dynamic>
{ {
inline static void run(Derived1 &dst, const Derived2 &src) inline static void run(Derived1 &dst, const Derived2 &src)
{ {
for(int j = 0; j < dst.cols(); j++) for(int j = 0; j < dst.cols(); ++j)
for(int i = 0; i <= j; i++) for(int i = 0; i <= j; ++i)
dst.copyCoeff(i, j, src); dst.copyCoeff(i, j, src);
} }
}; };
@ -230,8 +230,8 @@ struct ei_part_assignment_impl<Derived1, Derived2, Lower, Dynamic>
{ {
inline static void run(Derived1 &dst, const Derived2 &src) inline static void run(Derived1 &dst, const Derived2 &src)
{ {
for(int j = 0; j < dst.cols(); j++) for(int j = 0; j < dst.cols(); ++j)
for(int i = j; i < dst.rows(); i++) for(int i = j; i < dst.rows(); ++i)
dst.copyCoeff(i, j, src); dst.copyCoeff(i, j, src);
} }
}; };
@ -241,8 +241,8 @@ struct ei_part_assignment_impl<Derived1, Derived2, StrictlyUpper, Dynamic>
{ {
inline static void run(Derived1 &dst, const Derived2 &src) inline static void run(Derived1 &dst, const Derived2 &src)
{ {
for(int j = 0; j < dst.cols(); j++) for(int j = 0; j < dst.cols(); ++j)
for(int i = 0; i < j; i++) for(int i = 0; i < j; ++i)
dst.copyCoeff(i, j, src); dst.copyCoeff(i, j, src);
} }
}; };
@ -251,8 +251,8 @@ struct ei_part_assignment_impl<Derived1, Derived2, StrictlyLower, Dynamic>
{ {
inline static void run(Derived1 &dst, const Derived2 &src) inline static void run(Derived1 &dst, const Derived2 &src)
{ {
for(int j = 0; j < dst.cols(); j++) for(int j = 0; j < dst.cols(); ++j)
for(int i = j+1; i < dst.rows(); i++) for(int i = j+1; i < dst.rows(); ++i)
dst.copyCoeff(i, j, src); dst.copyCoeff(i, j, src);
} }
}; };
@ -261,9 +261,9 @@ struct ei_part_assignment_impl<Derived1, Derived2, SelfAdjoint, Dynamic>
{ {
inline static void run(Derived1 &dst, const Derived2 &src) inline static void run(Derived1 &dst, const Derived2 &src)
{ {
for(int j = 0; j < dst.cols(); j++) for(int j = 0; j < dst.cols(); ++j)
{ {
for(int i = 0; i < j; i++) for(int i = 0; i < j; ++i)
dst.coeffRef(j, i) = ei_conj(dst.coeffRef(i, j) = src.coeff(i, j)); dst.coeffRef(j, i) = ei_conj(dst.coeffRef(i, j) = src.coeff(i, j));
dst.coeffRef(j, j) = ei_real(src.coeff(j, j)); dst.coeffRef(j, j) = ei_real(src.coeff(j, j));
} }
@ -312,14 +312,14 @@ bool MatrixBase<Derived>::isUpper(RealScalar prec) const
{ {
if(cols() != rows()) return false; if(cols() != rows()) return false;
RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1); RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
for(int j = 0; j < cols(); j++) for(int j = 0; j < cols(); ++j)
for(int i = 0; i <= j; i++) for(int i = 0; i <= j; ++i)
{ {
RealScalar absValue = ei_abs(coeff(i,j)); RealScalar absValue = ei_abs(coeff(i,j));
if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue; if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
} }
for(int j = 0; j < cols()-1; j++) for(int j = 0; j < cols()-1; ++j)
for(int i = j+1; i < rows(); i++) for(int i = j+1; i < rows(); ++i)
if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperPart, prec)) return false; if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperPart, prec)) return false;
return true; return true;
} }
@ -334,14 +334,14 @@ bool MatrixBase<Derived>::isLower(RealScalar prec) const
{ {
if(cols() != rows()) return false; if(cols() != rows()) return false;
RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1); RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
for(int j = 0; j < cols(); j++) for(int j = 0; j < cols(); ++j)
for(int i = j; i < rows(); i++) for(int i = j; i < rows(); ++i)
{ {
RealScalar absValue = ei_abs(coeff(i,j)); RealScalar absValue = ei_abs(coeff(i,j));
if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue; if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
} }
for(int j = 1; j < cols(); j++) for(int j = 1; j < cols(); ++j)
for(int i = 0; i < j; i++) for(int i = 0; i < j; ++i)
if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerPart, prec)) return false; if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerPart, prec)) return false;
return true; return true;
} }

View File

@ -337,7 +337,7 @@ struct ei_product_coeff_impl<NoVectorization, Dynamic, Lhs, Rhs>
{ {
ei_assert(lhs.cols()>0 && "you are using a non initialized matrix"); ei_assert(lhs.cols()>0 && "you are using a non initialized matrix");
res = lhs.coeff(row, 0) * rhs.coeff(0, col); res = lhs.coeff(row, 0) * rhs.coeff(0, col);
for(int i = 1; i < lhs.cols(); i++) for(int i = 1; i < lhs.cols(); ++i)
res += lhs.coeff(row, i) * rhs.coeff(i, col); res += lhs.coeff(row, i) * rhs.coeff(i, col);
} }
}; };
@ -495,7 +495,7 @@ struct ei_product_packet_impl<RowMajor, Dynamic, Lhs, Rhs, PacketScalar, LoadMod
{ {
ei_assert(lhs.cols()>0 && "you are using a non initialized matrix"); ei_assert(lhs.cols()>0 && "you are using a non initialized matrix");
res = ei_pmul(ei_pset1(lhs.coeff(row, 0)),rhs.template packet<LoadMode>(0, col)); res = ei_pmul(ei_pset1(lhs.coeff(row, 0)),rhs.template packet<LoadMode>(0, col));
for(int i = 1; i < lhs.cols(); i++) for(int i = 1; i < lhs.cols(); ++i)
res = ei_pmadd(ei_pset1(lhs.coeff(row, i)), rhs.template packet<LoadMode>(i, col), res); res = ei_pmadd(ei_pset1(lhs.coeff(row, i)), rhs.template packet<LoadMode>(i, col), res);
} }
}; };
@ -507,7 +507,7 @@ struct ei_product_packet_impl<ColMajor, Dynamic, Lhs, Rhs, PacketScalar, LoadMod
{ {
ei_assert(lhs.cols()>0 && "you are using a non initialized matrix"); ei_assert(lhs.cols()>0 && "you are using a non initialized matrix");
res = ei_pmul(lhs.template packet<LoadMode>(row, 0), ei_pset1(rhs.coeff(0, col))); res = ei_pmul(lhs.template packet<LoadMode>(row, 0), ei_pset1(rhs.coeff(0, col)));
for(int i = 1; i < lhs.cols(); i++) for(int i = 1; i < lhs.cols(); ++i)
res = ei_pmadd(lhs.template packet<LoadMode>(row, i), ei_pset1(rhs.coeff(i, col)), res); res = ei_pmadd(lhs.template packet<LoadMode>(row, i), ei_pset1(rhs.coeff(i, col)), res);
} }
}; };

View File

@ -68,10 +68,10 @@ struct ei_redux_impl<BinaryOp, Derived, Start, Dynamic>
ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using a non initialized matrix"); ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using a non initialized matrix");
Scalar res; Scalar res;
res = mat.coeff(0,0); res = mat.coeff(0,0);
for(int i = 1; i < mat.rows(); i++) for(int i = 1; i < mat.rows(); ++i)
res = func(res, mat.coeff(i, 0)); res = func(res, mat.coeff(i, 0));
for(int j = 1; j < mat.cols(); j++) for(int j = 1; j < mat.cols(); ++j)
for(int i = 0; i < mat.rows(); i++) for(int i = 0; i < mat.rows(); ++i)
res = func(res, mat.coeff(i, j)); res = func(res, mat.coeff(i, j));
return res; return res;
} }

View File

@ -168,10 +168,10 @@ struct ei_sum_impl<Derived, NoVectorization, NoUnrolling>
ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using a non initialized matrix"); ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using a non initialized matrix");
Scalar res; Scalar res;
res = mat.coeff(0, 0); res = mat.coeff(0, 0);
for(int i = 1; i < mat.rows(); i++) for(int i = 1; i < mat.rows(); ++i)
res += mat.coeff(i, 0); res += mat.coeff(i, 0);
for(int j = 1; j < mat.cols(); j++) for(int j = 1; j < mat.cols(); ++j)
for(int i = 0; i < mat.rows(); i++) for(int i = 0; i < mat.rows(); ++i)
res += mat.coeff(i, j); res += mat.coeff(i, j);
return res; return res;
} }
@ -217,10 +217,10 @@ struct ei_sum_impl<Derived, LinearVectorization, NoUnrolling>
res = Scalar(0); res = Scalar(0);
} }
for(int index = 0; index < alignedStart; index++) for(int index = 0; index < alignedStart; ++index)
res += mat.coeff(index); res += mat.coeff(index);
for(int index = alignedEnd; index < size; index++) for(int index = alignedEnd; index < size; ++index)
res += mat.coeff(index); res += mat.coeff(index);
return res; return res;

View File

@ -55,10 +55,10 @@ struct ei_visitor_impl<Visitor, Derived, Dynamic>
inline static void run(const Derived& mat, Visitor& visitor) inline static void run(const Derived& mat, Visitor& visitor)
{ {
visitor.init(mat.coeff(0,0), 0, 0); visitor.init(mat.coeff(0,0), 0, 0);
for(int i = 1; i < mat.rows(); i++) for(int i = 1; i < mat.rows(); ++i)
visitor(mat.coeff(i, 0), i, 0); visitor(mat.coeff(i, 0), i, 0);
for(int j = 1; j < mat.cols(); j++) for(int j = 1; j < mat.cols(); ++j)
for(int i = 0; i < mat.rows(); i++) for(int i = 0; i < mat.rows(); ++i)
visitor(mat.coeff(i, j), i, j); visitor(mat.coeff(i, j), i, j);
} }
}; };

View File

@ -312,7 +312,7 @@ LU<MatrixType>::LU(const MatrixType& matrix)
int number_of_transpositions = 0; int number_of_transpositions = 0;
RealScalar biggest = RealScalar(0); RealScalar biggest = RealScalar(0);
for(int k = 0; k < size; k++) for(int k = 0; k < size; ++k)
{ {
int row_of_biggest_in_corner, col_of_biggest_in_corner; int row_of_biggest_in_corner, col_of_biggest_in_corner;
RealScalar biggest_in_corner; RealScalar biggest_in_corner;
@ -326,11 +326,11 @@ LU<MatrixType>::LU(const MatrixType& matrix)
cols_transpositions.coeffRef(k) = col_of_biggest_in_corner; cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
if(k != row_of_biggest_in_corner) { if(k != row_of_biggest_in_corner) {
m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner)); m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
number_of_transpositions++; ++number_of_transpositions;
} }
if(k != col_of_biggest_in_corner) { if(k != col_of_biggest_in_corner) {
m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner)); m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
number_of_transpositions++; ++number_of_transpositions;
} }
if(k==0) biggest = biggest_in_corner; if(k==0) biggest = biggest_in_corner;
@ -339,21 +339,21 @@ LU<MatrixType>::LU(const MatrixType& matrix)
if(k<rows-1) if(k<rows-1)
m_lu.col(k).end(rows-k-1) /= lu_k_k; m_lu.col(k).end(rows-k-1) /= lu_k_k;
if(k<size-1) if(k<size-1)
for( int col = k + 1; col < cols; col++ ) for(int col = k + 1; col < cols; ++col)
m_lu.col(col).end(rows-k-1) -= m_lu.col(k).end(rows-k-1) * m_lu.coeff(k,col); m_lu.col(col).end(rows-k-1) -= m_lu.col(k).end(rows-k-1) * m_lu.coeff(k,col);
} }
for(int k = 0; k < matrix.rows(); k++) m_p.coeffRef(k) = k; for(int k = 0; k < matrix.rows(); ++k) m_p.coeffRef(k) = k;
for(int k = size-1; k >= 0; k--) for(int k = size-1; k >= 0; k--)
std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k))); std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k)));
for(int k = 0; k < matrix.cols(); k++) m_q.coeffRef(k) = k; for(int k = 0; k < matrix.cols(); ++k) m_q.coeffRef(k) = k;
for(int k = 0; k < size; k++) for(int k = 0; k < size; ++k)
std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k))); std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k)));
m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_det_pq = (number_of_transpositions%2) ? -1 : 1;
for(m_rank = 0; m_rank < size; m_rank++) for(m_rank = 0; m_rank < size; ++m_rank)
if(ei_isMuchSmallerThan(m_lu.diagonal().coeff(m_rank), m_lu.diagonal().coeff(0))) if(ei_isMuchSmallerThan(m_lu.diagonal().coeff(m_rank), m_lu.diagonal().coeff(0)))
break; break;
} }
@ -395,10 +395,10 @@ void LU<MatrixType>::computeKernel(KernelResultType *result) const
.template marked<Upper>() .template marked<Upper>()
.solveTriangularInPlace(y); .solveTriangularInPlace(y);
for(int i = 0; i < m_rank; i++) for(int i = 0; i < m_rank; ++i)
result->row(m_q.coeff(i)) = y.row(i); result->row(m_q.coeff(i)) = y.row(i);
for(int i = m_rank; i < cols; i++) result->row(m_q.coeff(i)).setZero(); for(int i = m_rank; i < cols; ++i) result->row(m_q.coeff(i)).setZero();
for(int k = 0; k < dimker; k++) result->coeffRef(m_q.coeff(m_rank+k), k) = Scalar(1); for(int k = 0; k < dimker; ++k) result->coeffRef(m_q.coeff(m_rank+k), k) = Scalar(1);
} }
template<typename MatrixType> template<typename MatrixType>
@ -432,7 +432,7 @@ bool LU<MatrixType>::solve(
typename OtherDerived::Eval c(b.rows(), b.cols()); typename OtherDerived::Eval c(b.rows(), b.cols());
// Step 1 // Step 1
for(int i = 0; i < rows; i++) c.row(m_p.coeff(i)) = b.row(i); for(int i = 0; i < rows; ++i) c.row(m_p.coeff(i)) = b.row(i);
// Step 2 // Step 2
Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime, Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime,
@ -449,8 +449,8 @@ bool LU<MatrixType>::solve(
{ {
// is c is in the image of U ? // is c is in the image of U ?
RealScalar biggest_in_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff(); RealScalar biggest_in_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
for(int col = 0; col < c.cols(); col++) for(int col = 0; col < c.cols(); ++col)
for(int row = m_rank; row < c.rows(); row++) for(int row = m_rank; row < c.rows(); ++row)
if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c)) if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c))
return false; return false;
} }
@ -464,8 +464,8 @@ bool LU<MatrixType>::solve(
// Step 4 // Step 4
result->resize(m_lu.cols(), b.cols()); result->resize(m_lu.cols(), b.cols());
for(int i = 0; i < m_rank; i++) result->row(m_q.coeff(i)) = d.row(i); for(int i = 0; i < m_rank; ++i) result->row(m_q.coeff(i)) = d.row(i);
for(int i = m_rank; i < m_lu.cols(); i++) result->row(m_q.coeff(i)).setZero(); for(int i = m_rank; i < m_lu.cols(); ++i) result->row(m_q.coeff(i)).setZero();
return true; return true;
} }

View File

@ -122,7 +122,7 @@ MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
{ {
int n = m_eivec.cols(); int n = m_eivec.cols();
MatrixType matD = MatrixType::Zero(n,n); MatrixType matD = MatrixType::Zero(n,n);
for (int i=0; i<n; i++) for (int i=0; i<n; ++i)
{ {
if (ei_isMuchSmallerThan(ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i)))) if (ei_isMuchSmallerThan(ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i))))
matD.coeffRef(i,i) = ei_real(m_eivalues.coeff(i)); matD.coeffRef(i,i) = ei_real(m_eivalues.coeff(i));
@ -130,7 +130,7 @@ MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
{ {
matD.template block<2,2>(i,i) << ei_real(m_eivalues.coeff(i)), ei_imag(m_eivalues.coeff(i)), matD.template block<2,2>(i,i) << ei_real(m_eivalues.coeff(i)), ei_imag(m_eivalues.coeff(i)),
-ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i)); -ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i));
i++; ++i;
} }
} }
return matD; return matD;
@ -145,7 +145,7 @@ typename EigenSolver<MatrixType>::EigenvectorType EigenSolver<MatrixType>::eigen
{ {
int n = m_eivec.cols(); int n = m_eivec.cols();
EigenvectorType matV(n,n); EigenvectorType matV(n,n);
for (int j=0; j<n; j++) for (int j=0; j<n; ++j)
{ {
if (ei_isMuchSmallerThan(ei_abs(ei_imag(m_eivalues.coeff(j))), ei_abs(ei_real(m_eivalues.coeff(j))))) if (ei_isMuchSmallerThan(ei_abs(ei_imag(m_eivalues.coeff(j))), ei_abs(ei_real(m_eivalues.coeff(j)))))
{ {
@ -155,14 +155,14 @@ typename EigenSolver<MatrixType>::EigenvectorType EigenSolver<MatrixType>::eigen
else else
{ {
// we have a pair of complex eigen values // we have a pair of complex eigen values
for (int i=0; i<n; i++) for (int i=0; i<n; ++i)
{ {
matV.coeffRef(i,j) = Complex(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1)); matV.coeffRef(i,j) = Complex(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
matV.coeffRef(i,j+1) = Complex(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1)); matV.coeffRef(i,j+1) = Complex(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
} }
matV.col(j).normalize(); matV.col(j).normalize();
matV.col(j+1).normalize(); matV.col(j+1).normalize();
j++; ++j;
} }
} }
return matV; return matV;
@ -198,7 +198,7 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
int low = 0; int low = 0;
int high = n-1; int high = n-1;
for (int m = low+1; m <= high-1; m++) for (int m = low+1; m <= high-1; ++m)
{ {
// Scale column. // Scale column.
RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum(); RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum();
@ -290,7 +290,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
// FIXME to be efficient the following would requires a triangular reduxion code // FIXME to be efficient the following would requires a triangular reduxion code
// Scalar norm = matH.upper().cwise().abs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwise().abs().sum(); // Scalar norm = matH.upper().cwise().abs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwise().abs().sum();
Scalar norm = 0.0; Scalar norm = 0.0;
for (int j = 0; j < nn; j++) for (int j = 0; j < nn; ++j)
{ {
// FIXME what's the purpose of the following since the condition is always false // FIXME what's the purpose of the following since the condition is always false
if ((j < low) || (j > high)) if ((j < low) || (j > high))
@ -361,7 +361,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
q = q / r; q = q / r;
// Row modification // Row modification
for (int j = n-1; j < nn; j++) for (int j = n-1; j < nn; ++j)
{ {
z = matH.coeff(n-1,j); z = matH.coeff(n-1,j);
matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j); matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j);
@ -369,7 +369,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
} }
// Column modification // Column modification
for (int i = 0; i <= n; i++) for (int i = 0; i <= n; ++i)
{ {
z = matH.coeff(i,n-1); z = matH.coeff(i,n-1);
matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n); matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n);
@ -377,7 +377,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
} }
// Accumulate transformations // Accumulate transformations
for (int i = low; i <= high; i++) for (int i = low; i <= high; ++i)
{ {
z = m_eivec.coeff(i,n-1); z = m_eivec.coeff(i,n-1);
m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n); m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n);
@ -410,7 +410,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
if (iter == 10) if (iter == 10)
{ {
exshift += x; exshift += x;
for (int i = low; i <= n; i++) for (int i = low; i <= n; ++i)
matH.coeffRef(i,i) -= x; matH.coeffRef(i,i) -= x;
s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2)); s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2));
x = y = 0.75 * s; x = y = 0.75 * s;
@ -428,7 +428,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
if (y < x) if (y < x)
s = -s; s = -s;
s = x - w / ((y - x) / 2.0 + s); s = x - w / ((y - x) / 2.0 + s);
for (int i = low; i <= n; i++) for (int i = low; i <= n; ++i)
matH.coeffRef(i,i) -= s; matH.coeffRef(i,i) -= s;
exshift += s; exshift += s;
x = y = w = 0.964; x = y = w = 0.964;
@ -463,7 +463,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
m--; m--;
} }
for (int i = m+2; i <= n; i++) for (int i = m+2; i <= n; ++i)
{ {
matH.coeffRef(i,i-2) = 0.0; matH.coeffRef(i,i-2) = 0.0;
if (i > m+2) if (i > m+2)
@ -471,7 +471,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
} }
// Double QR step involving rows l:n and columns m:n // Double QR step involving rows l:n and columns m:n
for (int k = m; k <= n-1; k++) for (int k = m; k <= n-1; ++k)
{ {
int notlast = (k != n-1); int notlast = (k != n-1);
if (k != m) { if (k != m) {
@ -510,7 +510,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
r = r / p; r = r / p;
// Row modification // Row modification
for (int j = k; j < nn; j++) for (int j = k; j < nn; ++j)
{ {
p = matH.coeff(k,j) + q * matH.coeff(k+1,j); p = matH.coeff(k,j) + q * matH.coeff(k+1,j);
if (notlast) if (notlast)
@ -523,7 +523,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
} }
// Column modification // Column modification
for (int i = 0; i <= std::min(n,k+3); i++) for (int i = 0; i <= std::min(n,k+3); ++i)
{ {
p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1); p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1);
if (notlast) if (notlast)
@ -536,7 +536,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
} }
// Accumulate transformations // Accumulate transformations
for (int i = low; i <= high; i++) for (int i = low; i <= high; ++i)
{ {
p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1); p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1);
if (notlast) if (notlast)
@ -686,7 +686,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
} }
// Vectors of isolated roots // Vectors of isolated roots
for (int i = 0; i < nn; i++) for (int i = 0; i < nn; ++i)
{ {
// FIXME again what's the purpose of this test ? // FIXME again what's the purpose of this test ?
// in this algo low==0 and high==nn-1 !! // in this algo low==0 and high==nn-1 !!

View File

@ -87,7 +87,7 @@ void QR<MatrixType>::_compute(const MatrixType& matrix)
int rows = matrix.rows(); int rows = matrix.rows();
int cols = matrix.cols(); int cols = matrix.cols();
for (int k = 0; k < cols; k++) for (int k = 0; k < cols; ++k)
{ {
int remainingSize = rows-k; int remainingSize = rows-k;

View File

@ -202,7 +202,7 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool
// Sort eigenvalues and corresponding vectors. // Sort eigenvalues and corresponding vectors.
// TODO make the sort optional ? // TODO make the sort optional ?
// TODO use a better sort algorithm !! // TODO use a better sort algorithm !!
for (int i = 0; i < n-1; i++) for (int i = 0; i < n-1; ++i)
{ {
int k; int k;
m_eivalues.segment(i,n-i).minCoeff(&k); m_eivalues.segment(i,n-i).minCoeff(&k);

View File

@ -268,7 +268,7 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
* if we remove the specialization of Block for Matrix then it is even worse, much worse ! */ * if we remove the specialization of Block for Matrix then it is even worse, much worse ! */
#ifdef EIGEN_NEVER_DEFINED #ifdef EIGEN_NEVER_DEFINED
for (int j1=i+1; j1<n; ++j1) for (int j1=i+1; j1<n; ++j1)
for (int i1=j1; i1<n; i1++) for (int i1=j1; i1<n; ++i1)
matA.coeffRef(i1,j1) -= matA.coeff(i1,i)*ei_conj(hCoeffs.coeff(j1-1)) matA.coeffRef(i1,j1) -= matA.coeff(i1,i)*ei_conj(hCoeffs.coeff(j1-1))
+ hCoeffs.coeff(i1-1)*ei_conj(matA.coeff(j1,i)); + hCoeffs.coeff(i1-1)*ei_conj(matA.coeff(j1,i));
#endif #endif

View File

@ -111,18 +111,18 @@ void linearRegression(int numPoints,
Dynamic, VectorType::MaxSizeAtCompileTime, RowMajorBit> Dynamic, VectorType::MaxSizeAtCompileTime, RowMajorBit>
m(numPoints, size); m(numPoints, size);
if(funcOfOthers>0) if(funcOfOthers>0)
for(int i = 0; i < numPoints; i++) for(int i = 0; i < numPoints; ++i)
m.row(i).start(funcOfOthers) = points[i]->start(funcOfOthers); m.row(i).start(funcOfOthers) = points[i]->start(funcOfOthers);
if(funcOfOthers<size-1) if(funcOfOthers<size-1)
for(int i = 0; i < numPoints; i++) for(int i = 0; i < numPoints; ++i)
m.row(i).block(funcOfOthers, size-funcOfOthers-1) m.row(i).block(funcOfOthers, size-funcOfOthers-1)
= points[i]->end(size-funcOfOthers-1); = points[i]->end(size-funcOfOthers-1);
for(int i = 0; i < numPoints; i++) for(int i = 0; i < numPoints; ++i)
m.row(i).coeffRef(size-1) = Scalar(1); m.row(i).coeffRef(size-1) = Scalar(1);
VectorType v(size); VectorType v(size);
v.setZero(); v.setZero();
for(int i = 0; i < numPoints; i++) for(int i = 0; i < numPoints; ++i)
v += m.row(i).adjoint() * points[i]->coeff(funcOfOthers); v += m.row(i).adjoint() * points[i]->coeff(funcOfOthers);
ei_assert((m.adjoint()*m).lu().solve(v, result)); ei_assert((m.adjoint()*m).lu().solve(v, result));
@ -170,14 +170,14 @@ void fitHyperplane(int numPoints,
// compute the mean of the data // compute the mean of the data
VectorType mean = VectorType::Zero(size); VectorType mean = VectorType::Zero(size);
for(int i = 0; i < numPoints; i++) for(int i = 0; i < numPoints; ++i)
mean += *(points[i]); mean += *(points[i]);
mean /= numPoints; mean /= numPoints;
// compute the covariance matrix // compute the covariance matrix
CovMatrixType covMat = CovMatrixType::Zero(size, size); CovMatrixType covMat = CovMatrixType::Zero(size, size);
VectorType remean = VectorType::Zero(size); VectorType remean = VectorType::Zero(size);
for(int i = 0; i < numPoints; i++) for(int i = 0; i < numPoints; ++i)
{ {
VectorType diff = (*(points[i]) - mean).conjugate(); VectorType diff = (*(points[i]) - mean).conjugate();
covMat += diff * diff.adjoint(); covMat += diff * diff.adjoint();

View File

@ -115,7 +115,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
// in s and the super-diagonal elements in e. // in s and the super-diagonal elements in e.
int nct = std::min(m-1,n); int nct = std::min(m-1,n);
int nrt = std::max(0,std::min(n-2,m)); int nrt = std::max(0,std::min(n-2,m));
for (k = 0; k < std::max(nct,nrt); k++) for (k = 0; k < std::max(nct,nrt); ++k)
{ {
if (k < nct) if (k < nct)
{ {
@ -132,7 +132,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
m_sigma[k] = -m_sigma[k]; m_sigma[k] = -m_sigma[k];
} }
for (j = k+1; j < n; j++) for (j = k+1; j < n; ++j)
{ {
if ((k < nct) && (m_sigma[k] != 0.0)) if ((k < nct) && (m_sigma[k] != 0.0))
{ {
@ -168,7 +168,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
{ {
// Apply the transformation. // Apply the transformation.
work.end(m-k-1) = matA.corner(BottomRight,m-k-1,n-k-1) * e.end(n-k-1); work.end(m-k-1) = matA.corner(BottomRight,m-k-1,n-k-1) * e.end(n-k-1);
for (j = k+1; j < n; j++) for (j = k+1; j < n; ++j)
matA.col(j).end(m-k-1) += (-e[j]/e[k+1]) * work.end(m-k-1); matA.col(j).end(m-k-1) += (-e[j]/e[k+1]) * work.end(m-k-1);
} }
@ -192,7 +192,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
// If required, generate U. // If required, generate U.
if (wantu) if (wantu)
{ {
for (j = nct; j < nu; j++) for (j = nct; j < nu; ++j)
{ {
m_matU.col(j).setZero(); m_matU.col(j).setZero();
m_matU(j,j) = 1.0; m_matU(j,j) = 1.0;
@ -201,7 +201,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
{ {
if (m_sigma[k] != 0.0) if (m_sigma[k] != 0.0)
{ {
for (j = k+1; j < nu; j++) for (j = k+1; j < nu; ++j)
{ {
Scalar t = m_matU.col(k).end(m-k).dot(m_matU.col(j).end(m-k)); // FIXME is it really a dot product we want ? Scalar t = m_matU.col(k).end(m-k).dot(m_matU.col(j).end(m-k)); // FIXME is it really a dot product we want ?
t = -t/m_matU(k,k); t = -t/m_matU(k,k);
@ -227,7 +227,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
{ {
if ((k < nrt) & (e[k] != 0.0)) if ((k < nrt) & (e[k] != 0.0))
{ {
for (j = k+1; j < nu; j++) for (j = k+1; j < nu; ++j)
{ {
Scalar t = m_matV.col(k).end(n-k-1).dot(m_matV.col(j).end(n-k-1)); // FIXME is it really a dot product we want ? Scalar t = m_matV.col(k).end(n-k-1).dot(m_matV.col(j).end(n-k-1)); // FIXME is it really a dot product we want ?
t = -t/m_matV(k+1,k); t = -t/m_matV(k+1,k);
@ -302,7 +302,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
k = ks; k = ks;
} }
} }
k++; ++k;
// Perform the task indicated by kase. // Perform the task indicated by kase.
switch (kase) switch (kase)
@ -326,7 +326,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
} }
if (wantv) if (wantv)
{ {
for (i = 0; i < n; i++) for (i = 0; i < n; ++i)
{ {
t = cs*m_matV(i,j) + sn*m_matV(i,p-1); t = cs*m_matV(i,j) + sn*m_matV(i,p-1);
m_matV(i,p-1) = -sn*m_matV(i,j) + cs*m_matV(i,p-1); m_matV(i,p-1) = -sn*m_matV(i,j) + cs*m_matV(i,p-1);
@ -342,7 +342,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
{ {
Scalar f(e[k-1]); Scalar f(e[k-1]);
e[k-1] = 0.0; e[k-1] = 0.0;
for (j = k; j < p; j++) for (j = k; j < p; ++j)
{ {
Scalar t(hypot(m_sigma[j],f)); Scalar t(hypot(m_sigma[j],f));
Scalar cs( m_sigma[j]/t); Scalar cs( m_sigma[j]/t);
@ -352,7 +352,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
e[j] = cs*e[j]; e[j] = cs*e[j];
if (wantu) if (wantu)
{ {
for (i = 0; i < m; i++) for (i = 0; i < m; ++i)
{ {
t = cs*m_matU(i,j) + sn*m_matU(i,k-1); t = cs*m_matU(i,j) + sn*m_matU(i,k-1);
m_matU(i,k-1) = -sn*m_matU(i,j) + cs*m_matU(i,k-1); m_matU(i,k-1) = -sn*m_matU(i,j) + cs*m_matU(i,k-1);
@ -390,7 +390,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
// Chase zeros. // Chase zeros.
for (j = k; j < p-1; j++) for (j = k; j < p-1; ++j)
{ {
Scalar t = hypot(f,g); Scalar t = hypot(f,g);
Scalar cs = f/t; Scalar cs = f/t;
@ -403,7 +403,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
m_sigma[j+1] = cs*m_sigma[j+1]; m_sigma[j+1] = cs*m_sigma[j+1];
if (wantv) if (wantv)
{ {
for (i = 0; i < n; i++) for (i = 0; i < n; ++i)
{ {
t = cs*m_matV(i,j) + sn*m_matV(i,j+1); t = cs*m_matV(i,j) + sn*m_matV(i,j+1);
m_matV(i,j+1) = -sn*m_matV(i,j) + cs*m_matV(i,j+1); m_matV(i,j+1) = -sn*m_matV(i,j) + cs*m_matV(i,j+1);
@ -420,7 +420,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
e[j+1] = cs*e[j+1]; e[j+1] = cs*e[j+1];
if (wantu && (j < m-1)) if (wantu && (j < m-1))
{ {
for (i = 0; i < m; i++) for (i = 0; i < m; ++i)
{ {
t = cs*m_matU(i,j) + sn*m_matU(i,j+1); t = cs*m_matU(i,j) + sn*m_matU(i,j+1);
m_matU(i,j+1) = -sn*m_matU(i,j) + cs*m_matU(i,j+1); m_matU(i,j+1) = -sn*m_matU(i,j) + cs*m_matU(i,j+1);
@ -456,7 +456,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
m_matV.col(k).swap(m_matV.col(k+1)); m_matV.col(k).swap(m_matV.col(k+1));
if (wantu && (k < m-1)) if (wantu && (k < m-1))
m_matU.col(k).swap(m_matU.col(k+1)); m_matU.col(k).swap(m_matU.col(k+1));
k++; ++k;
} }
iter = 0; iter = 0;
p--; p--;
@ -473,12 +473,12 @@ SVD<MatrixType>& SVD<MatrixType>::sort()
int mv = m_matV.rows(); int mv = m_matV.rows();
int n = m_matU.cols(); int n = m_matU.cols();
for (int i=0; i<n; i++) for (int i=0; i<n; ++i)
{ {
int k = i; int k = i;
Scalar p = m_sigma.coeff(i); Scalar p = m_sigma.coeff(i);
for (int j=i+1; j<n; j++) for (int j=i+1; j<n; ++j)
{ {
if (m_sigma.coeff(j) > p) if (m_sigma.coeff(j) > p)
{ {
@ -520,7 +520,7 @@ bool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* resul
{ {
Matrix<Scalar,MatrixUType::RowsAtCompileTime,1> aux = m_matU.transpose() * b.col(j); Matrix<Scalar,MatrixUType::RowsAtCompileTime,1> aux = m_matU.transpose() * b.col(j);
for (int i = 0; i <m_matU.cols(); i++) for (int i = 0; i <m_matU.cols(); ++i)
{ {
Scalar si = m_sigma.coeff(i); Scalar si = m_sigma.coeff(i);
if (ei_isMuchSmallerThan(ei_abs(si),maxVal)) if (ei_isMuchSmallerThan(ei_abs(si),maxVal))

View File

@ -202,7 +202,7 @@ Scalar& AmbiVector<Scalar>::coeffRef(int i)
// this is the first element // this is the first element
m_llStart = 0; m_llStart = 0;
m_llCurrent = 0; m_llCurrent = 0;
m_llSize++; ++m_llSize;
llElements[0].value = Scalar(0); llElements[0].value = Scalar(0);
llElements[0].index = i; llElements[0].index = i;
llElements[0].next = -1; llElements[0].next = -1;
@ -216,7 +216,7 @@ Scalar& AmbiVector<Scalar>::coeffRef(int i)
el.index = i; el.index = i;
el.next = m_llStart; el.next = m_llStart;
m_llStart = m_llSize; m_llStart = m_llSize;
m_llSize++; ++m_llSize;
m_llCurrent = m_llStart; m_llCurrent = m_llStart;
return el.value; return el.value;
} }
@ -246,7 +246,7 @@ Scalar& AmbiVector<Scalar>::coeffRef(int i)
el.index = i; el.index = i;
el.next = llElements[m_llCurrent].next; el.next = llElements[m_llCurrent].next;
llElements[m_llCurrent].next = m_llSize; llElements[m_llCurrent].next = m_llSize;
m_llSize++; ++m_llSize;
return el.value; return el.value;
} }
} }
@ -332,7 +332,7 @@ class AmbiVector<_Scalar>::Iterator
if (m_isDense) if (m_isDense)
{ {
do { do {
m_cachedIndex++; ++m_cachedIndex;
} while (m_cachedIndex<m_vector.m_end && ei_abs(m_vector.m_buffer[m_cachedIndex])<m_epsilon); } while (m_cachedIndex<m_vector.m_end && ei_abs(m_vector.m_buffer[m_cachedIndex])<m_epsilon);
if (m_cachedIndex<m_vector.m_end) if (m_cachedIndex<m_vector.m_end)
m_cachedValue = m_vector.m_buffer[m_cachedIndex]; m_cachedValue = m_vector.m_buffer[m_cachedIndex];

View File

@ -140,7 +140,7 @@ class RandomSetter
m_keyBitsOffset = 0; m_keyBitsOffset = 0;
while (aux) while (aux)
{ {
m_keyBitsOffset++; ++m_keyBitsOffset;
aux = aux >> 1; aux = aux >> 1;
} }
KeyType ik = (1<<(OuterPacketBits+m_keyBitsOffset)); KeyType ik = (1<<(OuterPacketBits+m_keyBitsOffset));
@ -183,7 +183,7 @@ class RandomSetter
for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it) for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
{ {
const int outer = it->first & keyBitsMask; const int outer = it->first & keyBitsMask;
positions[outer]++; ++positions[outer];
} }
} }
// prefix sum // prefix sum

View File

@ -203,10 +203,10 @@ void SparseLDLT<MatrixType,Backend>::_symbolic(const MatrixType& a)
if (P) if (P)
{ {
/* If P is present then compute Pinv, the inverse of P */ /* If P is present then compute Pinv, the inverse of P */
for (int k = 0; k < size; k++) for (int k = 0; k < size; ++k)
Pinv[P[k]] = k; Pinv[P[k]] = k;
} }
for (int k = 0; k < size; k++) for (int k = 0; k < size; ++k)
{ {
/* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */ /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
m_parent[k] = -1; /* parent of k is not yet known */ m_parent[k] = -1; /* parent of k is not yet known */
@ -214,7 +214,7 @@ void SparseLDLT<MatrixType,Backend>::_symbolic(const MatrixType& a)
m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
int kk = P ? P[k] : k; /* kth original, or permuted, column */ int kk = P ? P[k] : k; /* kth original, or permuted, column */
int p2 = Ap[kk+1]; int p2 = Ap[kk+1];
for (int p = Ap[kk]; p < p2; p++) for (int p = Ap[kk]; p < p2; ++p)
{ {
/* A (i,k) is nonzero (original or permuted A) */ /* A (i,k) is nonzero (original or permuted A) */
int i = Pinv ? Pinv[Ai[p]] : Ai[p]; int i = Pinv ? Pinv[Ai[p]] : Ai[p];
@ -226,7 +226,7 @@ void SparseLDLT<MatrixType,Backend>::_symbolic(const MatrixType& a)
/* find parent of i if not yet determined */ /* find parent of i if not yet determined */
if (m_parent[i] == -1) if (m_parent[i] == -1)
m_parent[i] = k; m_parent[i] = k;
m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */ ++m_nonZerosPerCol[i]; /* L (k,i) is nonzero */
tags[i] = k; /* mark i as visited */ tags[i] = k; /* mark i as visited */
} }
} }
@ -234,7 +234,7 @@ void SparseLDLT<MatrixType,Backend>::_symbolic(const MatrixType& a)
} }
/* construct Lp index array from m_nonZerosPerCol column counts */ /* construct Lp index array from m_nonZerosPerCol column counts */
Lp[0] = 0; Lp[0] = 0;
for (int k = 0; k < size; k++) for (int k = 0; k < size; ++k)
Lp[k+1] = Lp[k] + m_nonZerosPerCol[k]; Lp[k+1] = Lp[k] + m_nonZerosPerCol[k];
m_matrix.resizeNonZeros(Lp[size]); m_matrix.resizeNonZeros(Lp[size]);
@ -265,7 +265,7 @@ bool SparseLDLT<MatrixType,Backend>::_numeric(const MatrixType& a)
const int* Pinv = 0; const int* Pinv = 0;
bool ok = true; bool ok = true;
for (int k = 0; k < size; k++) for (int k = 0; k < size; ++k)
{ {
/* compute nonzero pattern of kth row of L, in topological order */ /* compute nonzero pattern of kth row of L, in topological order */
y[k] = 0.0; /* Y(0:k) is now all zero */ y[k] = 0.0; /* Y(0:k) is now all zero */
@ -274,7 +274,7 @@ bool SparseLDLT<MatrixType,Backend>::_numeric(const MatrixType& a)
m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
int kk = (P) ? (P[k]) : (k); /* kth original, or permuted, column */ int kk = (P) ? (P[k]) : (k); /* kth original, or permuted, column */
int p2 = Ap[kk+1]; int p2 = Ap[kk+1];
for (int p = Ap[kk]; p < p2; p++) for (int p = Ap[kk]; p < p2; ++p)
{ {
int i = Pinv ? Pinv[Ai[p]] : Ai[p]; /* get A(i,k) */ int i = Pinv ? Pinv[Ai[p]] : Ai[p]; /* get A(i,k) */
if (i <= k) if (i <= k)
@ -293,20 +293,20 @@ bool SparseLDLT<MatrixType,Backend>::_numeric(const MatrixType& a)
/* compute numerical values kth row of L (a sparse triangular solve) */ /* compute numerical values kth row of L (a sparse triangular solve) */
m_diag[k] = y[k]; /* get D(k,k) and clear Y(k) */ m_diag[k] = y[k]; /* get D(k,k) and clear Y(k) */
y[k] = 0.0; y[k] = 0.0;
for (; top < size; top++) for (; top < size; ++top)
{ {
int i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */ int i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */
Scalar yi = y[i]; /* get and clear Y(i) */ Scalar yi = y[i]; /* get and clear Y(i) */
y[i] = 0.0; y[i] = 0.0;
int p2 = Lp[i] + m_nonZerosPerCol[i]; int p2 = Lp[i] + m_nonZerosPerCol[i];
int p; int p;
for (p = Lp[i]; p < p2; p++) for (p = Lp[i]; p < p2; ++p)
y[Li[p]] -= Lx[p] * yi; y[Li[p]] -= Lx[p] * yi;
Scalar l_ki = yi / m_diag[i]; /* the nonzero entry L(k,i) */ Scalar l_ki = yi / m_diag[i]; /* the nonzero entry L(k,i) */
m_diag[k] -= l_ki * yi; m_diag[k] -= l_ki * yi;
Li[p] = k; /* store L(k,i) in column form of L */ Li[p] = k; /* store L(k,i) in column form of L */
Lx[p] = l_ki; Lx[p] = l_ki;
m_nonZerosPerCol[i]++; /* increment count of nonzeros in col i */ ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */
} }
if (m_diag[k] == 0.0) if (m_diag[k] == 0.0)
{ {

View File

@ -163,7 +163,7 @@ class SparseMatrix
} }
assert(m_outerIndex[outer+1] == m_data.size()); assert(m_outerIndex[outer+1] == m_data.size());
int id = m_outerIndex[outer+1]; int id = m_outerIndex[outer+1];
m_outerIndex[outer+1]++; ++m_outerIndex[outer+1];
m_data.append(0, inner); m_data.append(0, inner);
return m_data.value(id); return m_data.value(id);
@ -192,7 +192,7 @@ class SparseMatrix
assert(m_outerIndex[outer+1] == m_data.size() && "invalid outer index"); assert(m_outerIndex[outer+1] == m_data.size() && "invalid outer index");
int startId = m_outerIndex[outer]; int startId = m_outerIndex[outer];
int id = m_outerIndex[outer+1]-1; int id = m_outerIndex[outer+1]-1;
m_outerIndex[outer+1]++; ++m_outerIndex[outer+1];
m_data.resize(id+2); m_data.resize(id+2);
while ( (id >= startId) && (m_data.index(id) > inner) ) while ( (id >= startId) && (m_data.index(id) > inner) )
@ -212,7 +212,7 @@ class SparseMatrix
// find the last filled column // find the last filled column
while (i>=0 && m_outerIndex[i]==0) while (i>=0 && m_outerIndex[i]==0)
--i; --i;
i++; ++i;
while (i<=m_outerSize) while (i<=m_outerSize)
{ {
m_outerIndex[i] = size; m_outerIndex[i] = size;
@ -299,7 +299,7 @@ class SparseMatrix
// FIXME the above copy could be merged with that pass // FIXME the above copy could be merged with that pass
for (int j=0; j<otherCopy.outerSize(); ++j) for (int j=0; j<otherCopy.outerSize(); ++j)
for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it) for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it)
m_outerIndex[it.index()]++; ++m_outerIndex[it.index()];
// prefix sum // prefix sum
int count = 0; int count = 0;

View File

@ -124,7 +124,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType; typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex; typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
RealScalar largerEps = 10*test_precision<RealScalar>(); // RealScalar largerEps = 10*test_precision<RealScalar>();
MatrixType a = MatrixType::Random(rows,cols); MatrixType a = MatrixType::Random(rows,cols);
MatrixType a1 = MatrixType::Random(rows,cols); MatrixType a1 = MatrixType::Random(rows,cols);

View File

@ -29,7 +29,7 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
double density = std::max(8./(rows*cols), 0.01); double density = std::max(8./(rows*cols), 0.01);
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector; typedef Matrix<Scalar,Dynamic,1> DenseVector;
Scalar eps = 1e-6; // Scalar eps = 1e-6;
DenseVector vec1 = DenseVector::Random(rows); DenseVector vec1 = DenseVector::Random(rows);
@ -128,7 +128,9 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
LU<DenseMatrix> refLu(refMat2); LU<DenseMatrix> refLu(refMat2);
refLu.solve(b, &refX); refLu.solve(b, &refX);
#if defined(EIGEN_SUPERLU_SUPPORT) || defined(EIGEN_UMFPACK_SUPPORT)
Scalar refDet = refLu.determinant(); Scalar refDet = refLu.determinant();
#endif
x.setZero(); x.setZero();
// // SparseLU<SparseMatrix<Scalar> > (m2).solve(b,&x); // // SparseLU<SparseMatrix<Scalar> > (m2).solve(b,&x);
// // VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: default"); // // VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: default");