1 # ifndef CPPAD_UTILITY_LU_FACTOR_HPP
2 # define CPPAD_UTILITY_LU_FACTOR_HPP
247 template <
typename Float>
248 inline bool AbsGeq(
const Float &x,
const Float &y)
250 if( xabs <= Float(0) )
253 if( yabs <= Float(0) )
258 const std::complex<double> &x,
259 const std::complex<double> &y)
260 {
double xsq = x.real() * x.real() + x.imag() * x.imag();
261 double ysq = y.real() * y.real() + y.imag() * y.imag();
266 const std::complex<float> &x,
267 const std::complex<float> &y)
268 {
float xsq = x.real() * x.real() + x.imag() * x.imag();
269 float ysq = y.real() * y.real() + y.imag() * y.imag();
275 template <
class SizeVector,
class FloatVector>
276 int LuFactor(SizeVector &ip, SizeVector &jp, FloatVector &LU)
282 CheckNumericType<Float>();
285 CheckSimpleVector<Float, FloatVector>();
286 CheckSimpleVector<size_t, SizeVector>();
289 const Float zero( 0 );
299 size_t n = ip.size();
301 size_t(jp.size()) == n,
302 "Error in LuFactor: jp must have size equal to n"
305 size_t(LU.size()) == n * n,
306 "Error in LuFactor: LU must have size equal to n * m"
311 for(i = 0; i < n; i++)
320 for(p = 0; p < n; p++)
325 for(i = p; i < n; i++)
326 {
for(j = p; j < n; j++)
328 (ip[i] < n) & (jp[j] < n)
330 etmp = LU[ ip[i] * n + jp[j] ];
341 (imax < n) & (jmax < n) ,
342 "LuFactor can't determine an element with "
343 "maximum absolute value.\n"
344 "Perhaps original matrix contains not a number or infinity.\n"
345 "Perhaps your specialization of AbsGeq is not correct."
362 pivot = LU[ ip[p] * n + jp[p] ];
374 for(j = p+1; j < n; j++)
375 LU[ ip[p] * n + jp[j] ] /= pivot;
381 for(i = p+1; i < n; i++ )
382 { etmp = LU[ ip[i] * n + jp[p] ];
383 for(j = p+1; j < n; j++)
384 { LU[ ip[i] * n + jp[j] ] -=
385 etmp * LU[ ip[p] * n + jp[j] ];
#define CPPAD_ASSERT_KNOWN(exp, msg)
Check that exp is true, if not print msg and terminate execution.
bool AbsGeq(const Float &x, const Float &y)
Define the CppAD error checking macros (all of which begin with CPPAD_ASSERT_)
int LuFactor(SizeVector &ip, SizeVector &jp, FloatVector &LU)
#define CPPAD_ASSERT_UNKNOWN(exp)
Check that exp is true, if not terminate execution.
std::complex< double > sign(const std::complex< double > &x)