Content |
LQ Factorization
In this example we compute the \(LQ\) factorization and use it for solving a system of linear equations. In this example we do not setup matrix \(Q\) explicitly.
Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
typedef double T;
int
main()
{
GeMatrix<FullStorage<double> > A(4,4);
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);
blas::sv(NoTrans, A.lower(), b);
lapack::ormlq(Left, Trans, A, tau, b);
// lapack::ormlq(Left, Trans, A, tau, b, work);
cout << "x = " << b << 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, work);
Solve \(L u = b\). Vector \(b\) gets overwritten with \(u\).
Compute \(x = Q^T u\). Vector \(b\) gets overwritten with \(x\).
// lapack::ormlq(Left, Trans, A, tau, b, work);
Compile
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -o lapack-gelqf lapack-gelqf.cc
Run
$shell> cd flens/examples $shell> ./lapack-gelqf A = 2 3 -1 0 -6 -5 0 2 2 -5 6 -6 4 6 2 -3 b = 20 -33 -43 49 x = 1 9 9 9
Example with Complex Numbers
Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
typedef double T;
int
main()
{
typedef complex<double> Complex;
const Complex I(0,1);
GeMatrix<FullStorage<Complex> > A(4,4);
DenseVector<Array<Complex> > b(4);
DenseVector<Array<Complex> > tau;
//DenseVector<Array<Complex> > work;
A = 2, 3, -1, 0,
-6, -5, 0, 2,
2, -5, 6, -6,
4, 6, 2, -3;
A *=I;
b = 20,
-33,
-43,
49;
b *= I;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
lapack::lqf(A, tau);
// lapack::lqf(A, tau, work);
blas::sv(NoTrans, A.lower(), b);
lapack::unmlq(Left, ConjTrans, A, tau, b);
// lapack::unmlq(Left, ConjTrans, A, tau, b, work);
cout << "x = " << b << 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, work);
Solve \(L u = b\). Vector \(b\) gets overwritten with \(u\).
Compute \(x = Q^T u\). Vector \(b\) gets overwritten with \(x\).
// lapack::unmlq(Left, ConjTrans, A, tau, b, work);
Compile
$shell> cd flens/examples $shell> g++ -DUSE_CXXLAPACK -framework vecLib -std=c++11 -Wall -I../.. -o lapack-complex-gelqf lapack-complex-gelqf.cc
Run
$shell> cd flens/examples $shell> ./lapack-complex-gelqf A = (0,2) (0,3) (-0,-1) (0,0) (-0,-6) (-0,-5) (0,0) (0,2) (0,2) (-0,-5) (0,6) (-0,-6) (0,4) (0,6) (0,2) (-0,-3) b = (0,20) (-0,-33) (-0,-43) (0,49) x = (1,8.88178e-16) (9,-1.30045e-15) (9,-1.94289e-15) (9,-2.02977e-16)