Content |
QR Factorization (Complex Variant)
In this example we again compute the \(QR\) factorization and use it for solving a system of linear equations. In this example we setup matrix \(Q\) explicitly. See lapack-unmqr 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;
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::qrf(A, tau);
//lapack::qrf(A, tau, work);
Q = A;
lapack::ungqr(Q, tau);
//lapack::orgqr(Q, tau, work);
cout << "Q = " << Q << endl;
DenseVector<Array<Complex> > x;
x = conjTrans(Q)*b;
const auto R = A.upper();
blas::sv(NoTrans, R, x);
cout << "x = " << x << 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), 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::qrf(A, tau);
//lapack::qrf(A, tau, work);
Q = A;
lapack::ungqr(Q, tau);
//lapack::orgqr(Q, tau, work);
cout << "Q = " << Q << endl;
DenseVector<Array<Complex> > x;
x = conjTrans(Q)*b;
const auto R = A.upper();
blas::sv(NoTrans, R, x);
cout << "x = " << x << 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);
Explicitly setup \(Q\).
Q = A;
lapack::ungqr(Q, tau);
//lapack::orgqr(Q, tau, work);
lapack::ungqr(Q, tau);
//lapack::orgqr(Q, tau, work);
Compute \(\tilde{b} = Q^T b\).
DenseVector<Array<Complex> > x;
x = conjTrans(Q)*b;
x = conjTrans(Q)*b;
Solve \(R x = \tilde{b}\). Vector \(b\) gets overwritten with \(x\).
const auto R = A.upper();
blas::sv(NoTrans, R, x);
blas::sv(NoTrans, R, x);
Compile
$shell> cd flens/examples $shell> clang++ -DUSE_CXXLAPACK -framework vecLib -std=c++11 -Wall -I../.. -o lapack-complex-ungqr lapack-complex-ungqr.cc
Run
$shell> cd flens/examples $shell> ./lapack-complex-ungqr 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.258199) (1.249e-16,0.182574) (-3.33067e-16,0.208237) (1.55431e-15,-0.925547) (2.77556e-17,0.774597) (8.32667e-17,1.38778e-17) (-2.77556e-17,-0.535468) (0,-0.336563) (0,-0.258199) (-1.38778e-17,-0.912871) (-6.245e-17,-0.267734) (4.16334e-16,-0.168281) (0,-0.516398) (-2.77556e-17,0.365148) (9.71445e-17,-0.773453) (-1.66533e-16,0.0420703) x = (1,8.95832e-15) (9,-4.39219e-16) (9,2.03785e-14) (9,2.38827e-14)