QR algorithm

Karaoulis Marios karaouli@geo.auth.gr
Tue May 30 08:05:00 GMT 2006

```Althought i understant that this is not the appropiate list to ask this,
but i have no other solution. I am currently using the qr function taken
from jama. This algorithm only produces the economy sized Q and R matrices.
If possible could somone please say the changes i must make in order to
produce the full ones?
The algorithm  is the above...

/*Based on JAMA 1.0.2 algorithm. Basically i only need the Q and R
matrices*/

/*NOTICE : This QR function is the economy sized one*/
/*  QR[mXn]     Q[mXn]  R[nXn]          */
void qr(double **QR,const int m,const int n,double **R,double **Q)
{

int k,i,j;
double nrm=0,s=0;
double hypot(double a,float b);

for (k = 0; k < n; k++) {
// Compute 2-norm of k-th column without under/overflow.
nrm = 0;
for (i = k; i < m; i++) {
nrm = hypot(nrm,QR[i][k]);
}

if (nrm != 0.0) {
// Form k-th Householder vector.
if (QR[k][k] < 0) {
nrm = -nrm;
}
for (i = k; i < m; i++) {
QR[i][k] /= nrm;
}
QR[k][k] += 1.0;

// Apply transformation to remaining columns.
for (j = k+1; j < n; j++) {
s = 0.0;
for ( i = k; i < m; i++) {
s += QR[i][k]*QR[i][j];
}
s = -s/QR[k][k];
for ( i = k; i < m; i++) {
QR[i][j] += s*QR[i][k];
}
}
}
Rdiag[k] = -nrm;
}

/*  Return the upper triangular factor*/
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
if (i < j) {
R[i][j] = QR[i][j];
} else if (i == j) {
R[i][j] = Rdiag[i];
} else {
R[i][j] = 0.0;
}
}
}

/* Generate and return the (economy-sized) orthogonal factor*/
for ( k = n-1; k >= 0; k--) {
for ( i = 0; i < m; i++) {
Q[i][k] = 0.0;
}
Q[k][k] = 1.0;
for ( j = k; j < n; j++) {
if (QR[k][k] != 0) {
s = 0.0;
for ( i = k; i < m; i++) {
s += QR[i][k]*Q[i][j];
}
s = -s/QR[k][k];
for ( i = k; i < m; i++) {
Q[i][j] += s*QR[i][k];
}
}
}
}

kill_float_vector(Rdiag,n+1); /*NOTICE : I don't need this anymore*/

return;
}

double hypot(double a,float b)
{
double r;
if (fabs(a) > fabs(b)) {
r = b/a;
r = fabs(a)*sqrt(1+r*r);
} else if (b != 0) {
r = a/b;
r = fabs(b)*sqrt(1+r*r);
} else {
r = 0.0;
}
return r;

}

--
*************************************
KARAOULIS MARIOS
Geophysical Laboratory
P.O. BOX 352-1
Aristotle University of Thessaloniki
MACEDONIA - GREECE
GR - 54124
------------------------------------------------