CppAD: A C++ Algorithmic Differentiation Package  20130102
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-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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines