CppAD: A C++ Algorithmic Differentiation Package 20110419
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-08 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common 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 $syntax%%sign% = LuFactor(%ip%, %jp%, %LU%)%$$
00048 
00049 
00050 $head Description$$
00051 Computes an LU factorization of the matrix $italic A$$ 
00052 where $italic 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 $italic sign$$ has prototype
00072 $syntax%
00073         int %sign%
00074 %$$
00075 If $italic A$$ is invertible, $italic sign$$ is plus or minus one
00076 and is the sign of the permutation corresponding to the row ordering
00077 $italic ip$$ and column ordering $italic jp$$.
00078 If $italic A$$ is not invertible, $italic sign$$ is zero.
00079 
00080 $head ip$$
00081 The argument $italic ip$$ has prototype
00082 $syntax%
00083         %SizeVector% &%ip%
00084 %$$
00085 (see description of $xref/LuFactor/SizeVector/SizeVector/$$ below).
00086 The size of $italic ip$$ is referred to as $italic n$$ in the
00087 specifications below.
00088 The input value of the elements of $italic ip$$ does not matter.
00089 The output value of the elements of $italic ip$$ determine
00090 the order of the rows in the permuted matrix.
00091 
00092 $head jp$$
00093 The argument $italic jp$$ has prototype
00094 $syntax%
00095         %SizeVector% &%jp%
00096 %$$
00097 (see description of $xref/LuFactor/SizeVector/SizeVector/$$ below).
00098 The size of $italic jp$$ must be equal to $italic n$$.
00099 The input value of the elements of $italic jp$$ does not matter.
00100 The output value of the elements of $italic jp$$ determine
00101 the order of the columns in the permuted matrix.
00102 
00103 $head LU$$
00104 The argument $italic LU$$ has the prototype
00105 $syntax%
00106         %FloatVector% &%LU%
00107 %$$
00108 and the size of $italic LU$$ must equal $latex n * n$$
00109 (see description of $xref/LuFactor/FloatVector/FloatVector/$$ below).
00110 
00111 $subhead A$$
00112 We define $italic A$$ as the matrix corresponding to the input 
00113 value of $italic LU$$.
00114 
00115 $subhead P$$
00116 We define the permuted matrix $italic P$$ in terms of $italic A$$ by
00117 $syntax%
00118         %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ]
00119 %$$
00120 
00121 $subhead L$$
00122 We define the lower triangular matrix $italic L$$ in terms of the 
00123 output value of $italic LU$$.
00124 The matrix $italic L$$ is zero above the diagonal
00125 and the rest of the elements are defined by
00126 $syntax%
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 $italic U$$ in terms of the
00133 output value of $italic LU$$.
00134 The matrix $italic U$$ is zero below the diagonal,
00135 one on the diagonal,
00136 and the rest of the elements are defined by
00137 $syntax%
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 $italic sign$$ is non-zero,
00144 $syntax%
00145         %L% * %U% = %P%
00146 %$$
00147 If the return value of $italic sign$$ is zero,
00148 the contents of $italic L$$ and $italic U$$ are not defined. 
00149 
00150 $subhead Determinant$$
00151 $index determinant$$
00152 If the return value $italic sign$$ is zero,
00153 the determinant of $italic A$$ is zero.
00154 If $italic sign$$ is non-zero,
00155 using the output value of $italic LU$$
00156 the determinant of the matrix $italic A$$ is equal to
00157 $syntax%
00158 %sign% * %LU%[%ip%[0], %jp%[0]] * %...% * %LU%[%ip%[%n%-1], %jp%[%n%-1]] 
00159 %$$
00160 
00161 $head SizeVector$$
00162 The type $italic SizeVector$$ must be a $xref/SimpleVector/$$ class with
00163 $xref/SimpleVector/Elements of Specified Type/elements of type size_t/$$.
00164 The routine $xref/CheckSimpleVector/$$ will generate an error message
00165 if this is not the case.
00166 
00167 $head FloatVector$$
00168 The type $italic FloatVector$$ must be a 
00169 $xref/SimpleVector//simple vector class/$$.
00170 The routine $xref/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 $italic FloatVector$$.
00176 The type $italic Float$$ must satisfy the conditions
00177 for a $xref/NumericType/$$ type.
00178 The routine $xref/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 $italic Float$$ objects $italic x$$ and $italic y$$:
00182 
00183 $table
00184 $bold Operation$$ $cnext $bold Description$$  $rnext
00185 $syntax%log(%x%)%$$ $cnext
00186         returns the logarithm of $italic x$$ as a $italic Float$$ object
00187 $tend
00188 
00189 $head AbsGeq$$
00190 Including the file $code lu_factor.hpp$$ defines the template function 
00191 $syntax%
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 $italic x$$ is greater than or equal the absolute value of $italic 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 $italic Float$$ objects. 
00201 If this operator is not defined for your use of $italic 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 $syntax%
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 $italic x$$ is greater than or equal the 
00218 sum of the square of the real and imaginary parts of $italic y$$. 
00219 
00220 $children%
00221         example/lu_factor.cpp%
00222         omh/lu_factor_hpp.omh
00223 %$$
00224 $head Example$$
00225 The file 
00226 $xref/LuFactor.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 $xref/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 PROGRAM
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                 jp.size() == n,
00309                 "Error in LuFactor: jp must have size equal to n"
00310         );
00311         CPPAD_ASSERT_KNOWN(
00312                 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 PROGRAM
00400 # endif