Prev Next Index-> contents reference index search external Up-> CppAD utility LuDetAndSolve LuInvert lu_invert.hpp utility-> ErrorHandler NearEqual speed_test SpeedTest time_test test_boolofvoid NumericType CheckNumericType SimpleVector CheckSimpleVector nan pow_int Poly LuDetAndSolve RombergOne RombergMul Runge45 Rosen34 OdeErrControl OdeGear OdeGearControl CppAD_vector thread_alloc index_sort to_string set_union sparse_rc sparse_rcv LuDetAndSolve-> LuSolve LuFactor LuInvert LuInvert-> lu_invert.cpp lu_invert.hpp lu_invert.hpp Headings

$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
Source: LuInvert

// LuInvert
template <typename SizeVector, typename FloatVector>
void LuInvert(
const SizeVector  &ip,
const SizeVector  &jp,
const FloatVector &LU,
FloatVector       &B )
{     size_t k; // column index in X
size_t p; // index along diagonal in LU
size_t i; // row index in LU and X

typedef typename FloatVector::value_type Float;

// check numeric type specifications
CheckNumericType<Float>();

// check simple vector class specifications
CheckSimpleVector<Float, FloatVector>();
CheckSimpleVector<size_t, SizeVector>();

Float etmp;

size_t n = ip.size();
size_t(jp.size()) == n,
"Error in LuInvert: jp must have size equal to n * n"
);
size_t(LU.size()) == n * n,
"Error in LuInvert: Lu must have size equal to n * m"
);
size_t m = size_t(B.size()) / n;
size_t(B.size()) == n * m,
"Error in LuSolve: B must have size equal to a multiple of n"
);

// temporary storage for reordered solution
FloatVector x(n);

// loop over equations
for(k = 0; k < m; k++)
{     // invert the equation c = L * b
for(p = 0; p < n; p++)
{     // solve for c[p]
etmp = B[ ip[p] * m + k ] / LU[ ip[p] * n + jp[p] ];
B[ ip[p] * m + k ] = etmp;
// subtract off effect on other variables
for(i = p+1; i < n; i++)
B[ ip[i] * m + k ] -=
etmp * LU[ ip[i] * n + jp[p] ];
}

// invert the equation x = U * c
p = n;
while( p > 0 )
{     --p;
etmp       = B[ ip[p] * m + k ];
x[ jp[p] ] = etmp;
for(i = 0; i < p; i++ )
B[ ip[i] * m + k ] -=
etmp * LU[ ip[i] * n + jp[p] ];
}

// copy reordered solution into B
for(i = 0; i < n; i++)
B[i * m + k] = x[i];
}
return;
}