QR algorithm

[Karaoulis Marios <karaouli at geo dot auth dot gr>]

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