CppAD: A C++ Algorithmic Differentiation Package  20130102
lu_solve.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_LU_SOLVE_INCLUDED
00003 # define CPPAD_LU_SOLVE_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 LuSolve$$
00018 $escape #$$
00019 $spell
00020      cppad.hpp
00021      det
00022      exp
00023      Leq
00024      typename
00025      bool
00026      const
00027      namespace
00028      std
00029      Geq
00030      Lu
00031      CppAD
00032      signdet
00033      logdet
00034 $$
00035 
00036 $index LuSolve$$
00037 $index linear, equation$$
00038 $index equation, linear$$
00039 $index determinant, Lu$$
00040 $index solve, linear equation$$
00041 
00042 $section Compute Determinant and Solve Linear Equations$$
00043 
00044 $pre
00045 $$
00046 
00047 $head Syntax$$ $code# include <cppad/lu_solve.hpp>$$
00048 $pre
00049 $$
00050 $icode%signdet% = LuSolve(%n%, %m%, %A%, %B%, %X%, %logdet%)%$$
00051 
00052 
00053 $head Description$$
00054 Use an LU factorization of the matrix $icode A$$ to
00055 compute its determinant 
00056 and solve for $icode X$$ in the linear of equation
00057 $latex \[
00058      A * X = B
00059 \] $$
00060 where $icode A$$ is an 
00061 $icode n$$ by $icode n$$ matrix,
00062 $icode X$$ is an 
00063 $icode n$$ by $icode m$$ matrix, and
00064 $icode B$$ is an $latex n x m$$ matrix.
00065 
00066 $head Include$$
00067 The file $code cppad/lu_solve.hpp$$ is included by $code cppad/cppad.hpp$$
00068 but it can also be included separately with out the rest of 
00069 the $code CppAD$$ routines.
00070 
00071 $head Factor and Invert$$
00072 This routine is an easy to user interface to
00073 $cref LuFactor$$ and $cref LuInvert$$ for computing determinants and
00074 solutions of linear equations.
00075 These separate routines should be used if
00076 one right hand side $icode B$$
00077 depends on the solution corresponding to another
00078 right hand side (with the same value of $icode A$$).
00079 In this case only one call to $code LuFactor$$ is required
00080 but there will be multiple calls to $code LuInvert$$.
00081 
00082 
00083 $head Matrix Storage$$
00084 All matrices are stored in row major order.
00085 To be specific, if $latex Y$$ is a vector
00086 that contains a $latex p$$ by $latex q$$ matrix,
00087 the size of $latex Y$$ must be equal to $latex  p * q $$ and for
00088 $latex i = 0 , \ldots , p-1$$,
00089 $latex j = 0 , \ldots , q-1$$,
00090 $latex \[
00091      Y_{i,j} = Y[ i * q + j ]
00092 \] $$
00093 
00094 $head signdet$$
00095 The return value $icode signdet$$ is a $code int$$ value
00096 that specifies the sign factor for the determinant of $icode A$$.
00097 This determinant of $icode A$$ is zero if and only if $icode signdet$$
00098 is zero. 
00099 
00100 $head n$$
00101 The argument $icode n$$ has type $code size_t$$ 
00102 and specifies the number of rows in the matrices
00103 $icode A$$,
00104 $icode X$$,
00105 and $icode B$$.
00106 The number of columns in $icode A$$ is also equal to $icode n$$.
00107 
00108 $head m$$
00109 The argument $icode m$$ has type $code size_t$$ 
00110 and specifies the number of columns in the matrices
00111 $icode X$$
00112 and $icode B$$.
00113 If $icode m$$ is zero,
00114 only the determinant of $icode A$$ is computed and
00115 the matrices $icode X$$ and $icode B$$ are not used.
00116 
00117 $head A$$
00118 The argument $icode A$$ has the prototype
00119 $codei%
00120      const %FloatVector% &%A%
00121 %$$
00122 and the size of $icode A$$ must equal $latex n * n$$
00123 (see description of $cref/FloatVector/LuSolve/FloatVector/$$ below).
00124 This is the $latex n$$ by $icode n$$ matrix that 
00125 we are computing the determinant of 
00126 and that defines the linear equation.
00127 
00128 $head B$$
00129 The argument $icode B$$ has the prototype
00130 $codei%
00131      const %FloatVector% &%B%
00132 %$$
00133 and the size of $icode B$$ must equal $latex n * m$$
00134 (see description of $cref/FloatVector/LuSolve/FloatVector/$$ below).
00135 This is the $latex n$$ by $icode m$$ matrix that 
00136 defines the right hand side of the linear equations.
00137 If $icode m$$ is zero, $icode B$$ is not used.
00138 
00139 $head X$$
00140 The argument $icode X$$ has the prototype
00141 $codei%
00142      %FloatVector% &%X%
00143 %$$
00144 and the size of $icode X$$ must equal $latex n * m$$
00145 (see description of $cref/FloatVector/LuSolve/FloatVector/$$ below).
00146 The input value of $icode X$$ does not matter.
00147 On output, the elements of $icode X$$ contain the solution
00148 of the equation we wish to solve
00149 (unless $icode signdet$$ is equal to zero).
00150 If $icode m$$ is zero, $icode X$$ is not used.
00151 
00152 $head logdet$$
00153 The argument $icode logdet$$ has prototype
00154 $codei%
00155      %Float% &%logdet%
00156 %$$
00157 On input, the value of $icode logdet$$ does not matter.
00158 On output, it has been set to the 
00159 log of the determinant of $icode A$$ 
00160 (but not quite).
00161 To be more specific,
00162 the determinant of $icode A$$ is given by the formula
00163 $codei%
00164      %det% = %signdet% * exp( %logdet% )
00165 %$$
00166 This enables $code LuSolve$$ to use logs of absolute values
00167 in the case where $icode Float$$ corresponds to a real number.
00168 
00169 $head Float$$
00170 The type $icode Float$$ must satisfy the conditions
00171 for a $cref NumericType$$ type.
00172 The routine $cref CheckNumericType$$ will generate an error message
00173 if this is not the case.
00174 In addition, the following operations must be defined for any pair
00175 of $icode Float$$ objects $icode x$$ and $icode y$$:
00176 
00177 $table
00178 $bold Operation$$ $cnext $bold Description$$  $rnext
00179 $codei%log(%x%)%$$ $cnext
00180      returns the logarithm of $icode x$$ as a $icode Float$$ object
00181 $tend
00182 
00183 $head FloatVector$$
00184 The type $icode FloatVector$$ must be a $cref SimpleVector$$ class with
00185 $cref/elements of type Float/SimpleVector/Elements of Specified Type/$$.
00186 The routine $cref CheckSimpleVector$$ will generate an error message
00187 if this is not the case.
00188 
00189 $head LeqZero$$
00190 Including the file $code lu_solve.hpp$$ defines the template function 
00191 $codei%
00192      template <typename %Float%>
00193      bool LeqZero<%Float%>(const %Float% &%x%)
00194 %$$
00195 in the $code CppAD$$ namespace.
00196 This function returns true if $icode x$$ is less than or equal to zero
00197 and false otherwise.
00198 It is used by $code LuSolve$$ to avoid taking the log of
00199 zero (or a negative number if $icode Float$$ corresponds to real numbers).
00200 This template function definition assumes that the operator 
00201 $code <=$$ is defined for $icode Float$$ objects. 
00202 If this operator is not defined for your use of $icode Float$$,
00203 you will need to specialize this template so that it works for your
00204 use of $code LuSolve$$.
00205 $pre
00206 
00207 $$
00208 Complex numbers do not have the operation or $code <=$$ defined.
00209 In addition, in the complex case, 
00210 one can take the log of a negative number.
00211 The specializations
00212 $codei%
00213      bool LeqZero< std::complex<float> > (const std::complex<float> &%x%)
00214      bool LeqZero< std::complex<double> >(const std::complex<double> &%x%)
00215 %$$ 
00216 are defined by including $code lu_solve.hpp$$.
00217 These return true if $icode x$$ is zero and false otherwise.
00218 
00219 $head AbsGeq$$
00220 Including the file $code lu_solve.hpp$$ defines the template function 
00221 $codei%
00222      template <typename %Float%>
00223      bool AbsGeq<%Float%>(const %Float% &%x%, const %Float% &%y%)
00224 %$$
00225 If the type $icode Float$$ does not support the $code <=$$ operation
00226 and it is not $code std::complex<float>$$ or $code std::complex<double>$$,
00227 see the documentation for $code AbsGeq$$ in $cref/LuFactor/LuFactor/AbsGeq/$$. 
00228 
00229 $children%
00230      example/lu_solve.cpp%
00231      omh/lu_solve_hpp.omh
00232 %$$
00233 $head Example$$
00234 The file 
00235 $cref lu_solve.cpp$$
00236 contains an example and test of using this routine.
00237 It returns true if it succeeds and false otherwise.
00238 
00239 $head Source$$
00240 The file $cref lu_solve.hpp$$ contains the
00241 current source code that implements these specifications.
00242 
00243 $end
00244 --------------------------------------------------------------------------
00245 */
00246 // BEGIN C++
00247 # include <complex>
00248 # include <vector>
00249 
00250 // link exp for float and double cases
00251 # include <cppad/base_require.hpp>
00252 
00253 # include <cppad/local/cppad_assert.hpp>
00254 # include <cppad/check_simple_vector.hpp>
00255 # include <cppad/check_numeric_type.hpp>
00256 # include <cppad/lu_factor.hpp>
00257 # include <cppad/lu_invert.hpp>
00258 
00259 namespace CppAD { // BEGIN CppAD namespace
00260 
00261 // LeqZero
00262 template <typename Float>
00263 inline bool LeqZero(const Float &x)
00264 {    return x <= Float(0); }
00265 inline bool LeqZero( const std::complex<double> &x )
00266 {    return x == std::complex<double>(0); }
00267 inline bool LeqZero( const std::complex<float> &x )
00268 {    return x == std::complex<float>(0); }
00269 
00270 // LuSolve
00271 template <typename Float, typename FloatVector>
00272 int LuSolve(
00273      size_t             n      ,
00274      size_t             m      , 
00275      const FloatVector &A      , 
00276      const FloatVector &B      , 
00277      FloatVector       &X      , 
00278      Float        &logdet      )
00279 {    
00280      // check numeric type specifications
00281      CheckNumericType<Float>();
00282 
00283      // check simple vector class specifications
00284      CheckSimpleVector<Float, FloatVector>();
00285 
00286      size_t        p;       // index of pivot element (diagonal of L)
00287      int     signdet;       // sign of the determinant
00288      Float     pivot;       // pivot element
00289 
00290      // the value zero
00291      const Float zero(0);
00292 
00293      // pivot row and column order in the matrix
00294      std::vector<size_t> ip(n);
00295      std::vector<size_t> jp(n);
00296 
00297      // -------------------------------------------------------
00298      CPPAD_ASSERT_KNOWN(
00299           size_t(A.size()) == n * n,
00300           "Error in LuSolve: A must have size equal to n * n"
00301      );
00302      CPPAD_ASSERT_KNOWN(
00303           size_t(B.size()) == n * m,
00304           "Error in LuSolve: B must have size equal to n * m"
00305      );
00306      CPPAD_ASSERT_KNOWN(
00307           size_t(X.size()) == n * m,
00308           "Error in LuSolve: X must have size equal to n * m"
00309      );
00310      // -------------------------------------------------------
00311 
00312      // copy A so that it does not change
00313      FloatVector Lu(A);
00314 
00315      // copy B so that it does not change
00316      X = B;
00317 
00318      // Lu factor the matrix A
00319      signdet = LuFactor(ip, jp, Lu);
00320 
00321      // compute the log of the determinant
00322      logdet  = Float(0);
00323      for(p = 0; p < n; p++)
00324      {    // pivot using the max absolute element
00325           pivot   = Lu[ ip[p] * n + jp[p] ];
00326 
00327           // check for determinant equal to zero
00328           if( pivot == zero )
00329           {    // abort the mission
00330                logdet = Float(0);
00331                return   0;
00332           }
00333 
00334           // update the determinant
00335           if( LeqZero ( pivot ) )
00336           {    logdet += log( - pivot );
00337                signdet = - signdet;
00338           }
00339           else logdet += log( pivot );
00340 
00341      }
00342 
00343      // solve the linear equations
00344      LuInvert(ip, jp, Lu, X);
00345 
00346      // return the sign factor for the determinant
00347      return signdet;
00348 }
00349 } // END CppAD namespace 
00350 // END C++
00351 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines