<?xml version='1.0'?>
<?xml-stylesheet type='text/xsl' href='pmathml.xsl'?>
<html xmlns='http://www.w3.org/1999/xhtml'>
<head>
<title>Source: LuFactor</title>
<meta name="description" id="description" content="Source: LuFactor"/>
<meta name="keywords" id="keywords" content=" Lufactor source "/>
<style type='text/css'>
body { color : black }
body { background-color : white }
A:link { color : blue }
A:visited { color : purple }
A:active { color : purple }
</style>
<script type='text/javascript' language='JavaScript' src='_lu_factor.hpp_xml.js'>
</script>
</head>
<body>
<table><tr>
<td>
<a href="http://www.coin-or.org/CppAD/" target="_top"><img border="0" src="_image.gif"/></a>
</td>
<td><a href="lufactor.cpp.xml" target="_top">Prev</a>
</td><td><a href="luinvert.xml" target="_top">Next</a>
</td><td>
<select onchange='choose_across0(this)'>
<option>Index-></option>
<option>contents</option>
<option>reference</option>
<option>index</option>
<option>search</option>
<option>external</option>
</select>
</td>
<td>
<select onchange='choose_up0(this)'>
<option>Up-></option>
<option>CppAD</option>
<option>library</option>
<option>LuDetAndSolve</option>
<option>LuFactor</option>
<option>lu_factor.hpp</option>
</select>
</td>
<td>
<select onchange='choose_down3(this)'>
<option>library-></option>
<option>ErrorHandler</option>
<option>NearEqual</option>
<option>speed_test</option>
<option>SpeedTest</option>
<option>NumericType</option>
<option>CheckNumericType</option>
<option>SimpleVector</option>
<option>CheckSimpleVector</option>
<option>nan</option>
<option>pow_int</option>
<option>Poly</option>
<option>LuDetAndSolve</option>
<option>RombergOne</option>
<option>RombergMul</option>
<option>Runge45</option>
<option>Rosen34</option>
<option>OdeErrControl</option>
<option>OdeGear</option>
<option>OdeGearControl</option>
<option>BenderQuad</option>
<option>LuRatio</option>
<option>std_math_unary</option>
<option>CppAD_vector</option>
<option>TrackNewDel</option>
</select>
</td>
<td>
<select onchange='choose_down2(this)'>
<option>LuDetAndSolve-></option>
<option>LuSolve</option>
<option>LuFactor</option>
<option>LuInvert</option>
</select>
</td>
<td>
<select onchange='choose_down1(this)'>
<option>LuFactor-></option>
<option>LuFactor.cpp</option>
<option>lu_factor.hpp</option>
</select>
</td>
<td>lu_factor.hpp</td>
<td>Headings</td>
</tr></table><br/>
<center><b><big><big>Source: LuFactor</big></big></b></center>
<code><font color="blue"># ifndef CPPAD_LU_FACTOR_INCLUDED
<code><span style='white-space: nowrap'><br/>
</span></code># define CPPAD_LU_FACTOR_INCLUDED
<pre style='display:inline'>
# include <complex>
# include <vector>
# include <cppad/local/cppad_assert.hpp>
# include <cppad/check_simple_vector.hpp>
# include <cppad/check_numeric_type.hpp>
namespace CppAD { // BEGIN CppAD namespace
// AbsGeq
template <typename Float>
inline bool AbsGeq(const Float &x, const Float &y)
{ Float xabs = x;
if( xabs <= Float(0) )
xabs = - xabs;
Float yabs = y;
if( yabs <= Float(0) )
yabs = - yabs;
return xabs >= yabs;
}
inline bool AbsGeq(
const std::complex<double> &x,
const std::complex<double> &y)
{ double xsq = x.real() * x.real() + x.imag() * x.imag();
double ysq = y.real() * y.real() + y.imag() * y.imag();
return xsq >= ysq;
}
inline bool AbsGeq(
const std::complex<float> &x,
const std::complex<float> &y)
{ float xsq = x.real() * x.real() + x.imag() * x.imag();
float ysq = y.real() * y.real() + y.imag() * y.imag();
return xsq >= ysq;
}
// Lines that are different from code in cppad/local/lu_ratio.hpp end with //
template <class SizeVector, class FloatVector> //
int LuFactor(SizeVector &ip, SizeVector &jp, FloatVector &LU) //
{
// type of the elements of LU //
typedef typename FloatVector::value_type Float; //
// check numeric type specifications
CheckNumericType<Float>();
// check simple vector class specifications
CheckSimpleVector<Float, FloatVector>();
CheckSimpleVector<size_t, SizeVector>();
size_t i, j; // some temporary indices
const Float zero( 0 ); // the value zero as a Float object
size_t imax; // row index of maximum element
size_t jmax; // column indx of maximum element
Float emax; // maximum absolute value
size_t p; // count pivots
int sign; // sign of the permutation
Float etmp; // temporary element
Float pivot; // pivot element
// -------------------------------------------------------
size_t n = ip.size();
CPPAD_ASSERT_KNOWN(
jp.size() == n,
"Error in LuFactor: jp must have size equal to n"
);
CPPAD_ASSERT_KNOWN(
LU.size() == n * n,
"Error in LuFactor: LU must have size equal to n * m"
);
// -------------------------------------------------------
// initialize row and column order in matrix not yet pivoted
for(i = 0; i < n; i++)
{ ip[i] = i;
jp[i] = i;
}
// initialize the sign of the permutation
sign = 1;
// ---------------------------------------------------------
// Reduce the matrix P to L * U using n pivots
for(p = 0; p < n; p++)
{ // determine row and column corresponding to element of
// maximum absolute value in remaining part of P
imax = jmax = n;
emax = zero;
for(i = p; i < n; i++)
{ for(j = p; j < n; j++)
{ CPPAD_ASSERT_UNKNOWN(
(ip[i] < n) & (jp[j] < n)
);
etmp = LU[ ip[i] * n + jp[j] ];
// check if maximum absolute value so far
if( AbsGeq (etmp, emax) )
{ imax = i;
jmax = j;
emax = etmp;
}
}
}
CPPAD_ASSERT_KNOWN(
(imax < n) & (jmax < n) ,
"LuFactor can't determine an element with "
"maximum absolute value.\n"
"Perhaps original matrix contains not a number or infinity.\n"
"Perhaps your specialization of AbsGeq is not correct."
);
if( imax != p )
{ // switch rows so max absolute element is in row p
i = ip[p];
ip[p] = ip[imax];
ip[imax] = i;
sign = -sign;
}
if( jmax != p )
{ // switch columns so max absolute element is in column p
j = jp[p];
jp[p] = jp[jmax];
jp[jmax] = j;
sign = -sign;
}
// pivot using the max absolute element
pivot = LU[ ip[p] * n + jp[p] ];
// check for determinant equal to zero
if( pivot == zero )
{ // abort the mission
return 0;
}
// Reduce U by the elementary transformations that maps
// LU( ip[p], jp[p] ) to one. Only need transform elements
// above the diagonal in U and LU( ip[p] , jp[p] ) is
// corresponding value below diagonal in L.
for(j = p+1; j < n; j++)
LU[ ip[p] * n + jp[j] ] /= pivot;
// Reduce U by the elementary transformations that maps
// LU( ip[i], jp[p] ) to zero. Only need transform elements
// above the diagonal in U and LU( ip[i], jp[p] ) is
// corresponding value below diagonal in L.
for(i = p+1; i < n; i++ )
{ etmp = LU[ ip[i] * n + jp[p] ];
for(j = p+1; j < n; j++)
{ LU[ ip[i] * n + jp[j] ] -=
etmp * LU[ ip[p] * n + jp[j] ];
}
}
}
return sign;
}
} // END CppAD namespace </pre>
# endif
</font></code>
<hr/>Input File: omh/lu_factor_hpp.omh
</body>
</html>
