00001 # ifndef CPPAD_DET_BY_LU_INCLUDED 00002 # define CPPAD_DET_BY_LU_INCLUDED 00003 00004 /* -------------------------------------------------------------------------- 00005 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell 00006 00007 CppAD is distributed under multiple licenses. This distribution is under 00008 the terms of the 00009 Common Public License Version 1.0. 00010 00011 A copy of this license is included in the COPYING file of this distribution. 00012 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00013 -------------------------------------------------------------------------- */ 00014 /* 00015 $begin det_by_lu$$ 00016 $spell 00017 cppad 00018 lu 00019 hpp 00020 typedef 00021 const 00022 hpp 00023 Det 00024 CppADvector 00025 namespace 00026 $$ 00027 00028 $section Determinant Using Expansion by Lu Factorization$$ 00029 00030 $index det_by_lu$$ 00031 $index determinant, lu factor$$ 00032 $index lu, factor determinant$$ 00033 $index factor, lu determinant$$ 00034 00035 $head Syntax$$ 00036 $syntax%# include <cppad/speed/det_by_lu.hpp> 00037 %$$ 00038 $syntax%det_by_lu<%Scalar%> %det%(%n%) 00039 %$$ 00040 $syntax%%d% = %det%(%matrix%) 00041 %$$ 00042 00043 $head Inclusion$$ 00044 The template class $code det_by_lu$$ is defined in the $code CppAD$$ 00045 namespace by including 00046 the file $code cppad/speed/det_by_lu.hpp$$ 00047 (relative to the CppAD distribution directory). 00048 It is only intended for example and testing purposes, 00049 so it is not automatically included by 00050 $cref/cppad.hpp/cppad/$$. 00051 00052 $head Constructor$$ 00053 The syntax 00054 $syntax% 00055 det_by_lu<%Scalar%> %det%(%n%) 00056 %$$ 00057 constructs the object $italic det$$ which can be used for 00058 evaluating the determinant of $italic n$$ by $italic n$$ matrices 00059 using LU factorization. 00060 00061 $head Scalar$$ 00062 The type $italic Scalar$$ can be any 00063 $cref/NumericType/$$ 00064 00065 $head n$$ 00066 The argument $italic n$$ has prototype 00067 $syntax% 00068 size_t %n% 00069 %$$ 00070 00071 $head det$$ 00072 The syntax 00073 $syntax% 00074 %d% = %det%(%matrix%) 00075 %$$ 00076 returns the determinant of $italic matrix$$ using LU factorization. 00077 The argument $italic matrix$$ has prototype 00078 $syntax% 00079 const %Vector% &%matrix% 00080 %$$ 00081 It must be a $italic Vector$$ with length $latex n * n$$ and with 00082 elements of type $italic Scalar$$. 00083 The return value $italic d$$ has prototype 00084 $syntax% 00085 %Scalar% %d% 00086 %$$ 00087 00088 $head Vector$$ 00089 If $italic y$$ is a $italic Vector$$ object, 00090 it must support the syntax 00091 $syntax% 00092 %y%[%i%] 00093 %$$ 00094 where $italic i$$ has type $code size_t$$ with value less than $latex n * n$$. 00095 This must return a $italic Scalar$$ value corresponding to the $th i$$ 00096 element of the vector $italic y$$. 00097 This is the only requirement of the type $italic Vector$$. 00098 00099 $children% 00100 speed/example/det_by_lu.cpp% 00101 omh/det_by_lu_hpp.omh 00102 %$$ 00103 00104 00105 $head Example$$ 00106 The file 00107 $xref/det_by_lu.cpp/$$ 00108 contains an example and test of $code det_by_lu.hpp$$. 00109 It returns true if it succeeds and false otherwise. 00110 00111 $head Source Code$$ 00112 The file 00113 $cref/det_by_lu.hpp/$$ 00114 contains the source for this template function. 00115 00116 00117 $end 00118 --------------------------------------------------------------------------- 00119 */ 00120 // BEGIN PROGRAM 00121 # include <cppad/cppad.hpp> 00122 # include <complex> 00123 00124 // BEGIN CppAD namespace 00125 namespace CppAD { 00126 00127 // The AD complex case is used by examples by not used by speed tests 00128 // Must define a specializatgion of LeqZero,AbsGeq for the ADComplex case 00129 typedef std::complex<double> Complex; 00130 typedef CppAD::AD<Complex> ADComplex; 00131 CPPAD_BOOL_UNARY(Complex, LeqZero ) 00132 CPPAD_BOOL_BINARY(Complex, AbsGeq ) 00133 00134 template <class Scalar> 00135 class det_by_lu { 00136 public: 00137 det_by_lu(size_t n_) : m(0), n(n_), A(n_ * n_) 00138 { } 00139 00140 template <class Vector> 00141 inline Scalar operator()(const Vector &x) 00142 { 00143 using CppAD::exp; 00144 00145 Scalar logdet; 00146 Scalar det; 00147 int signdet; 00148 size_t i; 00149 00150 // copy matrix so it is not overwritten 00151 for(i = 0; i < n * n; i++) 00152 A[i] = x[i]; 00153 00154 // comput log determinant 00155 signdet = CppAD::LuSolve( 00156 n, m, A, B, X, logdet); 00157 00158 // make sure the martix is not singular 00159 CPPAD_ASSERT_KNOWN( 00160 signdet != 0, 00161 "det_by_lu: matrix is singular" 00162 ); 00163 00164 // convert to determinant 00165 det = Scalar( signdet ) * exp( logdet ); 00166 00167 # ifdef FADBAD 00168 // Fadbad requires tempories to be set to constants 00169 for(i = 0; i < n * n; i++) 00170 A[i] = 0; 00171 # endif 00172 00173 return det; 00174 } 00175 private: 00176 const size_t m; 00177 const size_t n; 00178 CppADvector<Scalar> A; 00179 CppADvector<Scalar> B; 00180 CppADvector<Scalar> X; 00181 }; 00182 } // END CppAD namespace 00183 // END PROGRAM 00184 # endif