用拉帕克dgeqrf_求解线性系统
Solving a linear system with Lapack's dgeqrf_
我正试图用C++中的QR因子分解对矩阵进行因子分解，使用Lapack函数来求解线性方程组(Ax=b)
据我所知，dgeqrf计算QR因子分解并覆盖输入矩阵。输出中显然包含L(上三角)的值，但我如何获得Q？
我尝试了dormqr，据说它根据dgeqrf的输出计算Q，但结果与上一次调用中的矩阵相同。
这是我的完整代码：
boost::numeric::ublas::matrix<double> in_A(4, 3);
in_A(0, 0) = 1.0;
in_A(0, 1) = 2.0;
in_A(0, 2) = 3.0;
in_A(1, 1) = -3.0;
in_A(1, 2) = 2.0;
in_A(1, 3) = 1.0;
in_A(2, 1) = 2.0;
in_A(2, 2) = 0.0;
in_A(2, 3) = -1.0;
in_A(3, 1) = 3.0;
in_A(3, 2) = -1.0;
in_A(3, 3) = 2.0;
boost::numeric::ublas::vector<double> in_b(4);
in_b(0) = 2;
in_b(1) = 4;
in_b(2) = 6;
in_b(3) = 8;
int rows = in_A.size1();
int cols = in_A.size2();
double *A = (double *)malloc(rows*cols*sizeof(double));
double *b = (double *)malloc(in_b.size()*sizeof(double));
//Lapack has column-major order
for(size_t col=0; col<in_A.size2(); ++col)
{
for(size_t row = 0; row<in_A.size1(); ++row)
{
int D1_idx = col*in_A.size1() + row;
A[D1_idx] = in_A(row, col);
}
b[col] = in_b(col);
}
integer m = rows;
integer n = cols;
integer info = 0;
integer k = n; /* k = min(m,n); */
integer lda = m; /* lda = max(m,1); */
integer lwork = n; /* lwork = max(n,1); */
int max = lwork; /* max = max(lwork,1); */
double *work;
double *tau;
char *side = "L";
char *TR = "T";
integer one = 1;
int i;
double *vec;
work = (double *) malloc( max * sizeof( double ) );
tau = (double *) malloc( k * sizeof( double ) );
vec = (double *) malloc( m * sizeof( double ) );
memset(work, 0, max * sizeof(double));
memset(tau, 0, k * sizeof(double));
std::cout << std::endl;
for(size_t row = 0; row < rows; ++row)
{
for(size_t col = 0; col < cols; ++col)
{
size_t idx = col*rows + row;
std::cout << A[idx] << " ";
}
std::cout << std::endl;
}
dgeqrf_(&m, &n, A, &lda, tau, work, &lwork, &info);
//printf("tau[0] = %f tau[1] = %fn",tau[0],tau[1]);
std::cout << std::endl;
for(size_t row = 0; row < rows; ++row)
{
for(size_t col = 0; col < cols; ++col)
{
size_t idx = col*rows + row;
std::cout << A[idx] << " ";
}
std::cout << std::endl;
}
memset(vec, 0, m * sizeof(double));
vec[2] = 1.0;
dormqr_(side, TR, &m, &one, &k, A, &lda, tau, vec, &lda, work, &lwork, &info);
free(vec);
free(tau);
free(work);
我的代码出了什么问题？
如何分解矩阵并求解相应的线性方程组？
根据中的文档
(http://www.netlib.org/lapack/explore-html/da/d82/dormqr_8f.html)
你在vec中计算乘积Q^T*e3，其中e3是第三个规范基向量(0,0,1,0，…，0)。如果你想计算Q，那么vec应该包含一个矩阵大小的数组，用单位矩阵填充，TRANS应该是"N"。
dormqr (SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SIDE="L"，对于Q为左的正常QR分解，
TRANS="N"返回QC代替C
A在存储器中具有布局LDA x K，其中上部M x K块被使用并编码K个反射器
τ包含K个反射器的因子
C在存储器中具有布局LDCxM，其中上部MxN块将用于保存结果QC
为了让C在返回时保持Q，C必须是初始化为恒等式的平方MxM矩阵，即对角线项都为1。
您可以考虑使用为ublas提供的lapack数字绑定，如
(http://boost.2283326.n4.nabble.com/How-to-use-the-qr-decomposition-correctly-td2710159.html)
然而，这个项目现在可能已经失效或停止。
让我们从第一原则开始：目标是求解Ax=b，或者至少最小化|Ax-b|+|x|。为了使其一致，需要colsA=rowsx和rowsA=rowsb。
现在，为了使所讨论的代码工作，A必须是正方形或高矩形矩阵colsA<=rowsA，因此系统是超定的。
计算步骤
求解
Q*R=A：(http://www.netlib.no/netlib/lapack/double/dgeqrf.f)DGEQRF(行A、列A、A、行A、TAU、WORK、LWORK、INFO)
与
R*x=QT*b中的QT相乘得到QT*b(http://www.netlib.no/netlib/lapack/double/dormqr.f)DOMQR('L'，'T'，行A，1，列A，A，行A、TAU、b，行A和WORK，LWORK，INFO)
使用
A的右上部分进行反向替换(http://www.netlib.no/netlib/lapack/double/dtrtrs.f)DTRTRS('U'，'N'，'N'，列A，1，A，行A，b，行A和信息)
现在，
b的第一个colsA条目包含解向量x。索引colsA+1及其后的剩余项的欧氏范数是解的误差|A*x-b|。
备注：对于纯求解过程，没有理由显式计算"Q"或调用通用矩阵乘法DGEMM。这些应该保留给实验，以检查A-QR是否足够接近零。
备注：通过执行LWORK=-1的试运行来探索WORK阵列的最佳分配。
然而，为了总结一些有效的代码，ublas和lapack之间的连接似乎不是最佳的
#include "boost/numeric/ublas/matrix.hpp"
#include "boost/numeric/ublas/vector.hpp"
typedef boost::numeric::ublas::matrix<double> bmatrix;
typedef boost::numeric::ublas::vector<double> bvector;
namespace lapack {
extern "C" {
void dgeqrf_(int* M, int* N,
double* A, int* LDA, double* TAU,
double* WORK, int* LWORK, int* INFO );
void dormqr_(char* SIDE, char* TRANS,
int* M, int* N, int* K,
double* A, int* LDA, double* TAU,
double* C, int* LDC,
double* WORK, int* LWORK, int* INFO );
void dtrtrs_(char* UPLO, char* TRANS, char* DIAG,
int* N, int* NRHS,
double* A, int* LDA,
double* B, int* LDB,
int* INFO );
}
int geqrf(int m, int n,
double* A, int lda, double *tau) {
int info=0;
int lwork=-1;
double iwork;
dgeqrf_(&m, &n, A, &lda, tau,
&iwork, &lwork, &info);
lwork = (int)iwork;
double* work = new double[lwork];
dgeqrf_(&m, &n, A, &lda, tau,
work, &lwork, &info);
delete[] work;
return info;
}
int ormqr(char side, char trans, int m, int n, int k,
double *A, int lda, double *tau, double* C, int ldc) {
int info=0;
int lwork=-1;
double iwork;
dormqr_(&side, &trans, &m, &n, &k,
A, &lda, tau, C, &ldc, &iwork, &lwork, &info);
lwork = (int)iwork;
double* work = new double[lwork];
dormqr_(&side, &trans, &m, &n, &k,
A, &lda, tau, C, &ldc, work, &lwork, &info);
delete[] work;
return info;
}
int trtrs(char uplo, char trans, char diag,
int n, int nrhs,
double* A, int lda, double* B, int ldb
) {
int info = 0;
dtrtrs_(&uplo, &trans, &diag, &n, &nrhs,
A, &lda, B, &ldb, &info);
return info;
}
}
static void PrintMatrix(double A[], size_t rows, size_t cols) {
std::cout << std::endl;
for(size_t row = 0; row < rows; ++row)
{
for(size_t col = 0; col < cols; ++col)
{
// Lapack uses column major format
size_t idx = col*rows + row;
std::cout << A[idx] << " ";
}
std::cout << std::endl;
}
}
static int SolveQR(
const bmatrix &in_A, // IN
const bvector &in_b, // IN
bvector &out_x // OUT
) {
size_t rows = in_A.size1();
size_t cols = in_A.size2();
double *A = new double[rows*cols];
double *b = new double[in_b.size()];
//Lapack has column-major order
for(size_t col=0, D1_idx=0; col<cols; ++col)
{
for(size_t row = 0; row<rows; ++row)
{
// Lapack uses column major format
A[D1_idx++] = in_A(row, col);
}
b[col] = in_b(col);
}
for(size_t row = 0; row<rows; ++row)
{
b[row] = in_b(row);
}
// DGEQRF for Q*R=A, i.e., A and tau hold R and Householder reflectors
double* tau = new double[cols];
PrintMatrix(A, rows, cols);
lapack::geqrf(rows, cols, A, rows, tau);
PrintMatrix(A, rows, cols);
// DORMQR: to compute b := Q^T*b
lapack::ormqr('L', 'T', rows, 1, cols, A, rows, tau, b, rows);
PrintMatrix(b, rows, 1);
// DTRTRS: solve Rx=b by back substitution
lapack::trtrs('U', 'N', 'N', cols, 1, A, rows, b, rows);
for(size_t col=0; col<cols; col++) {
out_x(col)=b[col];
}
PrintMatrix(b,cols,1);
delete[] A;
delete[] b;
delete[] tau;
return 0;
}
int main() {
bmatrix in_A(4, 3);
in_A(0, 0) = 1.0; in_A(0, 1) = 2.0; in_A(0, 2) = 3.0;
in_A(1, 0) = -3.0; in_A(1, 1) = 2.0; in_A(1, 2) = 1.0;
in_A(2, 0) = 2.0; in_A(2, 1) = 0.0; in_A(2, 2) = -1.0;
in_A(3, 0) = 3.0; in_A(3, 1) = -1.0; in_A(3, 2) = 2.0;
bvector in_b(4);
in_b(0) = 2;
in_b(1) = 4;
in_b(2) = 6;
in_b(3) = 8;
bvector out_x(3);
SolveQR( in_A, in_b, out_x);
return 0;
}
虽然这是一个老问题，但如果您正在寻找一种使用QR和LAPACK使用dgels解决LLS的方法，它与上面的答案相同。
相关文章：
- C++ 本征线性系统求解，数值问题?
- 使用特征 3 线性系统求解器的错误结果
- 使用特征/英特尔 MKL 求解稀疏线性系统
- 犰狳库中线性系统的近似解
- 用SimplicialCholesky Eigen求解大型稀疏线性系统
- 特征 3 断言在求解线性系统时失败 - 据我所知，这是由于特征中的无效索引
- 在C++-GPU中解决稀疏线性系统的最佳方法可能吗
- 如何使用Z_2中的系数求解稀疏线性系统
- 以线性最小二乘方式求解系统 Ax=b，具有复元素和下三角形平方 A 矩阵
- 用拉帕克dgeqrf_求解线性系统
- 定维(N=9）对称半正定稠密线性系统的快速解
- 当仅恒定项之一变化时，线性系统AX = B的有效解决方案
- C++ 用于求解复杂线性系统 Ax=b 的库
- 使用指针和数组求解线性系统
- 无法访问特征线性系统解的元素
- 如何求解复杂的线性系统
- 改进稀疏线性系统的解
- c++中有没有免费的迭代线性系统求解器，可以让我输入任意的初始猜测
- 如何求解非常大的系统的稀疏线性系统
- 在CUDA中求解稀疏定正线性系统
