diff options
Diffstat (limited to 'Eigen/src/Core/SolveTriangular.h')
| -rw-r--r-- | Eigen/src/Core/SolveTriangular.h | 74 |
1 files changed, 37 insertions, 37 deletions
diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index b58dab01d..1c586b865 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -30,9 +30,9 @@ template<typename XprType, unsigned int Mode> struct ei_is_part<Part<XprType,Mod template<typename Lhs, typename Rhs, int TriangularPart = (int(Lhs::Flags) & LowerTriangularBit) - ? Lower + ? LowerTriangular : (int(Lhs::Flags) & UpperTriangularBit) - ? Upper + ? UpperTriangular : -1, int StorageOrder = ei_is_part<Lhs>::value ? -1 // this is to solve ambiguous specializations : int(Lhs::Flags) & (RowMajorBit|SparseBit) @@ -56,14 +56,14 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) { - const bool IsLower = (UpLo==Lower); + const bool IsLowerTriangular = (UpLo==LowerTriangular); const int size = lhs.cols(); /* We perform the inverse product per block of 4 rows such that we perfectly match * our optimized matrix * vector product. blockyStart represents the number of rows * we have process first using the non-block version. */ int blockyStart = (std::max(size-5,0)/4)*4; - if (IsLower) + if (IsLowerTriangular) blockyStart = size - blockyStart; else blockyStart -= 1; @@ -72,15 +72,15 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> // process first rows using the non block version if(!(Lhs::Flags & UnitDiagBit)) { - if (IsLower) + if (IsLowerTriangular) other.coeffRef(0,c) = other.coeff(0,c)/lhs.coeff(0, 0); else other.coeffRef(size-1,c) = other.coeff(size-1, c)/lhs.coeff(size-1, size-1); } - for(int i=(IsLower ? 1 : size-2); IsLower ? i<blockyStart : i>blockyStart; i += (IsLower ? 1 : -1) ) + for(int i=(IsLowerTriangular ? 1 : size-2); IsLowerTriangular ? i<blockyStart : i>blockyStart; i += (IsLowerTriangular ? 1 : -1) ) { Scalar tmp = other.coeff(i,c) - - (IsLower ? ((lhs.row(i).start(i)) * other.col(c).start(i)).coeff(0,0) + - (IsLowerTriangular ? ((lhs.row(i).start(i)) * other.col(c).start(i)).coeff(0,0) : ((lhs.row(i).end(size-i-1)) * other.col(c).end(size-i-1)).coeff(0,0)); if (Lhs::Flags & UnitDiagBit) other.coeffRef(i,c) = tmp; @@ -89,15 +89,15 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> } // now let's process the remaining rows 4 at once - for(int i=blockyStart; IsLower ? i<size : i>0; ) + for(int i=blockyStart; IsLowerTriangular ? i<size : i>0; ) { int startBlock = i; - int endBlock = startBlock + (IsLower ? 4 : -4); + int endBlock = startBlock + (IsLowerTriangular ? 4 : -4); /* Process the i cols times 4 rows block, and keep the result in a temporary vector */ // FIXME use fixed size block but take care to small fixed size matrices... Matrix<Scalar,Dynamic,1> btmp(4); - if (IsLower) + if (IsLowerTriangular) btmp = lhs.block(startBlock,0,4,i) * other.col(c).start(i); else btmp = lhs.block(i-3,i+1,4,size-1-i) * other.col(c).end(size-1-i); @@ -106,21 +106,21 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> * btmp stores the diagonal coefficients used to update the remaining part of the result. */ { - Scalar tmp = other.coeff(startBlock,c)-btmp.coeff(IsLower?0:3); + Scalar tmp = other.coeff(startBlock,c)-btmp.coeff(IsLowerTriangular?0:3); if (Lhs::Flags & UnitDiagBit) other.coeffRef(i,c) = tmp; else other.coeffRef(i,c) = tmp/lhs.coeff(i,i); } - i += IsLower ? 1 : -1; - for (;IsLower ? i<endBlock : i>endBlock; i += IsLower ? 1 : -1) + i += IsLowerTriangular ? 1 : -1; + for (;IsLowerTriangular ? i<endBlock : i>endBlock; i += IsLowerTriangular ? 1 : -1) { - int remainingSize = IsLower ? i-startBlock : startBlock-i; + int remainingSize = IsLowerTriangular ? i-startBlock : startBlock-i; Scalar tmp = other.coeff(i,c) - - btmp.coeff(IsLower ? remainingSize : 3-remainingSize) - - ( lhs.row(i).segment(IsLower ? startBlock : i+1, remainingSize) - * other.col(c).segment(IsLower ? startBlock : i+1, remainingSize)).coeff(0,0); + - btmp.coeff(IsLowerTriangular ? remainingSize : 3-remainingSize) + - ( lhs.row(i).segment(IsLowerTriangular ? startBlock : i+1, remainingSize) + * other.col(c).segment(IsLowerTriangular ? startBlock : i+1, remainingSize)).coeff(0,0); if (Lhs::Flags & UnitDiagBit) other.coeffRef(i,c) = tmp; @@ -133,10 +133,10 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> }; // Implements the following configurations: -// - inv(Lower, ColMajor) * Column vector -// - inv(Lower,UnitDiag,ColMajor) * Column vector -// - inv(Upper, ColMajor) * Column vector -// - inv(Upper,UnitDiag,ColMajor) * Column vector +// - inv(LowerTriangular, ColMajor) * Column vector +// - inv(LowerTriangular,UnitDiag,ColMajor) * Column vector +// - inv(UpperTriangular, ColMajor) * Column vector +// - inv(UpperTriangular,UnitDiag,ColMajor) * Column vector template<typename Lhs, typename Rhs, int UpLo> struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> { @@ -146,7 +146,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> static void run(const Lhs& lhs, Rhs& other) { - static const bool IsLower = (UpLo==Lower); + static const bool IsLowerTriangular = (UpLo==LowerTriangular); const int size = lhs.cols(); for(int c=0 ; c<other.cols() ; ++c) { @@ -155,27 +155,27 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> * we can process using the block version. */ int blockyEnd = (std::max(size-5,0)/4)*4; - if (!IsLower) + if (!IsLowerTriangular) blockyEnd = size-1 - blockyEnd; - for(int i=IsLower ? 0 : size-1; IsLower ? i<blockyEnd : i>blockyEnd;) + for(int i=IsLowerTriangular ? 0 : size-1; IsLowerTriangular ? i<blockyEnd : i>blockyEnd;) { /* Let's process the 4x4 sub-matrix as usual. * btmp stores the diagonal coefficients used to update the remaining part of the result. */ int startBlock = i; - int endBlock = startBlock + (IsLower ? 4 : -4); + int endBlock = startBlock + (IsLowerTriangular ? 4 : -4); Matrix<Scalar,4,1> btmp; - for (;IsLower ? i<endBlock : i>endBlock; - i += IsLower ? 1 : -1) + for (;IsLowerTriangular ? i<endBlock : i>endBlock; + i += IsLowerTriangular ? 1 : -1) { if(!(Lhs::Flags & UnitDiagBit)) other.coeffRef(i,c) /= lhs.coeff(i,i); - int remainingSize = IsLower ? endBlock-i-1 : i-endBlock-1; + int remainingSize = IsLowerTriangular ? endBlock-i-1 : i-endBlock-1; if (remainingSize>0) - other.col(c).segment((IsLower ? i : endBlock) + 1, remainingSize) -= + other.col(c).segment((IsLowerTriangular ? i : endBlock) + 1, remainingSize) -= other.coeffRef(i,c) - * Block<Lhs,Dynamic,1>(lhs, (IsLower ? i : endBlock) + 1, i, remainingSize, 1); - btmp.coeffRef(IsLower ? i-startBlock : remainingSize) = -other.coeffRef(i,c); + * Block<Lhs,Dynamic,1>(lhs, (IsLowerTriangular ? i : endBlock) + 1, i, remainingSize, 1); + btmp.coeffRef(IsLowerTriangular ? i-startBlock : remainingSize) = -other.coeffRef(i,c); } /* Now we can efficiently update the remaining part of the result as a matrix * vector product. @@ -187,11 +187,11 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> // FIXME this is cool but what about conjugate/adjoint expressions ? do we want to evaluate them ? // this is a more general problem though. ei_cache_friendly_product_colmajor_times_vector( - IsLower ? size-endBlock : endBlock+1, - &(lhs.const_cast_derived().coeffRef(IsLower ? endBlock : 0, IsLower ? startBlock : endBlock+1)), + IsLowerTriangular ? size-endBlock : endBlock+1, + &(lhs.const_cast_derived().coeffRef(IsLowerTriangular ? endBlock : 0, IsLowerTriangular ? startBlock : endBlock+1)), lhs.stride(), - btmp, &(other.coeffRef(IsLower ? endBlock : 0, c))); -// if (IsLower) + btmp, &(other.coeffRef(IsLowerTriangular ? endBlock : 0, c))); +// if (IsLowerTriangular) // other.col(c).end(size-endBlock) += (lhs.block(endBlock, startBlock, size-endBlock, endBlock-startBlock) // * other.col(c).block(startBlock,endBlock-startBlock)).lazy(); // else @@ -201,7 +201,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> /* Now we have to process the remaining part as usual */ int i; - for(i=blockyEnd; IsLower ? i<size-1 : i>0; i += (IsLower ? 1 : -1) ) + for(i=blockyEnd; IsLowerTriangular ? i<size-1 : i>0; i += (IsLowerTriangular ? 1 : -1) ) { if(!(Lhs::Flags & UnitDiagBit)) other.coeffRef(i,c) /= lhs.coeff(i,i); @@ -209,7 +209,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> /* NOTE we cannot use lhs.col(i).end(size-i-1) because Part::coeffRef gets called by .col() to * get the address of the start of the row */ - if(IsLower) + if(IsLowerTriangular) other.col(c).end(size-i-1) -= other.coeffRef(i,c) * Block<Lhs,Dynamic,1>(lhs, i+1,i, size-i-1,1); else other.col(c).start(i) -= other.coeffRef(i,c) * Block<Lhs,Dynamic,1>(lhs, 0,i, i, 1); |