CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_LU_INVERT_INCLUDED 00003 # define CPPAD_LU_INVERT_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 LuInvert$$ 00018 $escape #$$ 00019 $spell 00020 cppad.hpp 00021 Lu 00022 Cpp 00023 jp 00024 ip 00025 const 00026 namespace 00027 typename 00028 etmp 00029 $$ 00030 00031 $index LuInvert$$ 00032 $index linear, invert Lu equation$$ 00033 $index equation, Lu invert$$ 00034 00035 $section Invert an LU Factored Equation$$ 00036 00037 $pre 00038 $$ 00039 00040 $head Syntax$$ $code# include <cppad/lu_invert.hpp>$$ 00041 $pre 00042 $$ 00043 $codei%LuInvert(%ip%, %jp%, %LU%, %X%)%$$ 00044 00045 00046 $head Description$$ 00047 Solves the matrix equation $icode%A% * %X% = %B%$$ 00048 using an LU factorization computed by $cref LuFactor$$. 00049 00050 $head Include$$ 00051 The file $code cppad/lu_invert.hpp$$ is included by $code cppad/cppad.hpp$$ 00052 but it can also be included separately with out the rest of 00053 the $code CppAD$$ routines. 00054 00055 $head Matrix Storage$$ 00056 All matrices are stored in row major order. 00057 To be specific, if $latex Y$$ is a vector 00058 that contains a $latex p$$ by $latex q$$ matrix, 00059 the size of $latex Y$$ must be equal to $latex p * q $$ and for 00060 $latex i = 0 , \ldots , p-1$$, 00061 $latex j = 0 , \ldots , q-1$$, 00062 $latex \[ 00063 Y_{i,j} = Y[ i * q + j ] 00064 \] $$ 00065 00066 $head ip$$ 00067 The argument $icode ip$$ has prototype 00068 $codei% 00069 const %SizeVector% &%ip% 00070 %$$ 00071 (see description for $icode SizeVector$$ in 00072 $cref/LuFactor/LuFactor/SizeVector/$$ specifications). 00073 The size of $icode ip$$ is referred to as $icode n$$ in the 00074 specifications below. 00075 The elements of $icode ip$$ determine 00076 the order of the rows in the permuted matrix. 00077 00078 $head jp$$ 00079 The argument $icode jp$$ has prototype 00080 $codei% 00081 const %SizeVector% &%jp% 00082 %$$ 00083 (see description for $icode SizeVector$$ in 00084 $cref/LuFactor/LuFactor/SizeVector/$$ specifications). 00085 The size of $icode jp$$ must be equal to $icode n$$. 00086 The elements of $icode jp$$ determine 00087 the order of the columns in the permuted matrix. 00088 00089 $head LU$$ 00090 The argument $icode LU$$ has the prototype 00091 $codei% 00092 const %FloatVector% &%LU% 00093 %$$ 00094 and the size of $icode LU$$ must equal $latex n * n$$ 00095 (see description for $icode FloatVector$$ in 00096 $cref/LuFactor/LuFactor/FloatVector/$$ specifications). 00097 00098 $subhead L$$ 00099 We define the lower triangular matrix $icode L$$ in terms of $icode LU$$. 00100 The matrix $icode L$$ is zero above the diagonal 00101 and the rest of the elements are defined by 00102 $codei% 00103 %L%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ] 00104 %$$ 00105 for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , i$$. 00106 00107 $subhead U$$ 00108 We define the upper triangular matrix $icode U$$ in terms of $icode LU$$. 00109 The matrix $icode U$$ is zero below the diagonal, 00110 one on the diagonal, 00111 and the rest of the elements are defined by 00112 $codei% 00113 %U%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ] 00114 %$$ 00115 for $latex i = 0 , \ldots , n-2$$ and $latex j = i+1 , \ldots , n-1$$. 00116 00117 $subhead P$$ 00118 We define the permuted matrix $icode P$$ in terms of 00119 the matrix $icode L$$ and the matrix $icode U$$ 00120 by $icode%P% = %L% * %U%$$. 00121 00122 $subhead A$$ 00123 The matrix $icode A$$, 00124 which defines the linear equations that we are solving, is given by 00125 $codei% 00126 %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ] 00127 %$$ 00128 (Hence 00129 $icode LU$$ contains a permuted factorization of the matrix $icode A$$.) 00130 00131 00132 $head X$$ 00133 The argument $icode X$$ has prototype 00134 $codei% 00135 %FloatVector% &%X% 00136 %$$ 00137 (see description for $icode FloatVector$$ in 00138 $cref/LuFactor/LuFactor/FloatVector/$$ specifications). 00139 The matrix $icode X$$ 00140 must have the same number of rows as the matrix $icode A$$. 00141 The input value of $icode X$$ is the matrix $icode B$$ and the 00142 output value solves the matrix equation $icode%A% * %X% = %B%$$. 00143 00144 00145 $children% 00146 example/lu_invert.cpp% 00147 omh/lu_invert_hpp.omh 00148 %$$ 00149 $head Example$$ 00150 The file $cref lu_solve.hpp$$ is a good example usage of 00151 $code LuFactor$$ with $code LuInvert$$. 00152 The file 00153 $cref lu_invert.cpp$$ 00154 contains an example and test of using $code LuInvert$$ by itself. 00155 It returns true if it succeeds and false otherwise. 00156 00157 $head Source$$ 00158 The file $cref lu_invert.hpp$$ contains the 00159 current source code that implements these specifications. 00160 00161 $end 00162 -------------------------------------------------------------------------- 00163 */ 00164 // BEGIN C++ 00165 # include <cppad/local/cppad_assert.hpp> 00166 # include <cppad/check_simple_vector.hpp> 00167 # include <cppad/check_numeric_type.hpp> 00168 00169 namespace CppAD { // BEGIN CppAD namespace 00170 00171 // LuInvert 00172 template <typename SizeVector, typename FloatVector> 00173 void LuInvert( 00174 const SizeVector &ip, 00175 const SizeVector &jp, 00176 const FloatVector &LU, 00177 FloatVector &B ) 00178 { size_t k; // column index in X 00179 size_t p; // index along diagonal in LU 00180 size_t i; // row index in LU and X 00181 00182 typedef typename FloatVector::value_type Float; 00183 00184 // check numeric type specifications 00185 CheckNumericType<Float>(); 00186 00187 // check simple vector class specifications 00188 CheckSimpleVector<Float, FloatVector>(); 00189 CheckSimpleVector<size_t, SizeVector>(); 00190 00191 Float etmp; 00192 00193 size_t n = ip.size(); 00194 CPPAD_ASSERT_KNOWN( 00195 size_t(jp.size()) == n, 00196 "Error in LuInvert: jp must have size equal to n * n" 00197 ); 00198 CPPAD_ASSERT_KNOWN( 00199 size_t(LU.size()) == n * n, 00200 "Error in LuInvert: Lu must have size equal to n * m" 00201 ); 00202 size_t m = size_t(B.size()) / n; 00203 CPPAD_ASSERT_KNOWN( 00204 size_t(B.size()) == n * m, 00205 "Error in LuSolve: B must have size equal to a multiple of n" 00206 ); 00207 00208 // temporary storage for reordered solution 00209 FloatVector x(n); 00210 00211 // loop over equations 00212 for(k = 0; k < m; k++) 00213 { // invert the equation c = L * b 00214 for(p = 0; p < n; p++) 00215 { // solve for c[p] 00216 etmp = B[ ip[p] * m + k ] / LU[ ip[p] * n + jp[p] ]; 00217 B[ ip[p] * m + k ] = etmp; 00218 // subtract off effect on other variables 00219 for(i = p+1; i < n; i++) 00220 B[ ip[i] * m + k ] -= 00221 etmp * LU[ ip[i] * n + jp[p] ]; 00222 } 00223 00224 // invert the equation x = U * c 00225 p = n; 00226 while( p > 0 ) 00227 { --p; 00228 etmp = B[ ip[p] * m + k ]; 00229 x[ jp[p] ] = etmp; 00230 for(i = 0; i < p; i++ ) 00231 B[ ip[i] * m + k ] -= 00232 etmp * LU[ ip[i] * n + jp[p] ]; 00233 } 00234 00235 // copy reordered solution into B 00236 for(i = 0; i < n; i++) 00237 B[i * m + k] = x[i]; 00238 } 00239 return; 00240 } 00241 } // END CppAD namespace 00242 // END C++ 00243 # endif