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}} }@)@
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