<?xml version='1.0'?> <?xml-stylesheet type='text/xsl' href='pmathml.xsl'?> <html xmlns='http://www.w3.org/1999/xhtml'> <head> <title>LU Factorization of A Square Matrix</title> <meta name="description" id="description" content="LU Factorization of A Square Matrix"/> <meta name="keywords" id="keywords" content=" Lufactor linear Lu factor equation determinant solve "/> <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='_lufactor_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="lu_solve.hpp.xml" target="_top">Prev</a> </td><td><a href="lufactor.cpp.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> </select> </td> <td> <select onchange='choose_down3(this)'> <option>CppAD-></option> <option>Install</option> <option>Introduction</option> <option>AD</option> <option>ADFun</option> <option>library</option> <option>Example</option> <option>configure</option> <option>Appendix</option> </select> </td> <td> <select onchange='choose_down2(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_down1(this)'> <option>LuDetAndSolve-></option> <option>LuSolve</option> <option>LuFactor</option> <option>LuInvert</option> </select> </td> <td> <select onchange='choose_down0(this)'> <option>LuFactor-></option> <option>LuFactor.cpp</option> <option>lu_factor.hpp</option> </select> </td> <td> <select onchange='choose_current0(this)'> <option>Headings-></option> <option>Syntax</option> <option>Description</option> <option>Include</option> <option>Matrix Storage</option> <option>sign</option> <option>ip</option> <option>jp</option> <option>LU</option> <option>---..A</option> <option>---..P</option> <option>---..L</option> <option>---..U</option> <option>---..Factor</option> <option>---..Determinant</option> <option>SizeVector</option> <option>FloatVector</option> <option>Float</option> <option>AbsGeq</option> <option>Example</option> <option>Source</option> </select> </td> </tr></table><br/> <center><b><big><big>LU Factorization of A Square Matrix</big></big></b></center> <code><span style='white-space: nowrap'><br/> </span></code><b><big><a name="Syntax" id="Syntax">Syntax</a></big></b> <code><font color="blue"><br/> # include <cppad/lu_factor.hpp></font></code> <code><span style='white-space: nowrap'><br/> </span></code><code><font color="blue"></font></code><i><span style='white-space: nowrap'>sign</span></i><code><font color="blue"><span style='white-space: nowrap'> = LuFactor(</span></font></code><i><span style='white-space: nowrap'>ip</span></i><code><font color="blue"><span style='white-space: nowrap'>, </span></font></code><i><span style='white-space: nowrap'>jp</span></i><code><font color="blue"><span style='white-space: nowrap'>, </span></font></code><i><span style='white-space: nowrap'>LU</span></i><code><font color="blue"><span style='white-space: nowrap'>)</span></font></code> <br/> <br/> <b><big><a name="Description" id="Description">Description</a></big></b> <br/> Computes an LU factorization of the matrix <i>A</i> where <i>A</i> is a square matrix. <br/> <br/> <b><big><a name="Include" id="Include">Include</a></big></b> <br/> The file <code><font color="blue">cppad/lu_factor.hpp</font></code> is included by <code><font color="blue">cppad/cppad.hpp</font></code> but it can also be included separately with out the rest of the <code><font color="blue">CppAD</font></code> routines. <br/> <br/> <b><big><a name="Matrix Storage" id="Matrix Storage">Matrix Storage</a></big></b> <br/> All matrices are stored in row major order. To be specific, if <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>Y</mi> </mrow></math> is a vector that contains a <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>p</mi> </mrow></math> by <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>q</mi> </mrow></math> matrix, the size of <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>Y</mi> </mrow></math> must be equal to <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>p</mi> <mo stretchy="false">*</mo> <mi mathvariant='italic'>q</mi> </mrow></math> and for <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>i</mi> <mo stretchy="false">=</mo> <mn>0</mn> <mo stretchy="false">,</mo> <mo stretchy="false">…</mo> <mo stretchy="false">,</mo> <mi mathvariant='italic'>p</mi> <mn>-1</mn> </mrow></math> , <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>j</mi> <mo stretchy="false">=</mo> <mn>0</mn> <mo stretchy="false">,</mo> <mo stretchy="false">…</mo> <mo stretchy="false">,</mo> <mi mathvariant='italic'>q</mi> <mn>-1</mn> </mrow></math> , <math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><mrow> <msub><mi mathvariant='italic'>Y</mi> <mrow><mi mathvariant='italic'>i</mi> <mo stretchy="false">,</mo> <mi mathvariant='italic'>j</mi> </mrow> </msub> <mo stretchy="false">=</mo> <mi mathvariant='italic'>Y</mi> <mo stretchy="false">[</mo> <mi mathvariant='italic'>i</mi> <mo stretchy="false">*</mo> <mi mathvariant='italic'>q</mi> <mo stretchy="false">+</mo> <mi mathvariant='italic'>j</mi> <mo stretchy="false">]</mo> </mrow></math> <br/> <b><big><a name="sign" id="sign">sign</a></big></b> <br/> The return value <i>sign</i> has prototype <code><font color="blue"><span style='white-space: nowrap'><br/>      int </span></font></code><i><span style='white-space: nowrap'>sign</span></i><code><font color="blue"><span style='white-space: nowrap'><br/> </span></font></code>If <i>A</i> is invertible, <i>sign</i> is plus or minus one and is the sign of the permutation corresponding to the row ordering <i>ip</i> and column ordering <i>jp</i>. If <i>A</i> is not invertible, <i>sign</i> is zero. <br/> <br/> <b><big><a name="ip" id="ip">ip</a></big></b> <br/> The argument <i>ip</i> has prototype <code><font color="blue"><span style='white-space: nowrap'><br/>      </span></font></code><i><span style='white-space: nowrap'>SizeVector</span></i><code><font color="blue"><span style='white-space: nowrap'> &</span></font></code><i><span style='white-space: nowrap'>ip</span></i><code><font color="blue"><span style='white-space: nowrap'><br/> </span></font></code>(see description of <a href="lufactor.xml#SizeVector" target="_top"><span style='white-space: nowrap'>SizeVector</span></a> below). The size of <i>ip</i> is referred to as <i>n</i> in the specifications below. The input value of the elements of <i>ip</i> does not matter. The output value of the elements of <i>ip</i> determine the order of the rows in the permuted matrix. <br/> <br/> <b><big><a name="jp" id="jp">jp</a></big></b> <br/> The argument <i>jp</i> has prototype <code><font color="blue"><span style='white-space: nowrap'><br/>      </span></font></code><i><span style='white-space: nowrap'>SizeVector</span></i><code><font color="blue"><span style='white-space: nowrap'> &</span></font></code><i><span style='white-space: nowrap'>jp</span></i><code><font color="blue"><span style='white-space: nowrap'><br/> </span></font></code>(see description of <a href="lufactor.xml#SizeVector" target="_top"><span style='white-space: nowrap'>SizeVector</span></a> below). The size of <i>jp</i> must be equal to <i>n</i>. The input value of the elements of <i>jp</i> does not matter. The output value of the elements of <i>jp</i> determine the order of the columns in the permuted matrix. <br/> <br/> <b><big><a name="LU" id="LU">LU</a></big></b> <br/> The argument <i>LU</i> has the prototype <code><font color="blue"><span style='white-space: nowrap'><br/>      </span></font></code><i><span style='white-space: nowrap'>FloatVector</span></i><code><font color="blue"><span style='white-space: nowrap'> &</span></font></code><i><span style='white-space: nowrap'>LU</span></i><code><font color="blue"><span style='white-space: nowrap'><br/> </span></font></code>and the size of <i>LU</i> must equal <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>n</mi> <mo stretchy="false">*</mo> <mi mathvariant='italic'>n</mi> </mrow></math> (see description of <a href="lufactor.xml#FloatVector" target="_top"><span style='white-space: nowrap'>FloatVector</span></a> below). <br/> <br/> <b><a name="LU.A" id="LU.A">A</a></b> <br/> We define <i>A</i> as the matrix corresponding to the input value of <i>LU</i>. <br/> <br/> <b><a name="LU.P" id="LU.P">P</a></b> <br/> We define the permuted matrix <i>P</i> in terms of <i>A</i> by <code><font color="blue"><span style='white-space: nowrap'><br/>      </span></font></code><i><span style='white-space: nowrap'>P</span></i><code><font color="blue"><span style='white-space: nowrap'>(</span></font></code><i><span style='white-space: nowrap'>i</span></i><code><font color="blue"><span style='white-space: nowrap'>, </span></font></code><i><span style='white-space: nowrap'>j</span></i><code><font color="blue"><span style='white-space: nowrap'>) = </span></font></code><i><span style='white-space: nowrap'>A</span></i><code><font color="blue"><span style='white-space: nowrap'>[ </span></font></code><i><span style='white-space: nowrap'>ip</span></i><code><font color="blue"><span style='white-space: nowrap'>[</span></font></code><i><span style='white-space: nowrap'>i</span></i><code><font color="blue"><span style='white-space: nowrap'>] * </span></font></code><i><span style='white-space: nowrap'>n</span></i><code><font color="blue"><span style='white-space: nowrap'> + </span></font></code><i><span style='white-space: nowrap'>jp</span></i><code><font color="blue"><span style='white-space: nowrap'>[</span></font></code><i><span style='white-space: nowrap'>j</span></i><code><font color="blue"><span style='white-space: nowrap'>] ]<br/> </span></font></code><br/> <b><a name="LU.L" id="LU.L">L</a></b> <br/> We define the lower triangular matrix <i>L</i> in terms of the output value of <i>LU</i>. The matrix <i>L</i> is zero above the diagonal and the rest of the elements are defined by <code><font color="blue"><span style='white-space: nowrap'><br/>      </span></font></code><i><span style='white-space: nowrap'>L</span></i><code><font color="blue"><span style='white-space: nowrap'>(</span></font></code><i><span style='white-space: nowrap'>i</span></i><code><font color="blue"><span style='white-space: nowrap'>, </span></font></code><i><span style='white-space: nowrap'>j</span></i><code><font color="blue"><span style='white-space: nowrap'>) = </span></font></code><i><span style='white-space: nowrap'>LU</span></i><code><font color="blue"><span style='white-space: nowrap'>[ </span></font></code><i><span style='white-space: nowrap'>ip</span></i><code><font color="blue"><span style='white-space: nowrap'>[</span></font></code><i><span style='white-space: nowrap'>i</span></i><code><font color="blue"><span style='white-space: nowrap'>] * </span></font></code><i><span style='white-space: nowrap'>n</span></i><code><font color="blue"><span style='white-space: nowrap'> + </span></font></code><i><span style='white-space: nowrap'>jp</span></i><code><font color="blue"><span style='white-space: nowrap'>[</span></font></code><i><span style='white-space: nowrap'>j</span></i><code><font color="blue"><span style='white-space: nowrap'>] ]<br/> </span></font></code>for <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>i</mi> <mo stretchy="false">=</mo> <mn>0</mn> <mo stretchy="false">,</mo> <mo stretchy="false">…</mo> <mo stretchy="false">,</mo> <mi mathvariant='italic'>n</mi> <mn>-1</mn> </mrow></math> and <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>j</mi> <mo stretchy="false">=</mo> <mn>0</mn> <mo stretchy="false">,</mo> <mo stretchy="false">…</mo> <mo stretchy="false">,</mo> <mi mathvariant='italic'>i</mi> </mrow></math> . <br/> <br/> <b><a name="LU.U" id="LU.U">U</a></b> <br/> We define the upper triangular matrix <i>U</i> in terms of the output value of <i>LU</i>. The matrix <i>U</i> is zero below the diagonal, one on the diagonal, and the rest of the elements are defined by <code><font color="blue"><span style='white-space: nowrap'><br/>      </span></font></code><i><span style='white-space: nowrap'>U</span></i><code><font color="blue"><span style='white-space: nowrap'>(</span></font></code><i><span style='white-space: nowrap'>i</span></i><code><font color="blue"><span style='white-space: nowrap'>, </span></font></code><i><span style='white-space: nowrap'>j</span></i><code><font color="blue"><span style='white-space: nowrap'>) = </span></font></code><i><span style='white-space: nowrap'>LU</span></i><code><font color="blue"><span style='white-space: nowrap'>[ </span></font></code><i><span style='white-space: nowrap'>ip</span></i><code><font color="blue"><span style='white-space: nowrap'>[</span></font></code><i><span style='white-space: nowrap'>i</span></i><code><font color="blue"><span style='white-space: nowrap'>] * </span></font></code><i><span style='white-space: nowrap'>n</span></i><code><font color="blue"><span style='white-space: nowrap'> + </span></font></code><i><span style='white-space: nowrap'>jp</span></i><code><font color="blue"><span style='white-space: nowrap'>[</span></font></code><i><span style='white-space: nowrap'>j</span></i><code><font color="blue"><span style='white-space: nowrap'>] ]<br/> </span></font></code>for <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>i</mi> <mo stretchy="false">=</mo> <mn>0</mn> <mo stretchy="false">,</mo> <mo stretchy="false">…</mo> <mo stretchy="false">,</mo> <mi mathvariant='italic'>n</mi> <mn>-2</mn> </mrow></math> and <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mrow> <mi mathvariant='italic'>j</mi> <mo stretchy="false">=</mo> <mi mathvariant='italic'>i</mi> <mo stretchy="false">+</mo> <mn>1</mn> <mo stretchy="false">,</mo> <mo stretchy="false">…</mo> <mo stretchy="false">,</mo> <mi mathvariant='italic'>n</mi> <mn>-1</mn> </mrow></math> . <br/> <br/> <b><a name="LU.Factor" id="LU.Factor">Factor</a></b> <br/> If the return value <i>sign</i> is non-zero, <code><font color="blue"><span style='white-space: nowrap'><br/>      </span></font></code><i><span style='white-space: nowrap'>L</span></i><code><font color="blue"><span style='white-space: nowrap'> * </span></font></code><i><span style='white-space: nowrap'>U</span></i><code><font color="blue"><span style='white-space: nowrap'> = </span></font></code><i><span style='white-space: nowrap'>P</span></i><code><font color="blue"><span style='white-space: nowrap'><br/> </span></font></code>If the return value of <i>sign</i> is zero, the contents of <i>L</i> and <i>U</i> are not defined. <br/> <br/> <b><a name="LU.Determinant" id="LU.Determinant">Determinant</a></b> <br/> If the return value <i>sign</i> is zero, the determinant of <i>A</i> is zero. If <i>sign</i> is non-zero, using the output value of <i>LU</i> the determinant of the matrix <i>A</i> is equal to <code><font color="blue"><span style='white-space: nowrap'><br/> </span></font></code><i><span style='white-space: nowrap'>sign</span></i><code><font color="blue"><span style='white-space: nowrap'> * </span></font></code><i><span style='white-space: nowrap'>LU</span></i><code><font color="blue"><span style='white-space: nowrap'>[</span></font></code><i><span style='white-space: nowrap'>ip</span></i><code><font color="blue"><span style='white-space: nowrap'>[0], </span></font></code><i><span style='white-space: nowrap'>jp</span></i><code><font color="blue"><span style='white-space: nowrap'>[0]] * </span></font></code><i><span style='white-space: nowrap'>...</span></i><code><font color="blue"><span style='white-space: nowrap'> * </span></font></code><i><span style='white-space: nowrap'>LU</span></i><code><font color="blue"><span style='white-space: nowrap'>[</span></font></code><i><span style='white-space: nowrap'>ip</span></i><code><font color="blue"><span style='white-space: nowrap'>[</span></font></code><i><span style='white-space: nowrap'>n</span></i><code><font color="blue"><span style='white-space: nowrap'>-1], </span></font></code><i><span style='white-space: nowrap'>jp</span></i><code><font color="blue"><span style='white-space: nowrap'>[</span></font></code><i><span style='white-space: nowrap'>n</span></i><code><font color="blue"><span style='white-space: nowrap'>-1]] <br/> </span></font></code><br/> <b><big><a name="SizeVector" id="SizeVector">SizeVector</a></big></b> <br/> The type <i>SizeVector</i> must be a <a href="simplevector.xml" target="_top"><span style='white-space: nowrap'>SimpleVector</span></a> class with <a href="simplevector.xml#Elements of Specified Type" target="_top"><span style='white-space: nowrap'>elements of type size_t</span></a> . The routine <a href="checksimplevector.xml" target="_top"><span style='white-space: nowrap'>CheckSimpleVector</span></a> will generate an error message if this is not the case. <br/> <br/> <b><big><a name="FloatVector" id="FloatVector">FloatVector</a></big></b> <br/> The type <i>FloatVector</i> must be a <a href="simplevector.xml" target="_top"><span style='white-space: nowrap'>simple vector class</span></a> . The routine <a href="checksimplevector.xml" target="_top"><span style='white-space: nowrap'>CheckSimpleVector</span></a> will generate an error message if this is not the case. <br/> <br/> <b><big><a name="Float" id="Float">Float</a></big></b> <br/> This notation is used to denote the type corresponding to the elements of a <i>FloatVector</i>. The type <i>Float</i> must satisfy the conditions for a <a href="numerictype.xml" target="_top"><span style='white-space: nowrap'>NumericType</span></a> type. The routine <a href="checknumerictype.xml" target="_top"><span style='white-space: nowrap'>CheckNumericType</span></a> will generate an error message if this is not the case. In addition, the following operations must be defined for any pair of <i>Float</i> objects <i>x</i> and <i>y</i>: <table><tr><td align='left' valign='top'> <b>Operation</b> </td><td align='left' valign='top'> <b>Description</b> </td></tr><tr><td align='left' valign='top'> <code><font color="blue"><span style='white-space: nowrap'>log(</span></font></code><i><span style='white-space: nowrap'>x</span></i><code><font color="blue"><span style='white-space: nowrap'>)</span></font></code> </td><td align='left' valign='top'> returns the logarithm of <i>x</i> as a <i>Float</i> object </td></tr> </table> <br/> <b><big><a name="AbsGeq" id="AbsGeq">AbsGeq</a></big></b> <br/> Including the file <code><font color="blue">lu_factor.hpp</font></code> defines the template function <code><font color="blue"><span style='white-space: nowrap'><br/>      template <typename </span></font></code><i><span style='white-space: nowrap'>Float</span></i><code><font color="blue"><span style='white-space: nowrap'>><br/>      bool AbsGeq<</span></font></code><i><span style='white-space: nowrap'>Float</span></i><code><font color="blue"><span style='white-space: nowrap'>>(const </span></font></code><i><span style='white-space: nowrap'>Float</span></i><code><font color="blue"><span style='white-space: nowrap'> &</span></font></code><i><span style='white-space: nowrap'>x</span></i><code><font color="blue"><span style='white-space: nowrap'>, const </span></font></code><i><span style='white-space: nowrap'>Float</span></i><code><font color="blue"><span style='white-space: nowrap'> &</span></font></code><i><span style='white-space: nowrap'>y</span></i><code><font color="blue"><span style='white-space: nowrap'>)<br/> </span></font></code>in the <code><font color="blue">CppAD</font></code> namespace. This function returns true if the absolute value of <i>x</i> is greater than or equal the absolute value of <i>y</i>. It is used by <code><font color="blue">LuFactor</font></code> to choose the pivot elements. This template function definition uses the operator <code><font color="blue"><=</font></code> to obtain the absolute value for <i>Float</i> objects. If this operator is not defined for your use of <i>Float</i>, you will need to specialize this template so that it works for your use of <code><font color="blue">LuFactor</font></code>. <code><span style='white-space: nowrap'><br/> <br/> </span></code>Complex numbers do not have the operation <code><font color="blue"><=</font></code> defined. The specializations <code><font color="blue"><span style='white-space: nowrap'><br/> bool AbsGeq< std::complex<float> > <br/>      (const std::complex<float> &</span></font></code><i><span style='white-space: nowrap'>x</span></i><code><font color="blue"><span style='white-space: nowrap'>, const std::complex<float> &</span></font></code><i><span style='white-space: nowrap'>y</span></i><code><font color="blue"><span style='white-space: nowrap'>)<br/> bool AbsGeq< std::complex<double> > <br/>      (const std::complex<double> &</span></font></code><i><span style='white-space: nowrap'>x</span></i><code><font color="blue"><span style='white-space: nowrap'>, const std::complex<double> &</span></font></code><i><span style='white-space: nowrap'>y</span></i><code><font color="blue"><span style='white-space: nowrap'>)<br/> </span></font></code>are define by including <code><font color="blue">lu_factor.hpp</font></code> These return true if the sum of the square of the real and imaginary parts of <i>x</i> is greater than or equal the sum of the square of the real and imaginary parts of <i>y</i>. <br/> <br/> <b><big><a name="Example" id="Example">Example</a></big></b> <br/> The file <a href="lufactor.cpp.xml" target="_top"><span style='white-space: nowrap'>LuFactor.cpp</span></a> contains an example and test of using <code><font color="blue">LuFactor</font></code> by itself. It returns true if it succeeds and false otherwise. <code><span style='white-space: nowrap'><br/> <br/> </span></code>The file <a href="lu_solve.hpp.xml" target="_top"><span style='white-space: nowrap'>lu_solve.hpp</span></a> provides a useful example usage of <code><font color="blue">LuFactor</font></code> with <code><font color="blue">LuInvert</font></code>. <br/> <br/> <b><big><a name="Source" id="Source">Source</a></big></b> <br/> The file <a href="lu_factor.hpp.xml" target="_top"><span style='white-space: nowrap'>lu_factor.hpp</span></a> contains the current source code that implements these specifications. <hr/>Input File: cppad/lu_factor.hpp </body> </html>