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