使用Lapack的dgeqrf_求解线性系统 [英] Solving a linear system with Lapack's dgeqrf_
问题描述
我试图使用C ++中的 QR因式分解因子分解矩阵，使用Lapack的函数为了解决一个线性方程组（Ax = b）
根据我的理解， dgeqrf 计算QR因式分解并覆盖输入矩阵。输出清楚地包含L（上三角形）的值，但是如何获得Q？
我尝试了 dormqr
这是我完整的代码：
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有列主命令
for（size_t col = 0; col
for size_t row = 0; row< in_A.size1（）; ++ row）
{
int D1_idx = col * in_A.size1（）
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] =％f \\\
，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，& ;
free（vec）;
free（tau）;
免费（工作）;
我的代码有什么问题？
如何根据矩阵进行因式分解，并求解相应的线性方程组？
b
$ b
（ http：// www.netlib.org/lapack/explore-html/da/d82/dormqr_8f.html ）
您在计算产品Q ^ T * e3，其中e3是第三规范基本向量（0,0,1,0,0，...，0）。如果你想计算Q，那么vec应该包含一个用单位矩阵填充的矩阵大小的数组，TRANS应该是N。
dormqr（SIDE，TRANS，M，N，K，A，LDA，TAU，C，LDC，WORK，LWORK，INFO）
-
对于Q左侧的正常QR分解，SIDE =L，
/ li>
-
TRANS =N代替C
返回QC -
A具有布局LDA x K
> -
C在内存中具有布局LDC x M，其中上部M x N块将用于保存结果QC
-
对于C在返回时保持Q，C必须是初始化为identity的正方形M×M矩阵，即对角项为1。
您可以考虑使用为ublas提供的lapack数字绑定，例如
（ http：// boost .2283326.n4.nabble.com / How-to-use-the-qr-decomposition-correctly-td2710159.html ）
但是，此项目
是解决A * x = b，或至少最小化| A * xb | + | x |。为了保持一致，需要 colsA = rowsx 和 rowsA = rowsb 。 现在为讨论的代码工作 A 必须是正方形或高矩形矩阵，<$ c $
$ b
-
解决
Q * R = A：
（ http://www.netlib.no/netlib/lapack/double/dgeqrf.f ）
DGEQRF（rowsA，colsA，A，rowsA，TAU，WORK，LWORK，INFO）
-
乘以
QT以获得QT * b，如R * x = QT * b
（ http：//www.netlib .no / netlib / lapack / double / dormqr.f ）
DORMQR（'L'，'T'，rowsA，1，colsA，
-
使用
A 的右上角部分替换/ code>
（ http：//www.netlib。 no / netlib / lapack / double / dtrtrs.f ）
DTRTRS（'U'，'N'，'N'，colsA，1，A ，rowsA，b，rowsA，INFO）
-
现在第
colsAc> b 包含解决方案向量x。在索引colsA + 1和其后的剩余条目的欧几里得范数是误差| A * x-b |
备注：对于纯解决方案流程，没有理由明确地计算Q或调用通用矩阵乘法DGEMM。这些应保留用于实验，以检查 A-QR 是否足够接近零。
：通过执行LWORK = -1的干运行来探索WORK数组的最佳分配。
然而，总而言之，一些代码可以工作，ublas和lapack之间的连接似乎不是最佳的
#includeboost / numeric / matrix.hpp
#includeboost / numeric / ublas / vector.hpp
typedef boost :: numeric :: ublas :: matrix< double>巴马克斯
typedef boost :: numeric :: ublas :: vector< double> bvector;
命名空间lapack {
externC{
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 * 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使用列主格式
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_Asize2（）;
double * A = new double [rows * cols];
double * b = new double [in_b.size（）];
// Lapack有列主命令
for（size_t col = 0，D1_idx = 0; col< cols; ++ col）
{
for Lapack使用列主格式
A [D1_idx ++] = in_A（row，col）;
}
b [col] = in_b（col）;
}
for（size_t row = 0; row< rows; ++ row）
{
b [row] = in_b
}
//对于Q * R = A的DGEQRF，即A和tau保持R和Householder反射器
double * tau = new double [cols];
PrintMatrix（A，rows，cols）;
lapack :: geqrf（rows，cols，A，rows，tau）;
PrintMatrix（A，rows，cols）;
// DORMQR：计算b：= Q ^ T * b
lapack :: ormqr（'L'，'T'，rows，1，cols，A ，rows，tau，b，rows）;
PrintMatrix（b，rows，1）;
// DTRTRS：通过返回替换解决Rx = b
lapack :: trtrs（'U'，'N'，'N'，cols，1， row，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;
}
I am trying to factorize a matrix with the QR factorization in C++, using Lapack's functions in order to solve a system of linear equations (Ax=b)
As far as I understood, dgeqrf computes the QR factorization and overwrites the input matrix. The output clearly contains values for L (upper triangular), but how do I obtain Q?
I tried dormqr, which is said to calculate Q from dgeqrf's output, but the result is the same matrix as in the previous call.
Here's my complete code:
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] = %f\n",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);
What's wrong with my code?
How can I factorize a matrix and solve a corresponding system of linear equations?
According to the documentation in
(http://www.netlib.org/lapack/explore-html/da/d82/dormqr_8f.html)
you are computing in vec the product Q^T*e3, where e3 is the third canonical basis vector (0,0,1,0,0,...,0). If you want to compute Q, then vec should contain a matrix sized array filled with the unit matrix, and TRANS should be "N".
dormqr (SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SIDE = "L" for the normal QR decomposition with Q left,
TRANS = "N" to return QC in the place of C
A has layout LDA x K in memory, of which the upper M x K block is used and encodes K reflectors
tau contains the factors for the K reflectors
C has layout LDC x M in memory, of which the upper M x N block will be used to hold the result QC
For C to hold Q on return, C must be a square M x M matrix initialized as identity, i.e., with diagonal entries all 1.
You might consider to use the lapack numeric bindings provided for ublas, as in
(http://boost.2283326.n4.nabble.com/How-to-use-the-qr-decomposition-correctly-td2710159.html)
However, this project may be defunct or resting by now.
Lets start again from first principles: The aim is to solve A*x=b, or at least to minimize |A*x-b|+|x|. For that to be consistent one needs colsA=rowsx and rowsA=rowsb.
Now for the discussed code to work A has to be square or a tall rectangular matrix, colsA<=rowsA, so that the system is overdetermined.
Computation steps
Solve
Q*R=A: (http://www.netlib.no/netlib/lapack/double/dgeqrf.f)DGEQRF( rowsA, colsA, A, rowsA, TAU, WORK, LWORK, INFO )
Multiply by
QTto getQT*bas inR*x=QT*b(http://www.netlib.no/netlib/lapack/double/dormqr.f)DORMQR( 'L', 'T', rowsA, 1, colsA, A, rowsA, TAU, b, rowsA, WORK, LWORK, INFO )
Use back-substitution using the upper right part of
A(http://www.netlib.no/netlib/lapack/double/dtrtrs.f)DTRTRS( 'U', 'N', 'N', colsA, 1, A, rowsA, b, rowsA, INFO )
Now the first
colsAentries ofbcontain the solution vectorx. The euclidean norm of the remaining entries at index colsA+1 and thereafter is the error |A*x-b| of the solution.
Remark: For the pure solution process there is no reason to compute 'Q' explicitly or to invoke the generic matrix multiplication DGEMM. These should be reserved for experiments to check if A-QR is sufficiently close to zero.
Remark: Explore the optimal allocation of the WORK array by performing a dry run with LWORK=-1.
To conclude some code that works, however, the connection between ublas and lapack seems suboptimal
#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;
}
这篇关于使用Lapack的dgeqrf_求解线性系统的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
