Content |
LQ Factorization
In this example we again compute the \(LQ\) factorization and use it for solving a system of linear equations. However, in this example we do not setup matrix \(Q\) explicitly.
Example Code
#include <cxxstd/iostream.h>
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
int
main()
{
GeMatrix<FullStorage<double> > A(4,4), Q;
DenseVector<Array<double> > b(4);
DenseVector<Array<double> > tau;
//DenseVector<Array<double> > work;
A = 2, 3, -1, 0,
-6, -5, 0, 2,
2, -5, 6, -6,
4, 6, 2, -3;
b = 20,
-33,
-43,
49;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
lapack::lqf(A, tau);
//lapack::lqf(A, tau, work);
const auto L = A.lower();
blas::sv(NoTrans, L, b);
Q = A;
lapack::orglq(Q, tau);
//lapack::orglq(Q, tau, work);
cout << "Q = " << Q << endl;
DenseVector<Array<double> > x;
x = transpose(Q)*b;
cout << "x = " << x << endl;
}
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
int
main()
{
GeMatrix<FullStorage<double> > A(4,4), Q;
DenseVector<Array<double> > b(4);
DenseVector<Array<double> > tau;
//DenseVector<Array<double> > work;
A = 2, 3, -1, 0,
-6, -5, 0, 2,
2, -5, 6, -6,
4, 6, 2, -3;
b = 20,
-33,
-43,
49;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
lapack::lqf(A, tau);
//lapack::lqf(A, tau, work);
const auto L = A.lower();
blas::sv(NoTrans, L, b);
Q = A;
lapack::orglq(Q, tau);
//lapack::orglq(Q, tau, work);
cout << "Q = " << Q << endl;
DenseVector<Array<double> > x;
x = transpose(Q)*b;
cout << "x = " << x << endl;
}
Comments on Example Code
Compute the factorization \(A = LQ\). Note that the workspace gets created implicitly and temporarily. So you might not want to do this inside a loop.
lapack::lqf(A, tau);
//lapack::lqf(A, tau, work);
//lapack::lqf(A, tau, work);
Solve \(L u = b\). Vector \(b\) gets overwritten with \(u\).
const auto L = A.lower();
blas::sv(NoTrans, L, b);
blas::sv(NoTrans, L, b);
Explicitly setup \(Q\).
Q = A;
lapack::orglq(Q, tau);
//lapack::orglq(Q, tau, work);
lapack::orglq(Q, tau);
//lapack::orglq(Q, tau, work);
Compute \(x = Q^T u\).
DenseVector<Array<double> > x;
x = transpose(Q)*b;
x = transpose(Q)*b;
Compile
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -o lapack-orglq lapack-orglq.cc
Run
$shell> cd flens/examples $shell> ./lapack-orglq A = 2 3 -1 0 -6 -5 0 2 2 -5 6 -6 4 6 2 -3 b = 20 -33 -43 49 Q = -0.534522 -0.801784 0.267261 -0 0.595961 -0.218519 0.536365 -0.55623 0.564904 -0.386513 -0.0297318 0.728428 -0.2 0.4 0.8 0.4 x = 1 9 9 9