[CppAD] LuSolve nan

[Brad Bell]

In a previous reply I mentioned the LuVecAD alternative; see
    http://www.coin-or.org/CppAD/Doc/luvecad.xml
This reply contains an example of the retaping alternative.
The output from this alternative is
a = 2
{ 0.693147 }{ 0.5 }
a = -2
{ -Inf }{ -0.5 }
a = -2
{ 0.693147 }{ -0.5 }

Retaping Alternative Code:
==================
# include <CppAD/CppAD.h>
using namespace CppAD;
int main()
{
   // dimension of the matrix
   size_t n = 1;

   // Dummy arguments to LuSolve (do not need to be set because sequence
   // of operations does not depend on their values)
   vector<AD<double> > B(0);
   vector<AD<double> > X(0);

   // construct ADFun object F:A->log(abs(det(A))) corresponding to A = 1
   vector<AD<double> > A(n);
   vector<AD<double> > y(n);
   A[0]=1;
   Independent(A);
   AD<double> logdet;
   LuSolve(n, 0, A, B, X, logdet);
   y[0]=logdet;
   ADFun<double> F(A,y);

   // evaluate F(A) and F'(A) at A = 2
   vector<double> a(n);
   a[0]=2;
   std::cout << "a = 2\n";
   std::cout << F.Forward(0, a) ;
   std::cout << F.Jacobian(a) ;

   // evaluate F(A) and F'(A) at A = 2
   a[0]=-2;
   std::cout << "\na = -2\n";
   std::cout << F.Forward(0, a) ;
   std::cout << F.Jacobian(a) << "\n";

   // construct ADFun object G:A->log(abs(det(A))) corresponding to A = -1
   // is different from F because sequence of operations depends on A
   A[0]=-1;
   Independent(A);
   LuSolve(n, 0, A, B, X, logdet);
   y[0]=logdet;
   ADFun<double> G(A,y);

   // evaluate G(A) and G'(A) at A = -2
   a[0]=-2;
   std::cout << "\na = -2\n";
   std::cout << G.Forward(0, a) ;
   std::cout << G.Jacobian(a) << "\n";
}


Kasper Kristensen wrote:

>Hello all,
>
>I have attached a small code example which tapes the mapping
>A->log(abs(det(A))) for a 1x1 matrix A=1 using LuSolve.
>I was expecting the tape to be valid for all non-zero matrices A.
>However if I evaluate the tape with a matrix for which det(A)<0 the
>evaluation gives nan (but the Jacobian is correct!).
>
>Best regards
>Kasper
>
>PROGRAM OUTPUT:
>
>A=2
>{ 0.693147 }{ 0.5 }
>A=-2
>{ nan }{ -0.5 }
>
>CODE:
>
># include <CppAD/CppAD.h>
>using namespace CppAD;
>int main()
>{
>   // dimension of the matrix
>   size_t n = 1;
>
>   // Dummy arguments to LuSolve
>   vector<AD<double> > B(0);
>   vector<AD<double> > X(0);
>
>   // construct ADFun object F:A->log(abs(det(A)))
>   vector<AD<double> > A(n);
>   vector<AD<double> > y(n);
>   A[0]=1;
>   Independent(A);
>   AD<double> logdet;
>   LuSolve(n, 0, A, B, X, logdet);
>   y[0]=logdet;
>   ADFun<double> F(A,y);
>
>   // evaluate the tape
>   vector<double> A0(n);
>   A0[0]=2;
>   std::cout << "A=2\n";
>   std::cout << F.Forward(0,A0) ;
>   std::cout << F.Jacobian(A0) ;
>   A0[0]=-2;
>   std::cout << "\nA=-2\n";
>   std::cout << F.Forward(0,A0) ;
>   std::cout << F.Jacobian(A0) << "\n";
>
>}
>
>
>_______________________________________________
>CppAD mailing list
>CppAD at list.coin-or.org
>http://list.coin-or.org/mailman/listinfo/cppad
>
>
>  
>