Blender V2.61 - r43446

lu.h

Go to the documentation of this file.
00001 
00004 /*
00005 
00006 *
00007 * Template Numerical Toolkit (TNT): Linear Algebra Module
00008 *
00009 * Mathematical and Computational Sciences Division
00010 * National Institute of Technology,
00011 * Gaithersburg, MD USA
00012 *
00013 *
00014 * This software was developed at the National Institute of Standards and
00015 * Technology (NIST) by employees of the Federal Government in the course
00016 * of their official duties. Pursuant to title 17 Section 105 of the
00017 * United States Code, this software is not subject to copyright protection
00018 * and is in the public domain.  The Template Numerical Toolkit (TNT) is
00019 * an experimental system.  NIST assumes no responsibility whatsoever for
00020 * its use by other parties, and makes no guarantees, expressed or implied,
00021 * about its quality, reliability, or any other characteristic.
00022 *
00023 * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
00024 * see http://math.nist.gov/tnt for latest updates.
00025 *
00026 */
00027 
00028 
00029 
00030 #ifndef LU_H
00031 #define LU_H
00032 
00033 // Solve system of linear equations Ax = b.
00034 //
00035 //  Typical usage:
00036 //
00037 //    Matrix(double) A;
00038 //    Vector(Subscript) ipiv;
00039 //    Vector(double) b;
00040 //
00041 //    1)  LU_Factor(A,ipiv);
00042 //    2)  LU_Solve(A,ipiv,b);
00043 //
00044 //   Now b has the solution x.  Note that both A and b
00045 //   are overwritten.  If these values need to be preserved, 
00046 //   one can make temporary copies, as in 
00047 //
00048 //    O)  Matrix(double) T = A;
00049 //    1)  LU_Factor(T,ipiv);
00050 //    1a) Vector(double) x=b;
00051 //    2)  LU_Solve(T,ipiv,x);
00052 //
00053 //   See details below.
00054 //
00055 
00056 
00057 // for fabs() 
00058 //
00059 #include <cmath>
00060 
00061 // right-looking LU factorization algorithm (unblocked)
00062 //
00063 //   Factors matrix A into lower and upper  triangular matrices 
00064 //   (L and U respectively) in solving the linear equation Ax=b.  
00065 //
00066 //
00067 // Args:
00068 //
00069 // A        (input/output) Matrix(1:n, 1:n)  In input, matrix to be
00070 //                  factored.  On output, overwritten with lower and 
00071 //                  upper triangular factors.
00072 //
00073 // indx     (output) Vector(1:n)    Pivot vector. Describes how
00074 //                  the rows of A were reordered to increase
00075 //                  numerical stability.
00076 //
00077 // Return value:
00078 //
00079 // int      (0 if successful, 1 otherwise)
00080 //
00081 //
00082 
00083 
00084 namespace TNT
00085 {
00086 
00087 template <class MaTRiX, class VecToRSubscript>
00088 int LU_factor( MaTRiX &A, VecToRSubscript &indx)
00089 {
00090     assert(A.lbound() == 1);                // currently for 1-offset
00091     assert(indx.lbound() == 1);             // vectors and matrices
00092 
00093     Subscript M = A.num_rows();
00094     Subscript N = A.num_cols();
00095 
00096     if (M == 0 || N==0) return 0;
00097     if (indx.dim() != M)
00098         indx.newsize(M);
00099 
00100     Subscript i=0,j=0,k=0;
00101     Subscript jp=0;
00102 
00103     typename MaTRiX::element_type t;
00104 
00105     Subscript minMN =  (M < N ? M : N) ;        // min(M,N);
00106 
00107     for (j=1; j<= minMN; j++)
00108     {
00109 
00110         // find pivot in column j and  test for singularity.
00111 
00112         jp = j;
00113         t = fabs(A(j,j));
00114         for (i=j+1; i<=M; i++)
00115             if ( fabs(A(i,j)) > t)
00116             {
00117                 jp = i;
00118                 t = fabs(A(i,j));
00119             }
00120 
00121         indx(j) = jp;
00122 
00123         // jp now has the index of maximum element 
00124         // of column j, below the diagonal
00125 
00126         if ( A(jp,j) == 0 )                 
00127             return 1;       // factorization failed because of zero pivot
00128 
00129 
00130         if (jp != j)            // swap rows j and jp
00131             for (k=1; k<=N; k++)
00132             {
00133                 t = A(j,k);
00134                 A(j,k) = A(jp,k);
00135                 A(jp,k) =t;
00136             }
00137 
00138         if (j<M)                // compute elements j+1:M of jth column
00139         {
00140             // note A(j,j), was A(jp,p) previously which was
00141             // guarranteed not to be zero (Label #1)
00142             //
00143             typename MaTRiX::element_type recp =  1.0 / A(j,j);
00144 
00145             for (k=j+1; k<=M; k++)
00146                 A(k,j) *= recp;
00147         }
00148 
00149 
00150         if (j < minMN)
00151         {
00152             // rank-1 update to trailing submatrix:   E = E - x*y;
00153             //
00154             // E is the region A(j+1:M, j+1:N)
00155             // x is the column vector A(j+1:M,j)
00156             // y is row vector A(j,j+1:N)
00157 
00158             Subscript ii,jj;
00159 
00160             for (ii=j+1; ii<=M; ii++)
00161                 for (jj=j+1; jj<=N; jj++)
00162                     A(ii,jj) -= A(ii,j)*A(j,jj);
00163         }
00164     }
00165 
00166     return 0;
00167 }   
00168 
00169 
00170 
00171 
00172 template <class MaTRiX, class VecToR, class VecToRSubscripts>
00173 int LU_solve(const MaTRiX &A, const VecToRSubscripts &indx, VecToR &b)
00174 {
00175     assert(A.lbound() == 1);                // currently for 1-offset
00176     assert(indx.lbound() == 1);             // vectors and matrices
00177     assert(b.lbound() == 1);
00178 
00179     Subscript i,ii=0,ip,j;
00180     Subscript n = b.dim();
00181     typename MaTRiX::element_type sum = 0.0;
00182 
00183     for (i=1;i<=n;i++) 
00184     {
00185         ip=indx(i);
00186         sum=b(ip);
00187         b(ip)=b(i);
00188         if (ii)
00189             for (j=ii;j<=i-1;j++) 
00190                 sum -= A(i,j)*b(j);
00191         else if (sum) ii=i;
00192             b(i)=sum;
00193     }
00194     for (i=n;i>=1;i--) 
00195     {
00196         sum=b(i);
00197         for (j=i+1;j<=n;j++) 
00198             sum -= A(i,j)*b(j);
00199         b(i)=sum/A(i,i);
00200     }
00201 
00202     return 0;
00203 }
00204 
00205 } // namespace TNT
00206 
00207 #endif // LU_H
00208