CppAD: A C++ Algorithmic Differentiation Package 20110419
lu_invert.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_LU_INVERT_INCLUDED
00003 # define CPPAD_LU_INVERT_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-06 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 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 $syntax%LuInvert(%ip%, %jp%, %LU%, %X%)%$$
00044 
00045 
00046 $head Description$$
00047 Solves the matrix equation $syntax%%A% * %X% = %B%$$ 
00048 using an LU factorization computed by $xref/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 $italic ip$$ has prototype
00068 $syntax%
00069         const %SizeVector% &%ip%
00070 %$$
00071 (see description for $italic SizeVector$$ in
00072 $xref/LuFactor/SizeVector/LuFactor/$$ specifications).
00073 The size of $italic ip$$ is referred to as $italic n$$ in the
00074 specifications below.
00075 The elements of $italic ip$$ determine
00076 the order of the rows in the permuted matrix.
00077 
00078 $head jp$$
00079 The argument $italic jp$$ has prototype
00080 $syntax%
00081         const %SizeVector% &%jp%
00082 %$$
00083 (see description for $italic SizeVector$$ in
00084 $xref/LuFactor/SizeVector/LuFactor/$$ specifications).
00085 The size of $italic jp$$ must be equal to $italic n$$.
00086 The elements of $italic jp$$ determine
00087 the order of the columns in the permuted matrix.
00088 
00089 $head LU$$
00090 The argument $italic LU$$ has the prototype
00091 $syntax%
00092         const %FloatVector% &%LU%
00093 %$$
00094 and the size of $italic LU$$ must equal $latex n * n$$
00095 (see description for $italic FloatVector$$ in
00096 $xref/LuFactor/FloatVector/LuFactor/$$ specifications).
00097 
00098 $subhead L$$
00099 We define the lower triangular matrix $italic L$$ in terms of $italic LU$$.
00100 The matrix $italic L$$ is zero above the diagonal
00101 and the rest of the elements are defined by
00102 $syntax%
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 $italic U$$ in terms of $italic LU$$.
00109 The matrix $italic U$$ is zero below the diagonal,
00110 one on the diagonal,
00111 and the rest of the elements are defined by
00112 $syntax%
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 $italic P$$ in terms of 
00119 the matrix $italic L$$ and the matrix $italic U$$ 
00120 by $syntax%%P% = %L% * %U%$$.
00121 
00122 $subhead A$$
00123 The matrix $italic A$$, 
00124 which defines the linear equations that we are solving, is given by
00125 $syntax%
00126         %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ]
00127 %$$
00128 (Hence 
00129 $italic LU$$ contains a permuted factorization of the matrix $italic A$$.)
00130 
00131 
00132 $head X$$
00133 The argument $italic X$$ has prototype
00134 $syntax%
00135         %FloatVector% &%X%
00136 %$$
00137 (see description for $italic FloatVector$$ in
00138 $xref/LuFactor/FloatVector/LuFactor/$$ specifications).
00139 The matrix $italic X$$
00140 must have the same number of rows as the matrix $italic A$$.
00141 The input value of $italic X$$ is the matrix $italic B$$ and the 
00142 output value solves the matrix equation $syntax%%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 $xref/lu_solve.hpp/$$ is a good example usage of 
00151 $code LuFactor$$ with $code LuInvert$$.
00152 The file 
00153 $xref/LuInvert.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 PROGRAM
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                 jp.size() == n,
00196                 "Error in LuInvert: jp must have size equal to n * n"
00197         );
00198         CPPAD_ASSERT_KNOWN(
00199                 LU.size() == n * n,
00200                 "Error in LuInvert: Lu must have size equal to n * m"
00201         );
00202         size_t m = B.size() / n;
00203         CPPAD_ASSERT_KNOWN(
00204                 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 PROGRAM
00243 # endif