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

...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

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