[ADOL-C] advector example

[Antoine De Blois]

Hi Kshitij,

I started modifying my code based on your advice. I have a question before I start digging more into intense modifications!

In the stripped code posted below, the method GenerateResidualTape has the index i, j and k as arguments. Internally, there are many advector that use the index as adouble.  As you can see, the array W (the independent variables) changes as a function of i, j and k (GetPrimitives).  So far, so good. But how about all the other internal advectors that depend on the indices? How can the call to sparse_jac know about those changed indices (apart from W obviously)? For example, I have an advector of face areas. In GenerateResidualTape, I retrieve the areas of a given I,j,k. On the next iteration, I want the jacobian with respect to the new areas from the advector.

Thank you for your time,
Antoine


        int cnt=0;
        for (adouble k=2;k<mNmax[2]-2;k++) // i, j and k used to be int
        {
               for (adouble j=2;j<mNmax[1]-2;j++)
               {
                       for (adouble i=2;i<mNmax[0]-2;i++)
                       {
                               int ret;
                               if (cnt!=0)
                               {
                                      mpStencil->GetPrimitives(i,j,k,W);
                                      ret = sparse_jac(mTag, m, n, 0, W, &nnz, &rind, &cind, &values, options);
                                      free(rind); rind=NULL;
                                      free(cind); cind=NULL;
                                      free(values); values=NULL;
                               }

                               if ( (cnt==0) || (ret!=3))
                               {
                                      GenerateResidualTape(i,j,k);
                                      mpStencil->GetPrimitives(i,j,k,W);
                                      ret = sparse_jac(mTag, m, n, 0, W, &nnz, &rind, &cind, &values, options);
                                      free(rind); rind=NULL;
                                      free(cind); cind=NULL;
                                      free(values); values=NULL;
                               }
                               cnt++;
                       }
               }
        }



----------------------------------

Hello,

advector can indeed deal with changing indices in an array at a
different evaluation point.

The constructor advector(const std::vector<adouble>& v) will indeed copy
the elements. This is standard practice when dealing with STL classes.
If it did not copy then there would be the chance of some other part of
the code changing the memory contents from under the STL-container, and
there would be no warning for this.

The easiest way to do this efficiently is to construct the advector
using advector(size_t n) and then assigning the elements using the
operator [](size_t i) or using a dynamic_cast<vector<adouble>& >() to
access the std::vector inside the advector. This would avoid copying.
You anyway need to assign the elements at some point.

The index for accessing the elements, which may change at a different
evaluation point must however be an adouble or an expression with
adoubles and the indexing is done with operator [](const badouble& i)

Hope this helps
Kshitij

As on 2013-03-22 22:24, Antoine De Blois did write:
> Hi everyone,
>
> I was trying to find a way to avoid retaping when the index of an array changes  for a different taping point. I found today the advector, which sounds very promising for my application.
>
> I want to have the advector 's as members of a class and I want to efficiently construct the advector. I don't think that the
> advector(const std::vector<adouble>& v) will be efficient due to the involved copies. Is it done this way to avoid the lvalue issue?
>
> Is there any C++ code example someone is willing to share? I tried to search quickly on the web but I could not find any. This article confirms that ADOL-C can do what I want:
> "Computing derivatives in a meshless simulation using permuation in ADOL-C." Kulshreshtha
>
> Thank you for your time,
> A

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://list.coin-or.org/pipermail/adol-c/attachments/20130327/dbc26c8f/attachment.html>