Content |
Computing the QR Factorization
In this example we compute the \(QR\) factorization and use it for solving a system of linear equations.
Consider the system of linear equations \(Ax=b\). We store the coefficient matrix \(A\) and the right-hand side \(b\) in a single matrix \(Ab = (A,b)\). We then compute the \(QR\) factorization of \(Ab\) with lapack::qrf such that \(QR = (A,b).\) Hence, we have \(R = (Q^T A, Q^T b)\) where \(R\) is an upper trapezoidal matrix and \(Q^T A\) upper triangular. We finally use the triangular solver blas::sm to obtain the solution of \((Q^T A) x = (Q^T b)\).
Function lapack::qrf is the FLENS port of LAPACK's geqrf and blas::sm the port of the BLAS routine trsm.
Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
typedef double T;
int
main()
{
typedef GeMatrix<FullStorage<T> > GeMatrix;
typedef DenseVector<Array<T> > DenseVector;
typedef GeMatrix::IndexType IndexType;
const Underscore<IndexType> _;
const IndexType m = 4,
n = 5;
GeMatrix Ab(m, n);
Ab = 2, 3, -1, 0, 20,
-6, -5, 0, 2, -33,
2, -5, 6, -6, -43,
4, 6, 2, -3, 49;
cout << "Ab = " << Ab << endl;
DenseVector tau(std::min(m,n));
DenseVector work;
lapack::qrf(Ab, tau, work);
cout << "Ab = " << Ab << endl;
cout << "tau = " << tau << endl;
const auto A = Ab(_,_(1,m));
auto B = Ab(_,_(m+1,n));
blas::sm(Left, NoTrans, 1, A.upper(), B);
cout << "X = " << B << endl;
}
Comments on Example Code
With header flens.cxx all of FLENS gets included.
using namespace std;
using namespace flens;
typedef double T;
int
main()
{
Define some convenient typedefs for the matrix/vector types
typedef DenseVector<Array<T> > DenseVector;
Define an underscore operator for convenient matrix slicing
const Underscore<IndexType> _;
const IndexType m = 4,
n = 5;
GeMatrix Ab(m, n);
Ab = 2, 3, -1, 0, 20,
-6, -5, 0, 2, -33,
2, -5, 6, -6, -43,
4, 6, 2, -3, 49;
cout << "Ab = " << Ab << endl;
tau will contain the scalar factors of the elementary reflectors (see dgeqrf for details). Vector work is used as workspace and, if empty, gets resized by lapack::qrf to optimal size.
DenseVector work;
Compute the \(QR\) factorization of \(A\) with lapack::qrf the FLENS port of LAPACK's dgeqrf. Note: With lapack::qrf(Ab, tau) the workspace gets created internally and temporarily.
cout << "Ab = " << Ab << endl;
cout << "tau = " << tau << endl;
Solve the system of linear equation \(Ax =B\) using the triangular solver blas::sm which is the FLENS implementation of the BLAS routine trsm. Note that A.upper() creates a triangular view.
auto B = Ab(_,_(m+1,n));
blas::sm(Left, NoTrans, 1, A.upper(), B);
cout << "X = " << B << endl;
}
Compile
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -o lapack-geqrf lapack-geqrf.cc
Run
$shell> cd flens/examples $shell> ./lapack-geqrf Ab = 2 3 -1 0 20 -6 -5 0 2 -33 2 -5 6 -6 -43 4 6 2 -3 49 Ab = -7.74597 -6.45497 -2.32379 4.64758 -44.9266 -0.615639 -7.30297 4.9295 -4.38178 -60.7972 0.205213 -0.854313 -3.36155 2.85583 -4.55148 0.410426 0.260891 0.453851 0.210352 1.89316 tau = 1.2582 1.1124 1.6584 0 X = 1 9 9 9
Example with Implicit Workspace Creation
In this simplified variant of the above example the workspace gets created implicit. This simplifies the usage of FLENS-LAPACK routines. However, if you are doing this inside a loop it might lead to a considerable performance penalty.
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::qrf(A, tau);
//lapack::qrf(A, tau, work);
lapack::ormqr(Left, Trans, A, tau, b);
//lapack::ormqr(Left, Trans, A, tau, b, work);
blas::sv(NoTrans, A.upper(), 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, work);
Compute \(\tilde{b} = Q^H b\). Vector \(b\) gets overwritten with \(\tilde{b}\).
//lapack::ormqr(Left, Trans, A, tau, b, work);
Solve \(R x = \tilde{b}\). Vector \(b\) gets overwritten with \(x\).
Compile
$shell> cd flens/examples $shell> g++ -DUSE_CXXLAPACK -framework vecLib -std=c++11 -Wall -I../.. -o lapack-simple-geqrf lapack-simple-geqrf.cc
Run
$shell> cd flens/examples $shell> ./lapack-simple-geqrf 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;
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;
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::unmqr(Left, ConjTrans, A, tau, b);
blas::sv(NoTrans, A.upper(), b);
cout << "x = " << b << endl;
}
Comments on Example Code
Compute the factorization \(A = QR\).
Compute \(\tilde{b} = Q^H b\). Vector \(b\) gets overwritten with \(\tilde{b}\).
Solve \(R x = \tilde{b}\). Vector \(b\) gets overwritten with \(x\).
Compile
$shell> cd flens/examples $shell> g++ -DUSE_CXXLAPACK -framework vecLib -std=c++11 -Wall -I../.. -o lapack-complex-geqrf lapack-complex-geqrf.cc
Run
$shell> cd flens/examples $shell> ./lapack-complex-geqrf 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.3617e-14) (9,1.8825e-15) (9,-5.29762e-14) (9,-6.22797e-14)