PDA

View Full Version : Chapter 2 QR Decomposition


MPD78
02-15-2010, 09:20 AM
Hello all,

In looking at the code for the struct, QRdcmp, the line containing the constructor for the matrix A is:

QRdcmp(MatDoub_I &a);

There isn't a MatDoub_O for the output. Does this mean that there isn't a way to obtain just the QR decomposition of the matrix A? (Without modifying the existing code.)

Here is my simple example program. This is modeled after this post.

http://www.nr.com/forum/showthread.php?t=1389&highlight=cholesky post #2 by davekw7x

#include "qrdcmp.h"

int main()
{
// Loop counters
Int i,j;

Int pause; // dummy variable

// Matrix size
const Int N = 3;

// Array for the constructor
const Doub a_d[N*N] = {
2.0, 1.0, 3.0,
-1.0, 0.0, 7.0,
0.0, -1.0, -1.0
};

MatDoub a(N,N,a_d);

// Perform the QR decomposition
QRdcmp qr(a);

// Send original matrix and output to the screen
cout << "The original matrix" << endl << endl;
cout << scientific << setprecision(4) << showpos;
for(i=0;i<N;i++){
for(j=0;j<N;j++){
cout << setw(15) << a_d[i * N +j];
}
cout << endl;
}

// If possible, this is where the code for the
// output of just the QR decomposition would go.

cin >> pause; // dummy variable to keep the output screen visible

return 0;
}

If it is possible to obtain just the QR decomposition, any help on how to do that would be greatly appreciated.


Thanks
Matt

davekw7x
02-15-2010, 03:12 PM
...If it is possible to obtain just the QR decomposition...
Matrices q-transpose and r are calculated in the constructor.


#include "../code/nr3.h"
#include "../code/qrdcmp.h"

//
// Driver for routine qrdcmp
// davekw7x
//

int main()
{
int i, j, k;
const int n = 3;
Doub a_array[n*n] = {
2.0, 1.0, 3.0,
-1.0, 0.0, 7.0,
0.0, -1.0, -1.0
};

MatDoub a(n, n, a_array);
MatDoub q(n, n), r(n, n), x(n, n);

cout << fixed << showpos;

//
// Print out the a matrix for comparison with the product of
// the Q and R decomposition matrices.
//
cout << "Original matrix:" << endl;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << a[i][j];
}
cout << endl;
}
cout << endl;

//
// The constructor performs the decomposition.
//
QRdcmp qrdcmp(a);

if (qrdcmp.sing) {
cerr << "Singularity in QR decomposition." << endl;
exit(EXIT_FAILURE);
}

//
// Get the Q and R matrices.
// Note that the the object gives the transpose of q, so, for
// pedagogical purposes, I transpose it when retrieving the values.
//
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
r[i][j] = qrdcmp.r[i][j];
q[i][j] = qrdcmp.qt[j][i];
}
}

//
// Compute the product of the Q and R matrices for comparison
// with the original matrix.
//
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
x[i][j] = 0.0;
for (k = 0; k < n; k++) {
x[i][j] += q[i][k] * r[k][j];
}
}
}

cout << "Q matrix of the decomposition:" << endl;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << q[i][j];
}
cout << endl;
}
cout << endl;

cout << "R matrix of the decomposition:" << endl;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << r[i][j];
}
cout << endl;
}
cout << endl;

cout << "Product of Q and R matrices:" << endl;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << x[i][j];
}
cout << endl;
}
cout << endl;

//
// Just for kicks, I'll calculate the absolute value of the determinant.
// It's just the absolute value of the product of the diagonal
// elements of r.
//
Doub det = 1.0;
for (i = 0; i < n; i++) {
det *= r[i][i];
}
det = fabs(det);

cout << scientific;
cout << "The absolute value of the determinant is "
<< det << endl;

return 0;
}


Output:

Original matrix:
+2.000000 +1.000000 +3.000000
-1.000000 +0.000000 +7.000000
+0.000000 -1.000000 -1.000000

Q matrix of the decomposition:
-0.894427 -0.182574 +0.408248
+0.447214 -0.365148 +0.816497
+0.000000 +0.912871 +0.408248

R matrix of the decomposition:
-2.236068 -0.894427 +0.447214
+0.000000 -1.095445 -4.016632
+0.000000 +0.000000 +6.531973

Product of Q and R matrices:
+2.000000 +1.000000 +3.000000
-1.000000 -0.000000 +7.000000
+0.000000 -1.000000 -1.000000

The absolute value of the determinant is 1.600000e+01



Regards,

Dave

MPD78
02-16-2010, 07:06 AM
Thanks alot davekw7x, all is working fine.

I haven't done it yet, but I am assuming, in a similar fashion I can get the just LU decomposition from the LU decomposition algorithm.

Thanks
Matt