CppAD: A C++ Algorithmic Differentiation Package  20130918
lu_factor.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_LU_FACTOR_INCLUDED
00003 # define CPPAD_LU_FACTOR_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin LuFactor$$
00018 $escape #$$
00019 $spell
00020      cppad.hpp
00021      Cpp
00022      Geq
00023      Lu
00024      bool
00025      const
00026      ip
00027      jp
00028      namespace
00029      std
00030      typename
00031 $$
00032 
00033 $index LuFactor$$
00034 $index linear, Lu factor equation$$
00035 $index equation, Lu factor$$
00036 $index determinant, Lu factor$$
00037 $index solve, Lu factor$$
00038 
00039 $section LU Factorization of A Square Matrix$$
00040 
00041 $pre
00042 $$
00043 
00044 $head Syntax$$ $code# include <cppad/lu_factor.hpp>$$
00045 $pre
00046 $$
00047 $icode%sign% = LuFactor(%ip%, %jp%, %LU%)%$$
00048 
00049 
00050 $head Description$$
00051 Computes an LU factorization of the matrix $icode A$$ 
00052 where $icode A$$ is a square matrix.
00053 
00054 $head Include$$
00055 The file $code cppad/lu_factor.hpp$$ is included by $code cppad/cppad.hpp$$
00056 but it can also be included separately with out the rest of 
00057 the $code CppAD$$ routines.
00058 
00059 $head Matrix Storage$$
00060 All matrices are stored in row major order.
00061 To be specific, if $latex Y$$ is a vector
00062 that contains a $latex p$$ by $latex q$$ matrix,
00063 the size of $latex Y$$ must be equal to $latex  p * q $$ and for
00064 $latex i = 0 , \ldots , p-1$$,
00065 $latex j = 0 , \ldots , q-1$$,
00066 $latex \[
00067      Y_{i,j} = Y[ i * q + j ]
00068 \] $$
00069 
00070 $head sign$$
00071 The return value $icode sign$$ has prototype
00072 $codei%
00073      int %sign%
00074 %$$
00075 If $icode A$$ is invertible, $icode sign$$ is plus or minus one
00076 and is the sign of the permutation corresponding to the row ordering
00077 $icode ip$$ and column ordering $icode jp$$.
00078 If $icode A$$ is not invertible, $icode sign$$ is zero.
00079 
00080 $head ip$$
00081 The argument $icode ip$$ has prototype
00082 $codei%
00083      %SizeVector% &%ip%
00084 %$$
00085 (see description of $cref/SizeVector/LuFactor/SizeVector/$$ below).
00086 The size of $icode ip$$ is referred to as $icode n$$ in the
00087 specifications below.
00088 The input value of the elements of $icode ip$$ does not matter.
00089 The output value of the elements of $icode ip$$ determine
00090 the order of the rows in the permuted matrix.
00091 
00092 $head jp$$
00093 The argument $icode jp$$ has prototype
00094 $codei%
00095      %SizeVector% &%jp%
00096 %$$
00097 (see description of $cref/SizeVector/LuFactor/SizeVector/$$ below).
00098 The size of $icode jp$$ must be equal to $icode n$$.
00099 The input value of the elements of $icode jp$$ does not matter.
00100 The output value of the elements of $icode jp$$ determine
00101 the order of the columns in the permuted matrix.
00102 
00103 $head LU$$
00104 The argument $icode LU$$ has the prototype
00105 $codei%
00106      %FloatVector% &%LU%
00107 %$$
00108 and the size of $icode LU$$ must equal $latex n * n$$
00109 (see description of $cref/FloatVector/LuFactor/FloatVector/$$ below).
00110 
00111 $subhead A$$
00112 We define $icode A$$ as the matrix corresponding to the input 
00113 value of $icode LU$$.
00114 
00115 $subhead P$$
00116 We define the permuted matrix $icode P$$ in terms of $icode A$$ by
00117 $codei%
00118      %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ]
00119 %$$
00120 
00121 $subhead L$$
00122 We define the lower triangular matrix $icode L$$ in terms of the 
00123 output value of $icode LU$$.
00124 The matrix $icode L$$ is zero above the diagonal
00125 and the rest of the elements are defined by
00126 $codei%
00127      %L%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ]
00128 %$$
00129 for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , i$$.
00130 
00131 $subhead U$$
00132 We define the upper triangular matrix $icode U$$ in terms of the
00133 output value of $icode LU$$.
00134 The matrix $icode U$$ is zero below the diagonal,
00135 one on the diagonal,
00136 and the rest of the elements are defined by
00137 $codei%
00138      %U%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ]
00139 %$$
00140 for $latex i = 0 , \ldots , n-2$$ and $latex j = i+1 , \ldots , n-1$$.
00141 
00142 $subhead Factor$$
00143 If the return value $icode sign$$ is non-zero,
00144 $codei%
00145      %L% * %U% = %P%
00146 %$$
00147 If the return value of $icode sign$$ is zero,
00148 the contents of $icode L$$ and $icode U$$ are not defined. 
00149 
00150 $subhead Determinant$$
00151 $index determinant$$
00152 If the return value $icode sign$$ is zero,
00153 the determinant of $icode A$$ is zero.
00154 If $icode sign$$ is non-zero,
00155 using the output value of $icode LU$$
00156 the determinant of the matrix $icode A$$ is equal to
00157 $codei%
00158 %sign% * %LU%[%ip%[0], %jp%[0]] * %...% * %LU%[%ip%[%n%-1], %jp%[%n%-1]] 
00159 %$$
00160 
00161 $head SizeVector$$
00162 The type $icode SizeVector$$ must be a $cref SimpleVector$$ class with
00163 $cref/elements of type size_t/SimpleVector/Elements of Specified Type/$$.
00164 The routine $cref CheckSimpleVector$$ will generate an error message
00165 if this is not the case.
00166 
00167 $head FloatVector$$
00168 The type $icode FloatVector$$ must be a 
00169 $cref/simple vector class/SimpleVector/$$.
00170 The routine $cref CheckSimpleVector$$ will generate an error message
00171 if this is not the case.
00172 
00173 $head Float$$
00174 This notation is used to denote the type corresponding
00175 to the elements of a $icode FloatVector$$.
00176 The type $icode Float$$ must satisfy the conditions
00177 for a $cref NumericType$$ type.
00178 The routine $cref CheckNumericType$$ will generate an error message
00179 if this is not the case.
00180 In addition, the following operations must be defined for any pair
00181 of $icode Float$$ objects $icode x$$ and $icode y$$:
00182 
00183 $table
00184 $bold Operation$$ $cnext $bold Description$$  $rnext
00185 $codei%log(%x%)%$$ $cnext
00186      returns the logarithm of $icode x$$ as a $icode Float$$ object
00187 $tend
00188 
00189 $head AbsGeq$$
00190 Including the file $code lu_factor.hpp$$ defines the template function 
00191 $codei%
00192      template <typename %Float%>
00193      bool AbsGeq<%Float%>(const %Float% &%x%, const %Float% &%y%)
00194 %$$
00195 in the $code CppAD$$ namespace.
00196 This function returns true if the absolute value of 
00197 $icode x$$ is greater than or equal the absolute value of $icode y$$. 
00198 It is used by $code LuFactor$$ to choose the pivot elements.
00199 This template function definition uses the operator 
00200 $code <=$$ to obtain the absolute value for $icode Float$$ objects. 
00201 If this operator is not defined for your use of $icode Float$$,
00202 you will need to specialize this template so that it works for your
00203 use of $code LuFactor$$.
00204 $pre
00205 
00206 $$
00207 Complex numbers do not have the operation $code <=$$ defined.
00208 The specializations
00209 $codei%
00210 bool AbsGeq< std::complex<float> > 
00211      (const std::complex<float> &%x%, const std::complex<float> &%y%)
00212 bool AbsGeq< std::complex<double> > 
00213      (const std::complex<double> &%x%, const std::complex<double> &%y%)
00214 %$$ 
00215 are define by including $code lu_factor.hpp$$
00216 These return true if the sum of the square of the real and imaginary parts
00217 of $icode x$$ is greater than or equal the 
00218 sum of the square of the real and imaginary parts of $icode y$$. 
00219 
00220 $children%
00221      example/lu_factor.cpp%
00222      omh/lu_factor_hpp.omh
00223 %$$
00224 $head Example$$
00225 The file 
00226 $cref lu_factor.cpp$$
00227 contains an example and test of using $code LuFactor$$ by itself.
00228 It returns true if it succeeds and false otherwise.
00229 $pre
00230 
00231 $$
00232 The file $cref lu_solve.hpp$$ provides a useful example usage of 
00233 $code LuFactor$$ with $code LuInvert$$.
00234 
00235 $head Source$$
00236 The file $cref lu_factor.hpp$$ contains the 
00237 current source code that implements these specifications.
00238 
00239 $end
00240 --------------------------------------------------------------------------
00241 */
00242 // BEGIN C++
00243 
00244 # include <complex>
00245 # include <vector>
00246 
00247 # include <cppad/local/cppad_assert.hpp>
00248 # include <cppad/check_simple_vector.hpp>
00249 # include <cppad/check_numeric_type.hpp>
00250 
00251 namespace CppAD { // BEGIN CppAD namespace
00252 
00253 // AbsGeq
00254 template <typename Float>
00255 inline bool AbsGeq(const Float &x, const Float &y)
00256 {    Float xabs = x;
00257      if( xabs <= Float(0) )
00258           xabs = - xabs;
00259      Float yabs = y;
00260      if( yabs <= Float(0) )
00261           yabs = - yabs;
00262      return xabs >= yabs;
00263 }
00264 inline bool AbsGeq(
00265      const std::complex<double> &x, 
00266      const std::complex<double> &y)
00267 {    double xsq = x.real() * x.real() + x.imag() * x.imag();
00268      double ysq = y.real() * y.real() + y.imag() * y.imag();
00269 
00270      return xsq >= ysq;
00271 }
00272 inline bool AbsGeq(
00273      const std::complex<float> &x, 
00274      const std::complex<float> &y)
00275 {    float xsq = x.real() * x.real() + x.imag() * x.imag();
00276      float ysq = y.real() * y.real() + y.imag() * y.imag();
00277 
00278      return xsq >= ysq;
00279 }
00280 
00281 // Lines that are different from code in cppad/local/lu_ratio.hpp end with //
00282 template <class SizeVector, class FloatVector>                          //
00283 int LuFactor(SizeVector &ip, SizeVector &jp, FloatVector &LU)           //
00284 {    
00285      // type of the elements of LU                                   //
00286      typedef typename FloatVector::value_type Float;                 //
00287 
00288      // check numeric type specifications
00289      CheckNumericType<Float>();
00290 
00291      // check simple vector class specifications
00292      CheckSimpleVector<Float, FloatVector>();
00293      CheckSimpleVector<size_t, SizeVector>();
00294 
00295      size_t  i, j;          // some temporary indices
00296      const Float zero( 0 ); // the value zero as a Float object
00297      size_t  imax;          // row index of maximum element
00298      size_t  jmax;          // column indx of maximum element
00299      Float    emax;         // maximum absolute value
00300      size_t  p;             // count pivots
00301      int     sign;          // sign of the permutation
00302      Float   etmp;          // temporary element
00303      Float   pivot;         // pivot element
00304 
00305      // -------------------------------------------------------
00306      size_t n = ip.size();
00307      CPPAD_ASSERT_KNOWN(
00308           size_t(jp.size()) == n,
00309           "Error in LuFactor: jp must have size equal to n"
00310      );
00311      CPPAD_ASSERT_KNOWN(
00312           size_t(LU.size()) == n * n,
00313           "Error in LuFactor: LU must have size equal to n * m"
00314      );
00315      // -------------------------------------------------------
00316 
00317      // initialize row and column order in matrix not yet pivoted
00318      for(i = 0; i < n; i++)
00319      {    ip[i] = i;
00320           jp[i] = i;
00321      }
00322      // initialize the sign of the permutation
00323      sign = 1;
00324      // ---------------------------------------------------------
00325 
00326      // Reduce the matrix P to L * U using n pivots
00327      for(p = 0; p < n; p++)
00328      {    // determine row and column corresponding to element of 
00329           // maximum absolute value in remaining part of P
00330           imax = jmax = n;
00331           emax = zero;
00332           for(i = p; i < n; i++)
00333           {    for(j = p; j < n; j++)
00334                {    CPPAD_ASSERT_UNKNOWN(
00335                          (ip[i] < n) & (jp[j] < n)
00336                     );
00337                     etmp = LU[ ip[i] * n + jp[j] ];
00338 
00339                     // check if maximum absolute value so far
00340                     if( AbsGeq (etmp, emax) )
00341                     {    imax = i;
00342                          jmax = j;
00343                          emax = etmp;
00344                     }
00345                }
00346           }
00347           CPPAD_ASSERT_KNOWN( 
00348           (imax < n) & (jmax < n) ,
00349           "LuFactor can't determine an element with "
00350           "maximum absolute value.\n"
00351           "Perhaps original matrix contains not a number or infinity.\n" 
00352           "Perhaps your specialization of AbsGeq is not correct."
00353           );
00354           if( imax != p )
00355           {    // switch rows so max absolute element is in row p
00356                i        = ip[p];
00357                ip[p]    = ip[imax];
00358                ip[imax] = i;
00359                sign     = -sign;
00360           }
00361           if( jmax != p )
00362           {    // switch columns so max absolute element is in column p
00363                j        = jp[p];
00364                jp[p]    = jp[jmax];
00365                jp[jmax] = j;
00366                sign     = -sign;
00367           }
00368           // pivot using the max absolute element
00369           pivot   = LU[ ip[p] * n + jp[p] ];
00370 
00371           // check for determinant equal to zero
00372           if( pivot == zero )
00373           {    // abort the mission
00374                return   0;
00375           }
00376 
00377           // Reduce U by the elementary transformations that maps 
00378           // LU( ip[p], jp[p] ) to one.  Only need transform elements
00379           // above the diagonal in U and LU( ip[p] , jp[p] ) is
00380           // corresponding value below diagonal in L.
00381           for(j = p+1; j < n; j++)
00382                LU[ ip[p] * n + jp[j] ] /= pivot;
00383 
00384           // Reduce U by the elementary transformations that maps 
00385           // LU( ip[i], jp[p] ) to zero. Only need transform elements 
00386           // above the diagonal in U and LU( ip[i], jp[p] ) is 
00387           // corresponding value below diagonal in L.
00388           for(i = p+1; i < n; i++ )
00389           {    etmp = LU[ ip[i] * n + jp[p] ];
00390                for(j = p+1; j < n; j++)
00391                {    LU[ ip[i] * n + jp[j] ] -= 
00392                          etmp * LU[ ip[p] * n + jp[j] ];
00393                } 
00394           }
00395      }
00396      return sign;
00397 }
00398 } // END CppAD namespace 
00399 // END C++
00400 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines