C:/programs/SCPPNT/src/include/qr.h

Go to the documentation of this file.
00001 /*! \file qr.h
00002  \brief Definition of functions QR_factor, R_solve, and QR_solve.
00003  
00004  Contains modified version of the functions in qr.h from the
00005  Template Numerical Toolkit (TNT) using matrix iterators rather than
00006  matrix subscripting.
00007 
00008  */
00009 
00010 /*
00011 
00012  Simple C++ Numerical Toolkit (SCPPNT)
00013  http://www.smallwaters.com/software/cpp/scppnt.html
00014  This release updates original work contributed by 
00015  Brad Hanson (http://www.b-a-h.com/) in 2001.
00016 
00017  */
00018 
00019 // Modified version of qr.h from:
00020 /*
00021  *
00022  * Template Numerical Toolkit (TNT): Linear Algebra Module
00023  *
00024  * Mathematical and Computational Sciences Division
00025  * National Institute of Technology,
00026  * Gaithersburg, MD USA
00027  *
00028  *
00029  * This software was developed at the National Institute of Standards and
00030  * Technology (NIST) by employees of the Federal Government in the course
00031  * of their official duties. Pursuant to title 17 Section 105 of the
00032  * United States Code, this software is not subject to copyright protection
00033  * and is in the public domain.  The Template Numerical Toolkit (TNT) is
00034  * an experimental system.  NIST assumes no responsibility whatsoever for
00035  * its use by other parties, and makes no guarantees, expressed or implied,
00036  * about its quality, reliability, or any other characteristic.
00037  *
00038  * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
00039  * see http://math.nist.gov/tnt for latest updates.
00040  *
00041  */
00042 
00043 #ifndef SCPPNT_QR_H
00044 #define SCPPNT_QR_H
00045 
00046 #include <cmath>      //for sqrt() & fabs()
00047 #ifdef BOOST_NO_STDC_NAMESPACE
00048 namespace std
00049 { using ::fabs; using ::sqrt;}
00050 #endif
00051 
00052 #ifdef SCPPNT_NO_DIR_PREFIX
00053 #include "scppnt.h"
00054 #include "scppnt_error.h"
00055 #include "scppnt_math.h"  // for sign()
00056 #else
00057 #include "scppnt/scppnt.h"
00058 #include "scppnt/scppnt_error.h"
00059 #include "scppnt/scppnt_math.h"  // for sign()
00060 #endif
00061 
00062 namespace SCPPNT
00063 {
00064 
00065   /*!
00066    \brief Classical QR factorization, based on Stewart[1973].
00067    
00068    
00069    This algorithm computes the factorization of a matrix A
00070    into a product of an orthognal matrix (Q) and an upper triangular 
00071    matrix (R), such that QR = A.
00072    
00073    Returns:  
00074    
00075    0 if successful, 1 if A is detected to be singular
00076    
00077    Parameters:
00078    
00079    \param A   (in/out):  On input, A is square, Matrix(1:N, 1:N), that represents
00080    the matrix to be factored.
00081    On output, Q and R is encoded in the same Matrix(1:N,1:N)
00082    in the following manner:
00083    R is contained in the upper triangular section of A,
00084    except that R's main diagonal is in D.  The lower
00085    triangular section of A represents Q, where each
00086    column j is the vector  Qj = I - uj*uj'/pi_j.
00087    
00088    \param C  (output):    vector of Pi[j]
00089    
00090    \param D  (output):    main diagonal of R, i.e. D(i) is R(i,i)
00091    
00092    */
00093   template<class MaTRiX, class Vector> int QR_factor(MaTRiX &A, Vector& C, Vector &D)
00094   {
00095     //assert(A.lbound() == 1);        // ensure these are all
00096     if (A.lbound() != 1)
00097       throw BoundsError("SCPPNT::QR_factor");
00098 
00099     //assert(C.lbound() == 1);        // 1-based arrays and vectors
00100     if (C.lbound() != 1)
00101       throw BoundsError("SCPPNT::QR_factor");
00102 
00103     //assert(D.lbound() == 1);
00104     if (D.lbound() != 1)
00105       throw BoundsError("SCPPNT::R_factor");
00106 
00107     Subscript M = A.num_rows();
00108     Subscript N = A.num_columns();
00109 
00110     // assert(M == N);                 // make sure A is square
00111     if (M != N)
00112       throw BoundsError("QR_factor");
00113 
00114     Subscript i, j, k;
00115     typename MaTRiX::element_type eta, sigma, sum;
00116 
00117     // adjust the shape of C and D, if needed...
00118 
00119     if (N != C.size())
00120       C.newsize(N);
00121     if (N != D.size())
00122       D.newsize(N);
00123 
00124     typename MaTRiX::diag_iterator diagA = A.begin_diagonal(1, 1);
00125     typename Vector::iterator iC = C.begin();
00126     typename Vector::iterator iD = D.begin();
00127     typename MaTRiX::columns_iterator columnsA = A.begin_columns();
00128     typename Vector::iterator dn = D.begin() + N - 1;
00129     typename MaTRiX::row_iterator ann = A.begin_row(N) + N - 1;
00130 
00131     for (k=1; k<N; k++, ++columnsA, ++iC, ++iD)
00132     {
00133       // eta = max |M(i,k)|,  for k <= i <= n
00134       //
00135       eta = 0;
00136       typename MaTRiX::column_iterator columnA = *columnsA + k - 1;
00137 
00138       for (i=k; i<=N; i++, ++columnA)
00139       {
00140         double absA = std::fabs(*columnA); // double absA = std::fabs(A(i,k));
00141         eta = (absA > eta ? absA : eta );
00142       }
00143 
00144       if (eta == 0) // matrix is singular
00145       {
00146         // cerr << "QR: k=" << k << "\n";
00147         return 1;
00148       }
00149 
00150       // form Qk and premiltiply M by it
00151       //
00152       columnA = *columnsA + k - 1;
00153       for (i=k; i<=N; i++, ++columnA)
00154         *columnA /= eta; //A(i,k)  = A(i,k) / eta;
00155 
00156       sum = 0;
00157       columnA = *columnsA + k - 1;
00158       for (i=k; i<=N; i++, ++columnA)
00159         sum += *columnA * *columnA; //sum = sum + A(i,k)*A(i,k);
00160 
00161       sigma = sign(*diagA) * std::sqrt(sum); // sigma = sign(A(k,k)) *  sqrt(sum);
00162 
00163       *diagA += sigma; // A(k,k) = A(k,k) + sigma;
00164       *iC = sigma * *diagA; // C(k) = sigma * A(k,k);
00165       ++diagA;
00166       *iD = -eta * sigma; //D(k) = -eta * sigma;
00167 
00168       typename MaTRiX::columns_iterator columnsj = columnsA+1;
00169       for (j=k+1; j<=N; j++, ++columnsj)
00170       {
00171         sum = 0;
00172         columnA = *columnsA + k - 1;
00173         typename MaTRiX::column_iterator columnj = *columnsj + k - 1;
00174 
00175         for (i=k; i<=N; i++, ++columnj, ++columnA)
00176           sum += *columnA * *columnj; // sum = sum + A(i,k)*A(i,j);
00177 
00178         sum /= *iC; // sum = sum / C(k);
00179 
00180         columnA = *columnsA + k - 1;
00181         columnj = *columnsj + k - 1;
00182         for (i=k; i<=N; i++, ++columnA, ++columnj)
00183           *columnj -= sum * *columnA; // A(i,j) = A(i,j) - sum * A(i,k);
00184       }
00185 
00186       *dn = *ann; // D(N) = A(N,N);
00187     }
00188     return 0;
00189   }
00190 
00191   //! Modified form of upper triangular solve, except that the main diagonal of R (upper portion of A) is in D.
00192   template<class MaTRiX, class Vector> int R_solve(const MaTRiX &A, /*const*/Vector &D, Vector &b)
00193   {
00194     char *funcname = "SCPPNT::R_solve";
00195 
00196     // assert(A.lbound() == 1);        // ensure these are all 
00197     if (A.lbound() != 1)
00198       throw BoundsError(funcname);
00199 
00200     // assert(D.lbound() == 1);        // 1-based arrays and vectors
00201     if (D.lbound() != 1)
00202       throw BoundsError(funcname);
00203 
00204     // assert(b.lbound() == 1);
00205     if (b.lbound() != 1)
00206       throw BoundsError(funcname);
00207 
00208     Subscript i, j;
00209     Subscript N = A.num_rows();
00210 
00211     if (N != A.num_columns())
00212       throw BoundsError(funcname); // assert(N == A.num_columns());
00213     if (N != D.dim())
00214       throw BoundsError(funcname); // assert(N == D.dim());
00215     if (N != b.dim())
00216       throw BoundsError(funcname); // assert(N == b.dim());
00217 
00218     typename MaTRiX::element_type sum;
00219 
00220     if (D(N) == 0)
00221       return 1;
00222 
00223     b(N) = b(N) / D(N);
00224 
00225     typename MaTRiX::const_rows_iterator rowsA = A.begin_rows() + N - 2;
00226     typename Vector::iterator iD = D.begin() + N - 2;
00227     typename Vector::iterator ib = b.begin() + N - 2;
00228 
00229     for (i=N-1; i>=1; i--, --iD, --ib)
00230     {
00231       if (D(i) == 0)
00232         return 1;
00233 
00234       sum = 0;
00235       typename MaTRiX::const_row_iterator rowA = *rowsA + i;
00236       typename Vector::iterator jb = ib + 1;
00237 
00238       for (j=i+1; j<=N; j++, ++rowA, ++jb)
00239         sum += *rowA * *jb; // sum = sum + A(i,j)*b(j);
00240 
00241       --rowsA;
00242       *ib = (*ib - sum) / *iD; // b(i) = ( b(i) - sum ) / D(i);
00243     }
00244     return 0;
00245   }
00246 
00247   //! Solve Rx = Q'b, where A, c, and d are output from QR_factor
00248   template<class MaTRiX, class Vector> int QR_solve(const MaTRiX &A, const Vector &c, /*const*/
00249       Vector &d, Vector &b)
00250   {
00251     char *funcname = "SCPPNT::QR_solve";
00252 
00253     //assert(A.lbound() == 1);        // ensure these are all
00254     if (A.lbound() != 1)
00255       throw BoundsError(funcname);
00256     //assert(c.lbound() == 1);        // 1-based arrays and vectors
00257     if (c.lbound() != 1)
00258       throw BoundsError(funcname);
00259     if (d.lbound() != 1)
00260       throw BoundsError(funcname); // assert(d.lbound() == 1);
00261 
00262     Subscript N=A.num_rows();
00263 
00264     if (N != A.num_columns())
00265       throw BoundsError(funcname); //assert(N == A.num_columns());
00266     if (N != c.dim())
00267       throw BoundsError(funcname); //assert(N == c.dim());
00268     if (N != d.dim())
00269       throw BoundsError(funcname); //assert(N == d.dim());
00270     if (N != b.dim())
00271       throw BoundsError(funcname); //assert(N == b.dim());
00272 
00273     Subscript i, j;
00274     typename MaTRiX::element_type sum, tau;
00275     typename MaTRiX::const_columns_iterator columnsA = A.begin_columns();
00276     typename Vector::const_iterator cj = c.begin();
00277     for (j=1; j<N; j++, ++columnsA, ++cj)
00278     {
00279       // form Q'*b
00280       sum = 0;
00281       typename MaTRiX::const_column_iterator columnA = *columnsA + j - 1;
00282       typename Vector::iterator bi = b.begin() + j - 1;
00283 
00284       for (i=j; i<=N; i++, ++columnA, ++bi)
00285         sum += *columnA * *bi; // sum = sum + A(i,j)*b(i);
00286 
00287       if (*cj == 0) // if (c(j) == 0)
00288         return 1;
00289       tau = sum / *cj; // tau = sum / c(j);
00290       bi = b.begin() + j - 1;
00291       columnA = *columnsA + j - 1;
00292 
00293       for (i=j; i<=N; i++, ++bi, ++columnA)
00294         *bi -= tau * *columnA; // b(i) = b(i) - tau * A(i,j);
00295     }
00296     return R_solve(A, d, b); // solve Rx = Q'b
00297   }
00298 
00299 } // namespace SCPPNT
00300 
00301 #endif
00302 // SCPPNT_QR_H

Generated on Tue Dec 18 23:34:06 2007 for SCPPNT by  doxygen 1.5.4