Lines 232-313
template<typename VectorsType, typename
Link Here
|
232 |
{ |
232 |
{ |
233 |
return conjugate().setTrans(!m_trans); |
233 |
return conjugate().setTrans(!m_trans); |
234 |
} |
234 |
} |
235 |
|
235 |
|
236 |
/** \brief Inverse of the Householder sequence (equals the adjoint). */ |
236 |
/** \brief Inverse of the Householder sequence (equals the adjoint). */ |
237 |
ConjugateReturnType inverse() const { return adjoint(); } |
237 |
ConjugateReturnType inverse() const { return adjoint(); } |
238 |
|
238 |
|
239 |
/** \internal */ |
239 |
/** \internal */ |
240 |
template<typename DestType> void evalTo(DestType& dst) const |
240 |
template<typename DestType> inline void evalTo(DestType& dst) const |
241 |
{ |
241 |
{ |
|
|
242 |
Matrix<Scalar, DestType::RowsAtCompileTime, 1, |
243 |
AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows()); |
244 |
evalTo(dst, workspace); |
245 |
} |
246 |
|
247 |
/** \internal */ |
248 |
template<typename Dest, typename Workspace> |
249 |
void evalTo(Dest& dst, MatrixBase<Workspace>& workspace) const |
250 |
{ |
251 |
workspace.derived().resize(rows()); |
242 |
Index vecs = m_length; |
252 |
Index vecs = m_length; |
243 |
// FIXME find a way to pass this temporary if the user wants to |
253 |
if( internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value |
244 |
Matrix<Scalar, DestType::RowsAtCompileTime, 1, |
|
|
245 |
AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> temp(rows()); |
246 |
if( internal::is_same<typename internal::remove_all<VectorsType>::type,DestType>::value |
247 |
&& internal::extract_data(dst) == internal::extract_data(m_vectors)) |
254 |
&& internal::extract_data(dst) == internal::extract_data(m_vectors)) |
248 |
{ |
255 |
{ |
249 |
// in-place |
256 |
// in-place |
250 |
dst.diagonal().setOnes(); |
257 |
dst.diagonal().setOnes(); |
251 |
dst.template triangularView<StrictlyUpper>().setZero(); |
258 |
dst.template triangularView<StrictlyUpper>().setZero(); |
252 |
for(Index k = vecs-1; k >= 0; --k) |
259 |
for(Index k = vecs-1; k >= 0; --k) |
253 |
{ |
260 |
{ |
254 |
Index cornerSize = rows() - k - m_shift; |
261 |
Index cornerSize = rows() - k - m_shift; |
255 |
if(m_trans) |
262 |
if(m_trans) |
256 |
dst.bottomRightCorner(cornerSize, cornerSize) |
263 |
dst.bottomRightCorner(cornerSize, cornerSize) |
257 |
.applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0)); |
264 |
.applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); |
258 |
else |
265 |
else |
259 |
dst.bottomRightCorner(cornerSize, cornerSize) |
266 |
dst.bottomRightCorner(cornerSize, cornerSize) |
260 |
.applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0)); |
267 |
.applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); |
261 |
|
268 |
|
262 |
// clear the off diagonal vector |
269 |
// clear the off diagonal vector |
263 |
dst.col(k).tail(rows()-k-1).setZero(); |
270 |
dst.col(k).tail(rows()-k-1).setZero(); |
264 |
} |
271 |
} |
265 |
// clear the remaining columns if needed |
272 |
// clear the remaining columns if needed |
266 |
for(Index k = 0; k<cols()-vecs ; ++k) |
273 |
for(Index k = 0; k<cols()-vecs ; ++k) |
267 |
dst.col(k).tail(rows()-k-1).setZero(); |
274 |
dst.col(k).tail(rows()-k-1).setZero(); |
268 |
} |
275 |
} |
269 |
else |
276 |
else |
270 |
{ |
277 |
{ |
271 |
dst.setIdentity(rows(), rows()); |
278 |
dst.setIdentity(rows(), rows()); |
272 |
for(Index k = vecs-1; k >= 0; --k) |
279 |
for(Index k = vecs-1; k >= 0; --k) |
273 |
{ |
280 |
{ |
274 |
Index cornerSize = rows() - k - m_shift; |
281 |
Index cornerSize = rows() - k - m_shift; |
275 |
if(m_trans) |
282 |
if(m_trans) |
276 |
dst.bottomRightCorner(cornerSize, cornerSize) |
283 |
dst.bottomRightCorner(cornerSize, cornerSize) |
277 |
.applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0)); |
284 |
.applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); |
278 |
else |
285 |
else |
279 |
dst.bottomRightCorner(cornerSize, cornerSize) |
286 |
dst.bottomRightCorner(cornerSize, cornerSize) |
280 |
.applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0)); |
287 |
.applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); |
281 |
} |
288 |
} |
282 |
} |
289 |
} |
283 |
} |
290 |
} |
284 |
|
291 |
|
285 |
/** \internal */ |
292 |
/** \internal */ |
286 |
template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const |
293 |
template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const |
287 |
{ |
294 |
{ |
288 |
Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows()); |
295 |
Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows()); |
|
|
296 |
applyThisOnTheRight(dst, workspace); |
297 |
} |
298 |
|
299 |
/** \internal */ |
300 |
template<typename Dest, typename Workspace> |
301 |
inline void applyThisOnTheRight(Dest& dst, MatrixBase<Workspace>& workspace) const |
302 |
{ |
303 |
workspace.derived().resize(dst.rows()); |
289 |
for(Index k = 0; k < m_length; ++k) |
304 |
for(Index k = 0; k < m_length; ++k) |
290 |
{ |
305 |
{ |
291 |
Index actual_k = m_trans ? m_length-k-1 : k; |
306 |
Index actual_k = m_trans ? m_length-k-1 : k; |
292 |
dst.rightCols(rows()-m_shift-actual_k) |
307 |
dst.rightCols(rows()-m_shift-actual_k) |
293 |
.applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0)); |
308 |
.applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &workspace.coeffRef(0)); |
294 |
} |
309 |
} |
295 |
} |
310 |
} |
296 |
|
311 |
|
297 |
/** \internal */ |
312 |
/** \internal */ |
298 |
template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const |
313 |
template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const |
299 |
{ |
314 |
{ |
300 |
Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols()); |
315 |
Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols()); |
|
|
316 |
applyThisOnTheLeft(dst, workspace); |
317 |
} |
318 |
|
319 |
/** \internal */ |
320 |
template<typename Dest, typename Workspace> |
321 |
inline void applyThisOnTheLeft(Dest& dst, MatrixBase<Workspace>& workspace) const |
322 |
{ |
323 |
workspace.derived().resize(dst.cols()); |
301 |
for(Index k = 0; k < m_length; ++k) |
324 |
for(Index k = 0; k < m_length; ++k) |
302 |
{ |
325 |
{ |
303 |
Index actual_k = m_trans ? k : m_length-k-1; |
326 |
Index actual_k = m_trans ? k : m_length-k-1; |
304 |
dst.bottomRows(rows()-m_shift-actual_k) |
327 |
dst.bottomRows(rows()-m_shift-actual_k) |
305 |
.applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0)); |
328 |
.applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &workspace.coeffRef(0)); |
306 |
} |
329 |
} |
307 |
} |
330 |
} |
308 |
|
331 |
|
309 |
/** \brief Computes the product of a Householder sequence with a matrix. |
332 |
/** \brief Computes the product of a Householder sequence with a matrix. |
310 |
* \param[in] other %Matrix being multiplied. |
333 |
* \param[in] other %Matrix being multiplied. |
311 |
* \returns Expression object representing the product. |
334 |
* \returns Expression object representing the product. |
312 |
* |
335 |
* |
313 |
* This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this |
336 |
* This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this |