Hi Antoine,<div><br></div><div>Do you still have the same problem? Our group is maintaining ColPack and if there is any bug within ColPack or related with ADOL-C+ColPack, we are glad to help. So if it's fine, you can send the code and we will check it. <br>
<br>Sincerely,</div><div>Mu</div><div><br><div class="gmail_quote">On Wed, Mar 27, 2013 at 9:43 PM, Antoine De Blois <span dir="ltr"><<a href="mailto:antoine.deblois@aero.bombardier.com" target="_blank">antoine.deblois@aero.bombardier.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div lang="EN-US" link="blue" vlink="purple">
<div>
<p class="MsoNormal">Hi Kshitij,<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">I started modifying my code based on your advice. I have a question before I start digging more into intense modifications!
<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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.<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">Thank you for your time,<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">Antoine<u></u><u></u></span></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> int cnt=0;<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> for (adouble k=2;k<mNmax[2]-2;k++)
<b><span style="color:#4f6228">// i, j and k used to be int</span></b><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> {<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> for (adouble j=2;j<mNmax[1]-2;j++)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> {<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> for (adouble i=2;i<mNmax[0]-2;i++)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> {<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> int ret;<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> if (cnt!=0)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> {
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> mpStencil->GetPrimitives(i,j,k,W);<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> ret = sparse_jac(mTag, m, n, 0, W, &nnz, &rind, &cind, &values, options);<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> free(rind); rind=NULL;<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> free(cind); cind=NULL;<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> free(values); values=NULL;
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> }<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> if ( (cnt==0) || (ret!=3))<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> {
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> GenerateResidualTape(i,j,k);<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> mpStencil->GetPrimitives(i,j,k,W);
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> ret = sparse_jac(mTag, m, n, 0, W, &nnz, &rind, &cind, &values, options);
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> free(rind); rind=NULL;<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> free(cind); cind=NULL;<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> free(values); values=NULL;<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> }<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> cnt++;<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> }<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> }<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""> }<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">----------------------------------<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">Hello,<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">advector can indeed deal with changing indices in an array at a<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">different evaluation point.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">The constructor advector(const std::vector<adouble>& v) will indeed copy<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">the elements. This is standard practice when dealing with STL classes.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">If it did not copy then there would be the chance of some other part of<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">the code changing the memory contents from under the STL-container, and<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">there would be no warning for this.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">The easiest way to do this efficiently is to construct the advector<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">using advector(size_t n) and then assigning the elements using the<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">operator [](size_t i) or using a dynamic_cast<vector<adouble>& >() to<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">access the std::vector inside the advector. This would avoid copying.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">You anyway need to assign the elements at some point.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">The index for accessing the elements, which may change at a different<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">evaluation point must however be an adouble or an expression with<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">adoubles and the indexing is done with operator [](const badouble& i)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">Hope this helps<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">Kshitij<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New""><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">As on 2013-03-22 22:24, Antoine De Blois did write:<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i> Hi everyone,<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i>
<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i> 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.<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i>
<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i> 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<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i> advector(const std::vector<adouble>& v) will be efficient due to the involved copies. Is it done this way to avoid the lvalue issue?<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i>
<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i> 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:<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i> "Computing derivatives in a meshless simulation using permuation in ADOL-C." Kulshreshtha<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i>
<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i> Thank you for your time,<u></u><u></u></i></span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;font-family:"Courier New"">><i> A</i><u></u><u></u></span></p>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
</div>
<br>_______________________________________________<br>
ADOL-C mailing list<br>
<a href="mailto:ADOL-C@list.coin-or.org">ADOL-C@list.coin-or.org</a><br>
<a href="http://list.coin-or.org/mailman/listinfo/adol-c" target="_blank">http://list.coin-or.org/mailman/listinfo/adol-c</a><br></blockquote></div><br></div>