DGEQRF and SGEQRF from LAPACK 以打包格式返回 QR 分解的 Q 部分。解压它似乎需要 O(k^3)
步骤(k 个低等级产品),而且似乎不是很简单。另外,我不清楚进行 k
顺序乘法的数值稳定性。
LAPACK是否包含用于解包Q的子程序,如果没有,我该怎么做?
最佳答案
是的,LAPACK确实提供了一个从基本反射器(即DGEQRF返回的数据部分)检索Q的例程,它被称为DORGQR 。从描述来看:
* DORGQR generates an M-by-N real matrix Q with orthonormal columns,
* which is defined as the first N columns of a product of K elementary
* reflectors of order M
*
* Q = H(1) H(2) . . . H(k)
* as returned by DGEQRF.
使用 C
包装器 LAPACKE 从 A
完整计算 Q
和 R
可能如下所示(Fortran 改编应该是直接的):
void qr( double* const _Q, double* const _R, double* const _A, const size_t _m, const size_t _n) {
// Maximal rank is used by Lapacke
const size_t rank = std::min(_m, _n);
// Tmp Array for Lapacke
const std::unique_ptr<double[]> tau(new double[rank]);
// Calculate QR factorisations
LAPACKE_dgeqrf(LAPACK_ROW_MAJOR, (int) _m, (int) _n, _A, (int) _n, tau.get());
// Copy the upper triangular Matrix R (rank x _n) into position
for(size_t row =0; row < rank; ++row) {
memset(_R+row*_n, 0, row*sizeof(double)); // Set starting zeros
memcpy(_R+row*_n+row, _A+row*_n+row, (_n-row)*sizeof(double)); // Copy upper triangular part from Lapack result.
}
// Create orthogonal matrix Q (in tmpA)
LAPACKE_dorgqr(LAPACK_ROW_MAJOR, (int) _m, (int) rank, (int) rank, _A, (int) _n, tau.get());
//Copy Q (_m x rank) into position
if(_m == _n) {
memcpy(_Q, _A, sizeof(double)*(_m*_n));
} else {
for(size_t row =0; row < _m; ++row) {
memcpy(_Q+row*rank, _A+row*_n, sizeof(double)*(rank));
}
}
}
这是我自己的一段代码,我删除了所有检查以提高可读性。为了高效使用,您需要检查输入是否有效,并注意 LAPACK 调用的返回值。请注意,输入 A
已被销毁。
关于matrix - 如何从 QR 分解输出中获取 Q?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/29395461/