Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

tnt_qr.h

Go to the documentation of this file.
00001 /*
00002 *
00003 * Template Numerical Toolkit (SPUC): Linear Algebra Module
00004 *
00005 * Mathematical and Computational Sciences Division
00006 * National Institute of Technology,
00007 * Gaithersburg, MD USA
00008 *
00009 *
00010 * This software was developed at the National Institute of Standards and
00011 * Technology (NIST) by employees of the Federal Government in the course
00012 * of their official duties. Pursuant to title 17 Section 105 of the
00013 * United States Code, this software is not subject to copyright protection
00014 * and is in the public domain.  The Template Numerical Toolkit (SPUC) is
00015 * an experimental system.  NIST assumes no responsibility whatsoever for
00016 * its use by other parties, and makes no guarantees, expressed or implied,
00017 * about its quality, reliability, or any other characteristic.
00018 *
00019 * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
00020 * see http://math.nist.gov/tnt for latest updates.
00021 *
00022 */
00023 
00024 
00025 #ifndef QR_H
00026 #define QR_H
00027 
00028 // Classical QR factorization example, based on Stewart[1973].
00029 //
00030 //
00031 // This algorithm computes the factorization of a matrix A
00032 // into a product of an orthognal matrix (Q) and an upper triangular 
00033 // matrix (R), such that QR = A.
00034 //
00035 // Parameters:
00036 //
00037 //  A   (in):       Matrix(1:N, 1:N)
00038 //
00039 //  Q   (output):   Matrix(1:N, 1:N), collection of Householder
00040 //                      column vectors Q1, Q2, ... QN
00041 //
00042 //  R   (output):   upper triangular Matrix(1:N, 1:N)
00043 //
00044 // Returns:  
00045 //
00046 //  0 if successful, 1 if A is detected to be singular
00047 //
00048 
00049 
00050 #include <cmath>      //for sqrt() & fabs()
00051 #include "tntmath.h"  // for sign()
00052 
00053 // Classical QR factorization, based on Stewart[1973].
00054 //
00055 //
00056 // This algorithm computes the factorization of a matrix A
00057 // into a product of an orthognal matrix (Q) and an upper triangular 
00058 // matrix (R), such that QR = A.
00059 //
00060 // Parameters:
00061 //
00062 //  A   (in/out):  On input, A is square, Matrix(1:N, 1:N), that represents
00063 //                  the matrix to be factored.
00064 //
00065 //                 On output, Q and R is encoded in the same Matrix(1:N,1:N)
00066 //                 in the following manner:
00067 //
00068 //                  R is contained in the upper triangular section of A,
00069 //                  except that R's main diagonal is in D.  The lower
00070 //                  triangular section of A represents Q, where each
00071 //                  column j is the vector  Qj = I - uj*uj'/pi_j.
00072 //
00073 //  C  (output):    vector of Pi[j]
00074 //  D  (output):    main diagonal of R, i.e. D(i) is R(i,i)
00075 //
00076 // Returns:  
00077 //
00078 //  0 if successful, 1 if A is detected to be singular
00079 //
00080 
00081 namespace SPUC
00082 {
00083 
00084 template <class MaTRiX, class Vector>
00085 int QR_factor(MaTRiX &A, Vector& C, Vector &D)
00086 {
00087     assert(A.lbound() == 1);        // ensure these are all 
00088     assert(C.lbound() == 1);        // 1-based arrays and vectors
00089     assert(D.lbound() == 1);
00090 
00091     int M = A.dim1();
00092     int N = A.dim2(); 
00093 
00094     assert(M == N);                 // make sure A is square
00095 
00096     int i,j,k;
00097     typename MaTRiX::element_type eta, sigma, sum;
00098 
00099     // adjust the shape of C and D, if needed...
00100 
00101     if (N != C.size())  C.newsize(N);
00102     if (N != D.size())  D.newsize(N);
00103 
00104     for (k=1; k<N; k++)
00105     {
00106         // eta = max |M(i,k)|,  for k <= i <= n
00107         //
00108         eta = 0;
00109         for (i=k; i<=N; i++)
00110         {
00111             double absA = fabs(A(i,k));
00112             eta = ( absA >  eta ? absA : eta ); 
00113         }
00114 
00115         if (eta == 0)           // matrix is singular
00116         {
00117             cerr << "QR: k=" << k << "\n";
00118             return 1;
00119         }
00120 
00121         // form Qk and premiltiply M by it
00122         //
00123         for(i=k; i<=N; i++)
00124             A(i,k)  = A(i,k) / eta;
00125 
00126         sum = 0;
00127         for (i=k; i<=N; i++)
00128             sum = sum + A(i,k)*A(i,k);
00129         sigma = sign(A(k,k)) *  root(sum);
00130 
00131 
00132         A(k,k) = A(k,k) + sigma;
00133         C(k) = sigma * A(k,k);
00134         D(k) = -eta * sigma;
00135 
00136         for (j=k+1; j<=N; j++)
00137         {
00138             sum = 0;
00139             for (i=k; i<=N; i++)
00140                 sum = sum + A(i,k)*A(i,j);
00141             sum = sum / C(k);
00142 
00143             for (i=k; i<=N; i++)
00144                 A(i,j) = A(i,j) - sum * A(i,k);
00145         }
00146 
00147         D(N) = A(N,N);
00148     }
00149 
00150     return 0;
00151 }
00152 
00153 // modified form of upper triangular solve, except that the main diagonal
00154 // of R (upper portion of A) is in D.
00155 //
00156 template <class MaTRiX, class Vector>
00157 int R_solve(const MaTRiX &A, /*const*/ Vector &D, Vector &b)
00158 {
00159     assert(A.lbound() == 1);        // ensure these are all 
00160     assert(D.lbound() == 1);        // 1-based arrays and vectors
00161     assert(b.lbound() == 1);
00162 
00163     int i,j;
00164     int N = A.dim1();
00165 
00166     assert(N == A.dim2());
00167     assert(N == D.dim());
00168     assert(N == b.dim());
00169 
00170     typename MaTRiX::element_type sum;
00171 
00172     if (D(N) == 0)
00173         return 1;
00174 
00175     b(N) = b(N) / 
00176             D(N);
00177 
00178     for (i=N-1; i>=1; i--)
00179     {
00180         if (D(i) == 0)
00181             return 1;
00182         sum = 0;
00183         for (j=i+1; j<=N; j++)
00184             sum = sum + A(i,j)*b(j);
00185         b(i) = ( b(i) - sum ) / 
00186             D(i);
00187     }
00188 
00189     return 0;
00190 }
00191 
00192 
00193 template <class MaTRiX, class Vector>
00194 int QR_solve(const MaTRiX &A, const Vector &c, /*const*/ Vector &d, 
00195         Vector &b)
00196 {
00197     assert(A.lbound() == 1);        // ensure these are all 
00198     assert(c.lbound() == 1);        // 1-based arrays and vectors
00199     assert(d.lbound() == 1);
00200 
00201     int N=A.dim1();
00202 
00203     assert(N == A.dim2());
00204     assert(N == c.dim());
00205     assert(N == d.dim());
00206     assert(N == b.dim());
00207 
00208     int i,j;
00209     typename MaTRiX::element_type sum, tau;
00210 
00211     for (j=1; j<N; j++)
00212     {
00213         // form Q'*b
00214         sum = 0;
00215         for (i=j; i<=N; i++)
00216             sum = sum + A(i,j)*b(i);
00217         if (c(j) == 0)
00218             return 1;
00219         tau = sum / c(j);
00220        for (i=j; i<=N; i++)
00221             b(i) = b(i) - tau * A(i,j);
00222     }
00223     return R_solve(A, d, b);        // solve Rx = Q'b
00224 }
00225 
00226 } // namespace SPUC
00227 
00228 #endif
00229 // QR_H

Generated on Fri Sep 16 11:02:17 2005 for spuc by  doxygen 1.4.4