jama_qr.h

Go to the documentation of this file.

#ifndef JAMA_QR_H
#define JAMA_QR_H
 
#include "../tnt/tnt_array1d.h"
#include "../tnt/tnt_array2d.h"
#include "../tnt/tnt_math_utils.h"
 
namespace JAMA
{
 
template <class Real>
class QR {
 
 
   TNT::Array2D<Real> QR_;
 
   int m, n;
 
   TNT::Array1D<Real> Rdiag;
 
 
public:
    
    QR(const TNT::Array2D<Real> &A)     /* constructor */
    {
      QR_ = A.copy();
      m = A.dim1();
      n = A.dim2();
      Rdiag = TNT::Array1D<Real>(n);
      int i=0, j=0, k=0;
 
      // Main loop.
      for (k = 0; k < n; k++) {
         // Compute 2-norm of k-th column without under/overflow.
         Real nrm = 0;
         for (i = k; i < m; i++) {
            nrm = TNT::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++) {
               Real 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;
      }
   }
 
 
    int isFullRank() const      
    {
      for (int j = 0; j < n; j++) 
      {
         if (Rdiag[j] == 0)
            return 0;
      }
      return 1;
    }
    
    
 
 
   TNT::Array2D<Real> getHouseholder (void)  const
   {
      TNT::Array2D<Real> H(m,n);
 
      /* note: H is completely filled in by algorithm, so
         initializaiton of H is not necessary.
      */
      for (int i = 0; i < m; i++) 
      {
         for (int j = 0; j < n; j++) 
         {
            if (i >= j) {
               H[i][j] = QR_[i][j];
            } else {
               H[i][j] = 0.0;
            }
         }
      }
      return H;
   }
 
 
 
    TNT::Array2D<Real> getR() const
    {
      TNT::Array2D<Real> R(n,n);
      for (int i = 0; i < n; i++) {
         for (int 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;
            }
         }
      }
      return R;
   }
    
    
 
 
 
    TNT::Array2D<Real> getQ() const
    {
      int i=0, j=0, k=0;
 
      TNT::Array2D<Real> Q(m,n);
      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) {
               Real 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];
               }
            }
         }
      }
      return Q;
   }
 
 
   TNT::Array1D<Real> solve(const TNT::Array1D<Real> &b) const
   {
      if (b.dim1() != m)        /* arrays must be conformant */
        return TNT::Array1D<Real>();
 
      if ( !isFullRank() )      /* matrix is rank deficient */
      {
        return TNT::Array1D<Real>();
      }
 
      TNT::Array1D<Real> x = b.copy();
 
      // Compute Y = transpose(Q)*b
      for (int k = 0; k < n; k++) 
      {
            Real s = 0.0; 
            for (int i = k; i < m; i++) 
            {
               s += QR_[i][k]*x[i];
            }
            s = -s/QR_[k][k];
            for (int i = k; i < m; i++) 
            {
               x[i] += s*QR_[i][k];
            }
      }
      // Solve R*X = Y;
      for (int k = n-1; k >= 0; k--) 
      {
         x[k] /= Rdiag[k];
         for (int i = 0; i < k; i++) {
               x[i] -= x[k]*QR_[i][k];
         }
      }
 
 
      /* return n x nx portion of X */
      TNT::Array1D<Real> x_(n);
      for (int i=0; i<n; i++)
            x_[i] = x[i];
 
      return x_;
   }
 
   TNT::Array2D<Real> solve(const TNT::Array2D<Real> &B) const
   {
      if (B.dim1() != m)        /* arrays must be conformant */
        return TNT::Array2D<Real>(0,0);
 
      if ( !isFullRank() )      /* matrix is rank deficient */
      {
        return TNT::Array2D<Real>(0,0);
      }
 
      int nx = B.dim2(); 
      TNT::Array2D<Real> X = B.copy();
      int i=0, j=0, k=0;
 
      // Compute Y = transpose(Q)*B
      for (k = 0; k < n; k++) {
         for (j = 0; j < nx; j++) {
            Real s = 0.0; 
            for (i = k; i < m; i++) {
               s += QR_[i][k]*X[i][j];
            }
            s = -s/QR_[k][k];
            for (i = k; i < m; i++) {
               X[i][j] += s*QR_[i][k];
            }
         }
      }
      // Solve R*X = Y;
      for (k = n-1; k >= 0; k--) {
         for (j = 0; j < nx; j++) {
            X[k][j] /= Rdiag[k];
         }
         for (i = 0; i < k; i++) {
            for (j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*QR_[i][k];
            }
         }
      }
 
 
      /* return n x nx portion of X */
      TNT::Array2D<Real> X_(n,nx);
      for (i=0; i<n; i++)
        for (j=0; j<nx; j++)
            X_[i][j] = X[i][j];
 
      return X_;
   }
 
 
};
 
 
}
// namespace JAMA
 
#endif
// JAMA_QR__H