Lines 251-306
void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename
Link Here
|
251 |
} |
251 |
} |
252 |
|
252 |
|
253 |
/** \internal */ |
253 |
/** \internal */ |
254 |
template<typename MatrixQR, typename HCoeffs> |
254 |
template<typename MatrixQR, typename HCoeffs, |
255 |
void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs, |
255 |
typename MatrixQRScalar = typename MatrixQR::Scalar, |
256 |
typename MatrixQR::Index maxBlockSize=32, |
256 |
bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)> |
257 |
typename MatrixQR::Scalar* tempData = 0) |
257 |
struct householder_qr_inplace_blocked |
258 |
{ |
258 |
{ |
259 |
typedef typename MatrixQR::Index Index; |
259 |
// This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h |
260 |
typedef typename MatrixQR::Scalar Scalar; |
260 |
static void run(MatrixQR& mat, HCoeffs& hCoeffs, |
261 |
typedef Block<MatrixQR,Dynamic,Dynamic> BlockType; |
261 |
typename MatrixQR::Index maxBlockSize=32, |
262 |
|
262 |
typename MatrixQR::Scalar* tempData = 0) |
263 |
Index rows = mat.rows(); |
|
|
264 |
Index cols = mat.cols(); |
265 |
Index size = (std::min)(rows, cols); |
266 |
|
267 |
typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType; |
268 |
TempType tempVector; |
269 |
if(tempData==0) |
270 |
{ |
263 |
{ |
271 |
tempVector.resize(cols); |
264 |
typedef typename MatrixQR::Index Index; |
272 |
tempData = tempVector.data(); |
265 |
typedef typename MatrixQR::Scalar Scalar; |
273 |
} |
266 |
typedef Block<MatrixQR,Dynamic,Dynamic> BlockType; |
274 |
|
|
|
275 |
Index blockSize = (std::min)(maxBlockSize,size); |
276 |
|
277 |
Index k = 0; |
278 |
for (k = 0; k < size; k += blockSize) |
279 |
{ |
280 |
Index bs = (std::min)(size-k,blockSize); // actual size of the block |
281 |
Index tcols = cols - k - bs; // trailing columns |
282 |
Index brows = rows-k; // rows of the block |
283 |
|
267 |
|
284 |
// partition the matrix: |
268 |
Index rows = mat.rows(); |
285 |
// A00 | A01 | A02 |
269 |
Index cols = mat.cols(); |
286 |
// mat = A10 | A11 | A12 |
270 |
Index size = (std::min)(rows, cols); |
287 |
// A20 | A21 | A22 |
|
|
288 |
// and performs the qr dec of [A11^T A12^T]^T |
289 |
// and update [A21^T A22^T]^T using level 3 operations. |
290 |
// Finally, the algorithm continue on A22 |
291 |
|
271 |
|
292 |
BlockType A11_21 = mat.block(k,k,brows,bs); |
272 |
typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType; |
293 |
Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs); |
273 |
TempType tempVector; |
|
|
274 |
if(tempData==0) |
275 |
{ |
276 |
tempVector.resize(cols); |
277 |
tempData = tempVector.data(); |
278 |
} |
294 |
|
279 |
|
295 |
householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData); |
280 |
Index blockSize = (std::min)(maxBlockSize,size); |
296 |
|
281 |
|
297 |
if(tcols) |
282 |
Index k = 0; |
|
|
283 |
for (k = 0; k < size; k += blockSize) |
298 |
{ |
284 |
{ |
299 |
BlockType A21_22 = mat.block(k,k+bs,brows,tcols); |
285 |
Index bs = (std::min)(size-k,blockSize); // actual size of the block |
300 |
apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint()); |
286 |
Index tcols = cols - k - bs; // trailing columns |
|
|
287 |
Index brows = rows-k; // rows of the block |
288 |
|
289 |
// partition the matrix: |
290 |
// A00 | A01 | A02 |
291 |
// mat = A10 | A11 | A12 |
292 |
// A20 | A21 | A22 |
293 |
// and performs the qr dec of [A11^T A12^T]^T |
294 |
// and update [A21^T A22^T]^T using level 3 operations. |
295 |
// Finally, the algorithm continue on A22 |
296 |
|
297 |
BlockType A11_21 = mat.block(k,k,brows,bs); |
298 |
Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs); |
299 |
|
300 |
householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData); |
301 |
|
302 |
if(tcols) |
303 |
{ |
304 |
BlockType A21_22 = mat.block(k,k+bs,brows,tcols); |
305 |
apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint()); |
306 |
} |
301 |
} |
307 |
} |
302 |
} |
308 |
} |
303 |
} |
309 |
}; |
304 |
|
310 |
|
305 |
template<typename _MatrixType, typename Rhs> |
311 |
template<typename _MatrixType, typename Rhs> |
306 |
struct solve_retval<HouseholderQR<_MatrixType>, Rhs> |
312 |
struct solve_retval<HouseholderQR<_MatrixType>, Rhs> |
Lines 352-358
HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
Link Here
|
352 |
|
358 |
|
353 |
m_temp.resize(cols); |
359 |
m_temp.resize(cols); |
354 |
|
360 |
|
355 |
internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data()); |
361 |
internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data()); |
356 |
|
362 |
|
357 |
m_isInitialized = true; |
363 |
m_isInitialized = true; |
358 |
return *this; |
364 |
return *this; |