Content |
LQ Factorization (Complex Variant)
In this example we again compute the \(LQ\) factorization and use it for solving a system of linear equations. In this example we setup matrix \(Q\) explicitly. See lapack-gelqf for an example that avoids the explicit generation of \(Q\).
Example Code
#include <cxxstd/iostream.h>
#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), Q;
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(Q, tau, work);
blas::sv(NoTrans, A.lower(), b);
Q = A;
lapack::unglq(Q, tau);
//lapack::orglq(Q, tau, work);
cout << "Q = " << Q << endl;
DenseVector<Array<Complex> > x;
x = conjTrans(Q)*b;
cout << "x = " << x << endl;
}
#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), Q;
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(Q, tau, work);
blas::sv(NoTrans, A.lower(), b);
Q = A;
lapack::unglq(Q, tau);
//lapack::orglq(Q, tau, work);
cout << "Q = " << Q << endl;
DenseVector<Array<Complex> > x;
x = conjTrans(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(Q, tau, work);
// lapack::lqf(Q, tau, work);
Solve \(L u = b\). Vector \(b\) gets overwritten with \(u\).
blas::sv(NoTrans, A.lower(), b);
Explicitly setup \(Q\).
Q = A;
lapack::unglq(Q, tau);
//lapack::orglq(Q, tau, work);
lapack::unglq(Q, tau);
//lapack::orglq(Q, tau, work);
Compute \(x = Q^T u\). Vector \(b\) gets overwritten with \(x\).
DenseVector<Array<Complex> > x;
x = conjTrans(Q)*b;
x = conjTrans(Q)*b;
Compile
$shell> cd flens/examples $shell> clang++ -DUSE_CXXLAPACK -framework vecLib -std=c++11 -Wall -I../.. -o lapack-complex-unglq lapack-complex-unglq.cc
Run
$shell> cd flens/examples $shell> ./lapack-complex-unglq 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) Q = (0,-0.534522) (0,-0.801784) (0,0.267261) (-0,-0) (-5.55112e-17,-0.595961) (1.11022e-16,0.218519) (-1.11022e-16,-0.536365) (0,0.55623) (9.71445e-16,0.564904) (-8.88178e-16,-0.386513) (-2.22045e-16,-0.0297318) (8.60423e-16,0.728428) (-4.30211e-16,-0.2) (1.11022e-16,0.4) (-6.03684e-16,0.8) (-8.60423e-16,0.4) x = (1,-5.20006e-15) (9,1.6922e-15) (9,-5.23752e-15) (9,-8.81201e-15)