# QR Factorization with Column Pivoting

In this example we compute a $$QR$$ factorization with column pivoting and use it for solving a system of linear equations. Another (more interesting) of this method are rank deficient least square problems.

For a general $$m \times n$$ matrix $$A$$ the FLENS-LAPACK function lapack::qp3 the factorization $$AP = QR$$ where $$P$$ is a permutation matrix, $$Q$$ a orthogonal/unitary matrix and $$R$$ an upper triangular or trapezoidal matrix.

## Example Code

#include <iostream>
#include <flens/flens.cxx>

using namespace std;
using namespace flens;

int
main()
{
typedef GeMatrix<FullStorage<double> >   GeMatrix;
typedef typename GeMatrix::IndexType     IndexType;
typedef DenseVector<Array<IndexType> >   IndexVector;
typedef DenseVector<Array<double> >      DenseVector;

GeMatrix        A(4,4);
IndexVector     jPiv;
DenseVector     tau;
DenseVector     b(4);
//DenseVector   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::qp3(A, jPiv, tau);
//lapack::qp3(A, jPiv, tau, work);

lapack::ormqr(Left, Trans, A, tau, b);
//lapack::ormqr(Left, Trans, A, tau, b, work);

blas::sv(NoTrans, A.upper(), b);

DenseVector     x(4);
for (IndexType i=1; i<=x.length(); ++i) {
x(jPiv(i)) = b(i);
}

cout << "x = " << x << endl;

cout << "jPiv = " << jPiv << endl;
}

Compute the factorization $$AP = QR$$. As vector jPiv has initially length zero it gets internally resized to length $$n$$ and initialized with Zero. This means that all columns are free columns. Also note that the workspace gets created implicitly and temporarily. So you might not want to do this inside a loop.

lapack::qp3(A, jPiv, tau);
//lapack::qp3(A, jPiv, tau, work);

Compute $$\tilde{b} = Q^H b$$. Vector $$b$$ gets overwritten with $$\tilde{b}$$.

lapack::ormqr(Left, Trans, A, tau, b);
//lapack::ormqr(Left, Trans, A, tau, b, work);

Solve $$R u = \tilde{b}$$. Vector $$b$$ gets overwritten with $$u =P^T x$$.

blas::sv(NoTrans, A.upper(), b);

Compute $$x = Pu$$. Note that we need an extra vector here!

DenseVector     x(4);
for (IndexType i=1; i<=x.length(); ++i) {
x(jPiv(i)) = b(i);
}

cout << "jPiv = " << jPiv << endl;
}

## Compile

$shell> cd flens/examples$shell> g++ -std=c++11 -Wall -I../.. -o lapack-geqp3 lapack-geqp3.cc


## Run

$shell> cd flens/examples$shell> ./lapack-geqp3
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
jPiv =
2              4              1              3


## Example with Complex Numbers

### Example Code

#include <iostream>
#include <flens/flens.cxx>

using namespace std;
using namespace flens;

int
main()
{
typedef complex<double>                  Complex;
const Complex                            I(0,1);

typedef GeMatrix<FullStorage<Complex> >  GeMatrix;
typedef typename GeMatrix::IndexType     IndexType;
typedef DenseVector<Array<IndexType> >   IndexVector;
typedef DenseVector<Array<Complex> >     DenseVector;

GeMatrix        A(4,4);
IndexVector     jPiv;
DenseVector     tau;
DenseVector     b(4);
//DenseVector   work;

A =  2,   3,  -1,   0,
-6,  -5,   0,   2,
2,  -5,   6,  -6,
4,   6,   2,  -3;

b = 20,
-33,
-43,
49;

A *= I;
b *= 2.0*I;

cout << "A = " << A << endl;
cout << "b = " << b << endl;

lapack::qp3(A, jPiv, tau);
//lapack::qp3(A, jPiv, tau, work);

lapack::unmqr(Left, ConjTrans, A, tau, b);
//lapack::unmqr(Left, ConjTrans, A, tau, b, work);

blas::sv(NoTrans, A.upper(), b);

DenseVector     x(4);
for (IndexType i=1; i<=x.length(); ++i) {
x(jPiv(i)) = b(i);
}

cout << "x = " << x << endl;

cout << "jPiv = " << jPiv << endl;
}

Make that $$A$$ and $$b$$ somehow look complex ;-)

A *= I;
b *= 2.0*I;

Compute the factorization $$AP = QR$$. As vector jPiv has initially length zero it gets internally resized to length $$n$$ and initialized with Zero. This means that all columns are free columns. Also note that the workspace gets created implicitly and temporarily. So you might not want to do this inside a loop.

lapack::qp3(A, jPiv, tau);
//lapack::qp3(A, jPiv, tau, work);

Compute $$\tilde{b} = Q^H b$$. Vector $$b$$ gets overwritten with $$\tilde{b}$$.

lapack::unmqr(Left, ConjTrans, A, tau, b);
//lapack::unmqr(Left, ConjTrans, A, tau, b, work);

Solve $$R u = \tilde{b}$$. Vector $$b$$ gets overwritten with $$u =P^T x$$.

blas::sv(NoTrans, A.upper(), b);

Compute $$x = Pu$$. Note that we need an extra vector here!

DenseVector     x(4);
for (IndexType i=1; i<=x.length(); ++i) {
x(jPiv(i)) = b(i);
}

cout << "jPiv = " << jPiv << endl;
}

### Compile

$shell> cd flens/examples$shell> g++ -DUSE_CXXLAPACK -framework vecLib -std=c++11 -Wall -I../.. -o lapack-complex-geqp3 lapack-complex-geqp3.cc


### Run

$shell> cd flens/examples$shell> ./lapack-complex-geqp3
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,40)                      (-0,-66)                      (-0,-86)                        (0,98)
x =
(2,1.29776e-14)              (18,8.39706e-16)              (18,4.76582e-14)              (18,5.05304e-14)
jPiv =
2              4              1              3