Main Page   Namespace List   Compound List   File List   Compound Members

# jama_qr.h

Go to the documentation of this file.
```00001 #ifndef JAMA_QR_H
00002 #define JAMA_QR_H
00003
00004 #include "tnt_array1d.h"
00005 #include "tnt_array2d.h"
00006 #include "tnt_math_utils.h"
00007
00008 namespace JAMA
00009 {
00010
00033
00034 template <class Real>
00035 class QR {
00036
00037
00041
00042    TNT::Array2D<Real> QR_;
00043
00048    int m, n;
00049
00053    TNT::Array1D<Real> Rdiag;
00054
00055
00056 public:
00057
00063         QR(const TNT::Array2D<Real> &A)         /* constructor */
00064         {
00065       QR_ = A.copy();
00066       m = A.dim1();
00067       n = A.dim2();
00068       Rdiag = TNT::Array1D<Real>(n);
00069           int i=0, j=0, k=0;
00070
00071       // Main loop.
00072       for (k = 0; k < n; k++) {
00073          // Compute 2-norm of k-th column without under/overflow.
00074          Real nrm = 0;
00075          for (i = k; i < m; i++) {
00076             nrm = hypot(nrm,QR_[i][k]);
00077          }
00078
00079          if (nrm != 0.0) {
00080             // Form k-th Householder vector.
00081             if (QR_[k][k] < 0) {
00082                nrm = -nrm;
00083             }
00084             for (i = k; i < m; i++) {
00085                QR_[i][k] /= nrm;
00086             }
00087             QR_[k][k] += 1.0;
00088
00089             // Apply transformation to remaining columns.
00090             for (j = k+1; j < n; j++) {
00091                Real s = 0.0;
00092                for (i = k; i < m; i++) {
00093                   s += QR_[i][k]*QR_[i][j];
00094                }
00095                s = -s/QR_[k][k];
00096                for (i = k; i < m; i++) {
00097                   QR_[i][j] += s*QR_[i][k];
00098                }
00099             }
00100          }
00101          Rdiag[k] = -nrm;
00102       }
00103    }
00104
00105
00111         int isFullRank() const
00112         {
00113       for (int j = 0; j < n; j++)
00114           {
00115          if (Rdiag[j] == 0)
00116             return 0;
00117       }
00118       return 1;
00119         }
00120
00121
00122
00123
00129
00130    TNT::Array2D<Real> getHouseholder (void)  const
00131    {
00132           TNT::Array2D<Real> H(m,n);
00133
00134           /* note: H is completely filled in by algorithm, so
00135              initializaiton of H is not necessary.
00136           */
00137       for (int i = 0; i < m; i++)
00138           {
00139          for (int j = 0; j < n; j++)
00140                  {
00141             if (i >= j) {
00142                H[i][j] = QR_[i][j];
00143             } else {
00144                H[i][j] = 0.0;
00145             }
00146          }
00147       }
00148           return H;
00149    }
00150
00151
00152
00156
00157         TNT::Array2D<Real> getR() const
00158         {
00159       TNT::Array2D<Real> R(n,n);
00160       for (int i = 0; i < n; i++) {
00161          for (int j = 0; j < n; j++) {
00162             if (i < j) {
00163                R[i][j] = QR_[i][j];
00164             } else if (i == j) {
00165                R[i][j] = Rdiag[i];
00166             } else {
00167                R[i][j] = 0.0;
00168             }
00169          }
00170       }
00171           return R;
00172    }
00173
00174
00175
00176
00177
00182
00183         TNT::Array2D<Real> getQ() const
00184         {
00185           int i=0, j=0, k=0;
00186
00187           TNT::Array2D<Real> Q(m,n);
00188       for (k = n-1; k >= 0; k--) {
00189          for (i = 0; i < m; i++) {
00190             Q[i][k] = 0.0;
00191          }
00192          Q[k][k] = 1.0;
00193          for (j = k; j < n; j++) {
00194             if (QR_[k][k] != 0) {
00195                Real s = 0.0;
00196                for (i = k; i < m; i++) {
00197                   s += QR_[i][k]*Q[i][j];
00198                }
00199                s = -s/QR_[k][k];
00200                for (i = k; i < m; i++) {
00201                   Q[i][j] += s*QR_[i][k];
00202                }
00203             }
00204          }
00205       }
00206           return Q;
00207    }
00208
00209
00216
00217    TNT::Array1D<Real> solve(const TNT::Array1D<Real> &b) const
00218    {
00219           if (b.dim1() != m)            /* arrays must be conformant */
00220                 return TNT::Array1D<Real>();
00221
00222           if ( !isFullRank() )          /* matrix is rank deficient */
00223           {
00224                 return TNT::Array1D<Real>();
00225           }
00226
00227           TNT::Array1D<Real> x = b;
00228           int i=0, j=0, k=0;
00229
00230       // Compute Y = transpose(Q)*b
00231       for (k = 0; k < n; k++)
00232           {
00233             Real s = 0.0;
00234             for (i = k; i < m; i++)
00235                         {
00236                s += QR_[i][k]*x[i];
00237             }
00238             s = -s/QR_[k][k];
00239             for (i = k; i < m; i++)
00240                         {
00241                x[i] += s*QR_[i][k];
00242             }
00243       }
00244       // Solve R*X = Y;
00245       for (k = n-1; k >= 0; k--)
00246           {
00247          x[k] /= Rdiag[k];
00248          for (i = 0; i < k; i++) {
00249                x[i] -= x[k]*QR_[i][k];
00250          }
00251       }
00252
00253
00254           /* return n x nx portion of X */
00255           TNT::Array1D<Real> x_(n);
00256           for (i=0; i<n; i++)
00257                         x_[i] = x[i];
00258
00259           return x_;
00260    }
00261
00268
00269    TNT::Array2D<Real> solve(const TNT::Array2D<Real> &B) const
00270    {
00271           if (B.dim1() != m)            /* arrays must be conformant */
00272                 return TNT::Array2D<Real>(0,0);
00273
00274           if ( !isFullRank() )          /* matrix is rank deficient */
00275           {
00276                 return TNT::Array2D<Real>(0,0);
00277           }
00278
00279       int nx = B.dim2();
00280           TNT::Array2D<Real> X = B;
00281           int i=0, j=0, k=0;
00282
00283       // Compute Y = transpose(Q)*B
00284       for (k = 0; k < n; k++) {
00285          for (j = 0; j < nx; j++) {
00286             Real s = 0.0;
00287             for (i = k; i < m; i++) {
00288                s += QR_[i][k]*X[i][j];
00289             }
00290             s = -s/QR_[k][k];
00291             for (i = k; i < m; i++) {
00292                X[i][j] += s*QR_[i][k];
00293             }
00294          }
00295       }
00296       // Solve R*X = Y;
00297       for (k = n-1; k >= 0; k--) {
00298          for (j = 0; j < nx; j++) {
00299             X[k][j] /= Rdiag[k];
00300          }
00301          for (i = 0; i < k; i++) {
00302             for (j = 0; j < nx; j++) {
00303                X[i][j] -= X[k][j]*QR_[i][k];
00304             }
00305          }
00306       }
00307
00308
00309           /* return n x nx portion of X */
00310           TNT::Array2D<Real> X_(n,nx);
00311           for (i=0; i<n; i++)
00312                 for (j=0; j<nx; j++)
00313                         X_[i][j] = X[i][j];
00314
00315           return X_;
00316    }
00317
00318
00319 };
00320
00321
00322 }
00323 // namespace JAMA
00324
00325 #endif
00326 // JAMA_QR__H
00327
```

Generated at Mon Jan 20 07:47:17 2003 for JAMA/C++ by 1.2.5 written by Dimitri van Heesch, © 1997-2001