Prev Next Index-> contents reference index search external Up-> CppAD utility LuDetAndSolve LuFactor lu_factor.cpp 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 LuFactor-> lu_factor.cpp lu_factor.hpp lu_factor.cpp 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}} }$
LuFactor: Example and Test
# include <cstdlib> // for rand function # include <cppad/utility/lu_factor.hpp> // for CppAD::LuFactor # include <cppad/utility/near_equal.hpp> // for CppAD::NearEqual # include <cppad/utility/vector.hpp> // for CppAD::vector bool LuFactor(void) { bool ok = true; # ifndef _MSC_VER using std::rand; using std::srand; # endif using CppAD::NearEqual; double eps99 = 99.0 * std::numeric_limits<double>::epsilon(); size_t n = 5; // number rows in A double rand_max = double(RAND_MAX); // maximum rand value double sum; // element of L * U double pij; // element of permuted A size_t i, j, k; // temporary indices // A is an n by n matrix CppAD::vector<double> A(n*n), LU(n*n), L(n*n), U(n*n); // set A equal to an n by n random matrix for(i = 0; i < n; i++) for(j = 0; j < n; j++) A[i * n + j] = rand() / rand_max; // pivot vectors CppAD::vector<size_t> ip(n); CppAD::vector<size_t> jp(n); // factor the matrix A LU = A; CppAD::LuFactor(ip, jp, LU); // check that ip and jp are permutations of the indices 0, ... , n-1 for(i = 0; i < n; i++) { ok &= (ip[i] < n); ok &= (jp[i] < n); for(j = 0; j < n; j++) { if( i != j ) { ok &= (ip[i] != ip[j]); ok &= (jp[i] != jp[j]); } } } // Extract L from LU for(i = 0; i < n; i++) { // elements along and below the diagonal for(j = 0; j <= i; j++) L[i * n + j] = LU[ ip[i] * n + jp[j] ]; // elements above the diagonal for(j = i+1; j < n; j++) L[i * n + j] = 0.; } // Extract U from LU for(i = 0; i < n; i++) { // elements below the diagonal for(j = 0; j < i; j++) U[i * n + j] = 0.; // elements along the diagonal U[i * n + i] = 1.; // elements above the diagonal for(j = i+1; j < n; j++) U[i * n + j] = LU[ ip[i] * n + jp[j] ]; } // Compute L * U for(i = 0; i < n; i++) { for(j = 0; j < n; j++) { // compute element (i,j) entry in L * U sum = 0.; for(k = 0; k < n; k++) sum += L[i * n + k] * U[k * n + j]; // element (i,j) in permuted version of A pij = A[ ip[i] * n + jp[j] ]; // compare ok &= NearEqual(pij, sum, eps99, eps99); } } return ok; } 
Input File: example/utility/lu_factor.cpp