Prev Next 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}} }@)@ This is cppad-20221105 documentation. Here is a link to its current documentation .
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