This is the mail archive of the
gsl-discuss@sourceware.org
mailing list for the GSL project.
QR algorithm
- From: Karaoulis Marios <karaouli at geo dot auth dot gr>
- To: gsl <gsl-discuss at sourceware dot org>
- Date: Mon, 29 May 2006 19:15:33 +0300
- Subject: QR algorithm
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;
float *Rdiag=float_vector(n,"Rdiadg");
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
------------------------------------------------
e-mail address karaouli@geo.auth.gr
web site http://base_06.geo.auth.gr/
------------------------------------------------
Work: +30 2310 991409
Fax: +30 2310 998528
*************************************