ADMB Documentation  11.1.2403
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2prd.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2prd.cpp 2350 2014-09-15 23:19:51Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 #ifndef OPT_LIB
00013   #include <cassert>
00014   #include <climits>
00015 #endif
00016 
00021   df1b2variable operator * (const df1b2variable& _x,const df1b2variable& _y)
00022   {
00023     ADUNCONST(df1b2variable,x)
00024     ADUNCONST(df1b2variable,y)
00025     df1b2variable z;
00026     double * xd=x.get_u_dot();
00027     double * yd=y.get_u_dot();
00028     double * zd=z.get_u_dot();
00029     double xu=*x.get_u();
00030     double yu=*y.get_u();
00031 
00032     *z.get_u()=xu*yu;
00033 
00034     for (int i=0;i<df1b2variable::nvar;i++)
00035     {
00036       *zd++ = yu * *xd++ + xu * *yd++;
00037     }
00038 
00039     // WRITE WHATEVER ON TAPE
00040     if (!df1b2_gradlist::no_derivatives)
00041       f1b2gradlist->write_pass1_prod(&x,&y,&z);
00042     return z;
00043   }
00044 
00045 void ad_read_pass2_prod(void);
00046 
00051  int df1b2_gradlist::write_pass1_prod(const df1b2variable * _px,
00052    const df1b2variable * _py,df1b2variable * pz)
00053  {
00054    ADUNCONST(df1b2variable*,px)
00055    ADUNCONST(df1b2variable*,py)
00056    ncount++;
00057 #if defined(CHECK_COUNT)
00058   if (ncount >= ncount_check)
00059     cout << ncount << endl;
00060 #endif
00061    int nvar=df1b2variable::nvar;
00062 
00063    //int total_bytes=3*sizeof(df1b2_header)+sizeof(char*)
00064    //  +2*(nvar+1)*sizeof(double);
00065    size_t total_bytes=3*sizeof(df1b2_header)
00066      +2*(nvar+1)*sizeof(double);
00067 // string identifier debug stuff
00068 #if defined(SAFE_ALL)
00069   char ids[]="DL";
00070   size_t slen=strlen(ids);
00071   total_bytes+=slen;
00072 #endif
00073 
00074 #ifndef OPT_LIB
00075   assert(total_bytes <= INT_MAX);
00076 #endif
00077   list.check_buffer_size((int)total_bytes);
00078   void * tmpptr=list.bptr;
00079 #if defined(SAFE_ALL)
00080   memcpy(list,ids,slen);
00081 #endif
00082 // end of string identifier debug stuff
00083 
00084    memcpy(list,(df1b2_header*)(px),sizeof(df1b2_header));
00085    memcpy(list,(df1b2_header*)(py),sizeof(df1b2_header));
00086    memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00087    //memcpy(list,&pf,sizeof(char *));
00088    //*(char**)(list.bptr)=(char*)pf;
00089   const int sizeofdouble = sizeof(double);
00090    memcpy(list,px->get_u(),sizeofdouble);
00091    memcpy(list,py->get_u(),sizeofdouble);
00092    memcpy(list,px->get_u_dot(),nvar*sizeofdouble);
00093    memcpy(list,py->get_u_dot(),nvar*sizeofdouble);
00094    // ***** write  record size
00095    nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00096    nlist.bptr->pf=(ADrfptr)(&ad_read_pass2_prod);
00097       ++nlist;
00098    return 0;
00099  }
00100 
00101 
00102 void read_pass2_1_prod(void);
00103 void read_pass2_2_prod(void);
00104 void read_pass2_3_prod(void);
00105 
00110 void ad_read_pass2_prod(void)
00111 {
00112   switch(df1b2variable::passnumber)
00113   {
00114   case 1:
00115     read_pass2_1_prod();
00116     break;
00117   case 2:
00118     read_pass2_2_prod();
00119     break;
00120   case 3:
00121     read_pass2_3_prod();
00122     break;
00123   default:
00124     cerr << "illegal value for df1b2variable::pass = "
00125          << df1b2variable::passnumber << endl;
00126     exit(1);
00127   }
00128 }
00129 
00134 void read_pass2_1_prod(void)
00135 {
00136   // We are going backword for bptr and nbptr
00137   // and  forward for bptr2 and nbptr2
00138   // the current entry+2 in bptr is the size of the record i.e
00139   // points to the next record
00140   //char * bptr=f1b2gradlist->bptr;
00141   //char * bptr2=f1b2gradlist2->bptr;
00142   int nvar=df1b2variable::nvar;
00143   test_smartlist& list=f1b2gradlist->list;
00144   //f1b2gradlist->nlist-=sizeof(int);
00145   int num_bytes=f1b2gradlist->nlist.bptr->numbytes;
00146   list-=num_bytes;
00147   list.saveposition(); // save pointer to beginning of record;
00148   double xu,yu;
00149   //ad_dstar xdot,ydot;
00150   //df1b2function2 * pf;
00151 
00152   // get info from tape1
00153 #if defined(SAFE_ALL)
00154   checkidentiferstring("DL",f1b2gradlist->list);
00155 #endif
00156   char * bptr=f1b2gradlist->list.bptr;
00157   df1b2_header * px=(df1b2_header *) bptr;
00158   bptr+=sizeof(df1b2_header);
00159   df1b2_header * py=(df1b2_header *) bptr;
00160   bptr+=sizeof(df1b2_header);
00161   df1b2_header * pz=(df1b2_header *) bptr;
00162   bptr+=sizeof(df1b2_header);
00163   //pf=*(df1b2function2 **) bptr;
00164   //bptr+=sizeof(char*);
00165   memcpy(&xu,bptr,sizeof(double));
00166   bptr+=sizeof(double);
00167   memcpy(&yu,bptr,sizeof(double));
00168   bptr+=sizeof(double);
00169   double * xdot=(double*)bptr;
00170   bptr+=nvar*sizeof(double);
00171   double * ydot=(double*)bptr;
00172 
00173   list.restoreposition(); // save pointer to beginning of record;
00174 
00175   // ****************************************************************
00176   // turn this off if no third derivatives are calculated
00177   // if (!no_third_derivatives)
00178   // {
00179   // save for second reverse pass
00180   // save identifier 1
00181      test_smartlist & list2 = f1b2gradlist->list2;
00182 
00183   size_t total_bytes=2*nvar*sizeof(double);
00184 // string identifier debug stuff
00185 #if defined(SAFE_ALL)
00186   char ids[]="QK";
00187   size_t slen=strlen(ids);
00188   total_bytes+=slen;
00189 #endif
00190 
00191 #ifndef OPT_LIB
00192   assert(total_bytes <= INT_MAX);
00193 #endif
00194 
00195   list2.check_buffer_size((int)total_bytes);
00196 
00197   void * tmpptr=list2.bptr;
00198 #if defined(SAFE_ALL)
00199   memcpy(list2,ids,slen);
00200 #endif
00201 
00202   fixed_smartlist2 & nlist2 = f1b2gradlist->nlist2;
00203   const int sizeofdouble = sizeof(double);
00204   memcpy(list2,pz->get_u_bar(),nvar*sizeofdouble);
00205   memcpy(list2,pz->get_u_dot_bar(),nvar*sizeofdouble);
00206   *nlist2.bptr=adptr_diff(list2.bptr,tmpptr);
00207   ++nlist2;
00208 
00209   // Do first reverse pass calculations
00210   for (int i=0;i<nvar;i++)
00211   {
00212    //px->u_bar[i]+=(pf->df1)(xu,yu)*pz->u_bar[i];
00213     px->u_bar[i]+=yu*pz->u_bar[i];
00214   }
00215   for (int i=0;i<nvar;i++)
00216   {
00217     //py->u_bar[i]+=(pf->df2)(xu,yu)*pz->u_bar[i];
00218     py->u_bar[i]+=xu*pz->u_bar[i];
00219   }
00220 
00221   for (int i=0;i<nvar;i++)
00222   {
00223     //px->u_bar[i]+=(pf->d2f11)(xu,yu)*xdot[i]*pz->u_dot_bar[i];
00224     //px->u_bar[i]+=(pf->d2f12)(xu,yu)*ydot[i]*pz->u_dot_bar[i];
00225     px->u_bar[i]+=ydot[i]*pz->u_dot_bar[i];
00226   }
00227 
00228   for (int i=0;i<nvar;i++)
00229   {
00230     //py->u_bar[i]+=(pf->d2f22)(xu,yu)*ydot[i]*pz->u_dot_bar[i];
00231     //py->u_bar[i]+=(pf->d2f12)(xu,yu)*xdot[i]*pz->u_dot_bar[i];
00232     py->u_bar[i]+=xdot[i]*pz->u_dot_bar[i];
00233   }
00234   for (int i=0;i<nvar;i++)
00235   {
00236     //px->u_dot_bar[i]+=(pf->df1)(xu,yu)*pz->u_dot_bar[i];
00237     px->u_dot_bar[i]+=yu*pz->u_dot_bar[i];
00238   }
00239   for (int i=0;i<nvar;i++)
00240   {
00241     //py->u_dot_bar[i]+=(pf->df2)(xu,yu)*pz->u_dot_bar[i];
00242     py->u_dot_bar[i]+=xu*pz->u_dot_bar[i];
00243   }
00244 
00245   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00246   for (int i=0;i<nvar;i++)
00247   {
00248     pz->u_bar[i]=0;
00249   }
00250   for (int i=0;i<nvar;i++)
00251   {
00252     pz->u_dot_bar[i]=0;
00253   }
00254 }
00255 
00260 void read_pass2_2_prod(void)
00261 {
00262   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00263   // We are going forward for bptr and backword for bptr2
00264   //
00265   // list 1
00266   //
00267   int nvar=df1b2variable::nvar;
00268   test_smartlist & list=f1b2gradlist->list;
00269 
00270   //int total_bytes=3*sizeof(df1b2_header)+sizeof(char*)
00271   //  +2*(nvar+1)*sizeof(double);
00272   size_t total_bytes=3*sizeof(df1b2_header)
00273     +2*(nvar+1)*sizeof(double);
00274 // string identifier debug stuff
00275 #if defined(SAFE_ALL)
00276   char ids[]="DL";
00277   size_t slen=strlen(ids);
00278   total_bytes+=slen;
00279 #endif
00280 
00281 #ifndef OPT_LIB
00282   assert(total_bytes <= INT_MAX);
00283 #endif
00284 
00285   list.check_buffer_size((int)total_bytes);
00286 // end of string identifier debug stuff
00287 
00288   list.saveposition(); // save pointer to beginning of record;
00289   fixed_smartlist & nlist=f1b2gradlist->nlist;
00290    // nlist-=sizeof(int);
00291   // get record size
00292   int num_bytes=nlist.bptr->numbytes;
00293   //
00294   // list 2
00295   //
00296   test_smartlist & list2=f1b2gradlist->list2;
00297   fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00298   // get record size
00299   int num_bytes2=*nlist2.bptr;
00300   --nlist2;
00301   // backup the size of the record
00302   list2-=num_bytes2;
00303   list2.saveposition(); // save pointer to beginning of record;
00304   // save the pointer to the beginning of the record
00305   // bptr and bptr2 now both point to the beginning of their records
00306 
00307   double xu,yu;
00308   double * xdot;
00309   double * ydot;
00310   //df1b2_header x,z;
00311   //df1b2function2 * pf;
00312 
00313   // get info from tape1
00314   // get info from tape1
00315 #if defined(SAFE_ALL)
00316   checkidentiferstring("DL",list);
00317   checkidentiferstring("QK",list2);
00318 #endif
00319   df1b2_header * px=(df1b2_header *) list.bptr;
00320   list.bptr+=sizeof(df1b2_header);
00321   df1b2_header * py=(df1b2_header *) list.bptr;
00322   list.bptr+=sizeof(df1b2_header);
00323   df1b2_header * pz=(df1b2_header *) list.bptr;
00324   list.bptr+=sizeof(df1b2_header);
00325   //pf=*(df1b2function2 **) list.bptr;
00326   //list.bptr+=sizeof(char*);
00327   memcpy(&xu,list.bptr,sizeof(double));
00328   list.bptr+=sizeof(double);
00329   memcpy(&yu,list.bptr,sizeof(double));
00330   list.bptr+=sizeof(double);
00331   xdot=(double*)list.bptr;
00332   list.bptr+=nvar*sizeof(double);
00333   ydot=(double*)list.bptr;
00334   list.restoreposition(num_bytes); // save pointer to beginning of record;
00335 
00336   double * zbar;
00337   double * zdotbar;
00338 
00339   zbar=(double*)list2.bptr;
00340   zdotbar=(double*)(list2.bptr+nvar*sizeof(double));
00341   list2.restoreposition(); // save pointer to beginning of record;
00342 
00343   double * x_tilde=px->get_u_tilde();
00344   double * x_dot_tilde=px->get_u_dot_tilde();
00345   double * x_bar_tilde=px->get_u_bar_tilde();
00346   double * x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00347   double * y_tilde=py->get_u_tilde();
00348   double * y_dot_tilde=py->get_u_dot_tilde();
00349   double * y_bar_tilde=py->get_u_bar_tilde();
00350   double * y_dot_bar_tilde=py->get_u_dot_bar_tilde();
00351   double * z_bar_tilde=pz->get_u_bar_tilde();
00352   double * z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00353   // Do second "reverse-reverse" pass calculations
00354 
00355   for (int i=0;i<nvar;i++)
00356   {
00357     z_bar_tilde[i]=0;
00358     z_dot_bar_tilde[i]=0;
00359   }
00360 
00361   // start with x and add y
00362   for (int i=0;i<nvar;i++)
00363   {
00364     //*x_tilde+=(pf->d2f11)(xu,yu)*zbar[i]*x_bar_tilde[i];
00365     //z_bar_tilde[i]+=(pf->df1)(xu,yu)*x_bar_tilde[i];
00366     //*y_tilde+=(pf->d2f12)(xu,yu)*zbar[i]*x_bar_tilde[i];
00367     z_bar_tilde[i]+=yu*x_bar_tilde[i];
00368     *y_tilde+=zbar[i]*x_bar_tilde[i];
00369   }
00370 
00371   for (int i=0;i<nvar;i++)
00372   {
00373     //*x_tilde+=(pf->d2f11)(xu,yu)*zdotbar[i]*x_dot_bar_tilde[i];
00374     //*y_tilde+=(pf->d2f12)(xu,yu)*zdotbar[i]*x_dot_bar_tilde[i];
00375     //z_dot_bar_tilde[i]+=(pf->df1)(xu,yu)*x_dot_bar_tilde[i];
00376     *y_tilde+=zdotbar[i]*x_dot_bar_tilde[i];
00377     z_dot_bar_tilde[i]+=yu*x_dot_bar_tilde[i];
00378   }
00379 
00380   /*
00381   for (i=0;i<nvar;i++)
00382   {
00383     x_dot_tilde[i]+=(pf->d2f11)(xu,yu)*zdotbar[i]*x_bar_tilde[i];
00384     z_dot_bar_tilde[i]+=(pf->d2f11)(xu,yu)*xdot[i]*x_bar_tilde[i];
00385     *x_tilde+=(pf->d3f111)(xu,yu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00386     *y_tilde+=(pf->d3f112)(xu,yu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00387   }
00388   */
00389   // start with y and add x
00390   for (int i=0;i<nvar;i++)
00391   {
00392     //*y_tilde+=(pf->d2f22)(xu,yu)*zbar[i]*y_bar_tilde[i];
00393     //*x_tilde+=(pf->d2f12)(xu,yu)*zbar[i]*y_bar_tilde[i];
00394     //z_bar_tilde[i]+=(pf->df2)(xu,yu)*y_bar_tilde[i];
00395     *x_tilde+=zbar[i]*y_bar_tilde[i];
00396     z_bar_tilde[i]+=xu*y_bar_tilde[i];
00397   }
00398 
00399   for (int i=0;i<nvar;i++)
00400   {
00401     //*y_tilde+=(pf->d2f22)(xu,yu)*zdotbar[i]*y_dot_bar_tilde[i];
00402     //*x_tilde+=(pf->d2f12)(xu,yu)*zdotbar[i]*y_dot_bar_tilde[i];
00403     //z_dot_bar_tilde[i]+=(pf->df2)(xu,yu)*y_dot_bar_tilde[i];
00404     *x_tilde+=zdotbar[i]*y_dot_bar_tilde[i];
00405     z_dot_bar_tilde[i]+=xu*y_dot_bar_tilde[i];
00406   }
00407 
00408   /*
00409   for (i=0;i<nvar;i++)
00410   {
00411     y_dot_tilde[i]+=(pf->d2f22)(xu,yu)*zdotbar[i]*y_bar_tilde[i];
00412     z_dot_bar_tilde[i]+=(pf->d2f22)(xu,yu)*ydot[i]*y_bar_tilde[i];
00413     *y_tilde+=(pf->d3f222)(xu,yu)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00414     *x_tilde+=(pf->d3f122)(xu,yu)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00415   }
00416   */
00417   for (int i=0;i<nvar;i++)
00418   {
00419     //*x_tilde+=(pf->d3f112)(xu,yu)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00420     //*y_tilde+=(pf->d3f122)(xu,yu)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00421     //y_dot_tilde[i]+=(pf->d2f12)(xu,yu)*zdotbar[i]*x_bar_tilde[i];
00422     //z_dot_bar_tilde[i]+=(pf->d2f12)(xu,yu)*ydot[i]*x_bar_tilde[i];
00423     y_dot_tilde[i]+=zdotbar[i]*x_bar_tilde[i];
00424     z_dot_bar_tilde[i]+=ydot[i]*x_bar_tilde[i];
00425   }
00426   for (int i=0;i<nvar;i++)
00427   {
00428     //*x_tilde+=(pf->d3f112)(xu,yu)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00429     //*y_tilde+=(pf->d3f122)(xu,yu)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00430     //x_dot_tilde[i]+=(pf->d2f12)(xu,yu)*zdotbar[i]*y_bar_tilde[i];
00431     //z_dot_bar_tilde[i]+=(pf->d2f12)(xu,yu)*xdot[i]*y_bar_tilde[i];
00432     x_dot_tilde[i]+=zdotbar[i]*y_bar_tilde[i];
00433     z_dot_bar_tilde[i]+=xdot[i]*y_bar_tilde[i];
00434   }
00435 }
00436 
00441 void read_pass2_3_prod(void)
00442 {
00443   // We are going backword for bptr and forward for bptr2
00444   // the current entry+2 in bptr is the size of the record i.e
00445   // points to the next record
00446   int nvar=df1b2variable::nvar;
00447   fixed_smartlist & nlist=f1b2gradlist->nlist;
00448   test_smartlist& list=f1b2gradlist->list;
00449    // nlist-=sizeof(int);
00450   // get record size
00451   int num_bytes=nlist.bptr->numbytes;
00452   // backup the size of the record
00453   list-=num_bytes;
00454   list.saveposition(); // save pointer to beginning of record;
00455   // save the pointer to the beginning of the record
00456   double xu;
00457   double yu;
00458   double * xdot;
00459   double * ydot;
00460   //df1b2_header x,z;
00461   //df1b2function2 * pf;
00462 
00463   // get info from tape1
00464   // get info from tape1
00465 #if defined(SAFE_ALL)
00466   checkidentiferstring("DL",list);
00467 #endif
00468   df1b2_header * px=(df1b2_header *) list.bptr;
00469   list.bptr+=sizeof(df1b2_header);
00470   df1b2_header * py=(df1b2_header *) list.bptr;
00471   list.bptr+=sizeof(df1b2_header);
00472   df1b2_header * pz=(df1b2_header *) list.bptr;
00473   list.bptr+=sizeof(df1b2_header);
00474   //pf=*(df1b2function2 **) list.bptr;
00475   //list.bptr+=sizeof(char*);
00476   memcpy(&xu,list.bptr,sizeof(double));
00477   list.bptr+=sizeof(double);
00478   memcpy(&yu,list.bptr,sizeof(double));
00479   list.bptr+=sizeof(double);
00480   xdot=(double*)list.bptr;
00481   list.bptr+=nvar*sizeof(double);
00482   ydot=(double*)list.bptr;
00483   list.restoreposition(); // save pointer to beginning of record;
00484 
00485   //*(px->u_tilde)+=(pf->df1)(xu,yu)* *(pz->u_tilde);
00486   //*(py->u_tilde)+=(pf->df2)(xu,yu)* *(pz->u_tilde);
00487   *(px->u_tilde)+=yu* *(pz->u_tilde);
00488   *(py->u_tilde)+=xu* *(pz->u_tilde);
00489   for (int i=0;i<nvar;i++)
00490   {
00491     //*(px->u_tilde)+=(pf->d2f11)(xu,yu)*xdot[i]*pz->u_dot_tilde[i];
00492     //*(py->u_tilde)+=(pf->d2f12)(xu,yu)*xdot[i]*pz->u_dot_tilde[i];
00493     *(py->u_tilde)+=xdot[i]*pz->u_dot_tilde[i];
00494     //*(py->u_tilde)+=(pf->d2f22)(xu,yu)*ydot[i]*pz->u_dot_tilde[i];
00495     //*(px->u_tilde)+=(pf->d2f12)(xu,yu)*ydot[i]*pz->u_dot_tilde[i];
00496     *(px->u_tilde)+=ydot[i]*pz->u_dot_tilde[i];
00497   }
00498   for (int i=0;i<nvar;i++)
00499   {
00500     //px->u_dot_tilde[i]+=(pf->df1)(xu,yu)*pz->u_dot_tilde[i];
00501     //py->u_dot_tilde[i]+=(pf->df2)(xu,yu)*pz->u_dot_tilde[i];
00502     px->u_dot_tilde[i]+=yu*pz->u_dot_tilde[i];
00503     py->u_dot_tilde[i]+=xu*pz->u_dot_tilde[i];
00504   }
00505   *(pz->u_tilde)=0;
00506   for (int i=0;i<nvar;i++)
00507   {
00508     pz->u_dot_tilde[i]=0;
00509   }
00510 }