我正试图对一个矩阵进行因式分解,并将其与 QR因子化 在C++语言中,使用Lapack函数来求解线性方程组(Ax=b)
据我了解。dgeqrf 计算QR分解,并覆盖输入矩阵。输出的结果显然包含了L(上三角)的值,但我如何获得Q?
Ax-b dormqr,据说它可以从 dgeqrf
的输出,但结果和之前调用的矩阵是一样的。
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);
这是我的完整代码。
我的代码有什么问题?
我如何对一个矩阵进行分解,并解出相应的线性方程组?我试图用C++中的QR分解法对一个矩阵进行因子化,使用Lapack的函数来求解一个线性方程组(Ax=b)据我了解,dgeqrf计算QR ...
根据文档中的 (http:/www.netlib.orglapackexplore-htmldad82dormqr_8f.html
)
dormqr (SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
你在vec中计算的是乘积Q^T*e3,其中e3是第三规范基向量(0,0,1,0,0,...,0)。如果你想计算Q,那么vec应该包含一个矩阵大小的数组,里面填充的是单位矩阵,trans应该是 "N"。
A*x-b
TRANS = "N",以返回QC代替C。
A在存储器中具有布局LDA×K,其中上层M×K块被使用,并编码K个反射器。
tau包含了K反射器的系数。
C在内存中的布局为LDC×M,其中上层的M×N块将用于存放结果QC。
为了让C在返回时保持Q,C必须是一个M x M的正方形矩阵,初始化为identity,即对角线条目都是1。
你可以考虑使用为ublas提供的lapack数字绑定,如在 (http:/boost.2283326.n4.nabble.com如何正确使用qr分解-td2710159.html。
)
然而,这个项目现在可能已经不存在了,或者已经休息了。让我们从第一原则重新开始。 我们的目标是解决Ax=b,或者至少要尽量减少/www.netlib.nonetliblapackdoubledgeqrf.f。colsA=rowsx
)rowsA=rowsb
DGEQRF( rowsA, colsA, A, rowsA, TAU, WORK, LWORK, INFO )
乘以 A
以获得 colsA<=rowsA
如同
(
http:/www.netlib.nonetliblapackdoubledormqr.fQ*R=A
)DORMQR( 'L', 'T', rowsA, 1, colsA, A, rowsA, TAU, b, rowsA, WORK, LWORK, INFO )使用右上角的反向替换。
(
http:/www.netlib.nonetliblapackdoubledtrtrs.fQT
)QT*b
DTRTRS( 'U', 'N', 'N', colsA, 1, A, rowsA, b, rowsA, INFO )R*x=QT*b
现在是第一个 的条目 包含解决方案向量
. 其余条目在索引colsA+1处和之后的欧氏法线为误差。