QR algorithm

QR algorithm

[Karaoulis Marios]

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