Content |
QR Factorization (Complex Variant)
In this example we again compute the \(QR\) 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()
{
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::qrf(A, tau);
//lapack::qrf(A, tau, work);
lapack::unmqr(Left, ConjTrans, A, tau, b);
// lapack::ungqr(A, tau, b, work);
const auto R = A.upper();
blas::sv(NoTrans, R, b);
cout << "x = " << b << endl;
}
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
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::qrf(A, tau);
//lapack::qrf(A, tau, work);
lapack::unmqr(Left, ConjTrans, A, tau, b);
// lapack::ungqr(A, tau, b, work);
const auto R = A.upper();
blas::sv(NoTrans, R, b);
cout << "x = " << b << endl;
}
Comments on Example Code
Compute the factorization \(A = QR\). Note that the workspace gets created implicitly and temporarily. So you might not want to do this inside a loop.
lapack::qrf(A, tau);
//lapack::qrf(A, tau, work);
//lapack::qrf(A, tau, work);
Compute \(\tilde{b} = Q^T b\). Note that we do not setup \(A\) explicitly.
lapack::unmqr(Left, ConjTrans, A, tau, b);
// lapack::ungqr(A, tau, b, work);
// lapack::ungqr(A, tau, b, work);
Solve \(R x = \tilde{b}\). Vector \(b\) gets overwritten with \(x\).
const auto R = A.upper();
blas::sv(NoTrans, R, b);
blas::sv(NoTrans, R, b);
Compile
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -o lapack-complex-unmqr lapack-complex-unmqr.cc
Run
$shell> cd flens/examples $shell> ./lapack-complex-unmqr 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,-2.11698e-15) (9,1.96902e-16) (9,-4.7151e-15) (9,-5.27794e-15)