ADMB Documentation  11.1x.2717
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2prd.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2prd.cpp 2601 2014-11-08 22:34:58Z 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 (unsigned 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    unsigned 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   list.check_buffer_size(total_bytes);
00075   void * tmpptr=list.bptr;
00076 #if defined(SAFE_ALL)
00077   memcpy(list,ids,slen);
00078 #endif
00079 // end of string identifier debug stuff
00080 
00081    memcpy(list,(df1b2_header*)(px),sizeof(df1b2_header));
00082    memcpy(list,(df1b2_header*)(py),sizeof(df1b2_header));
00083    memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00084    //memcpy(list,&pf,sizeof(char *));
00085    //*(char**)(list.bptr)=(char*)pf;
00086   const int sizeofdouble = sizeof(double);
00087    memcpy(list,px->get_u(),sizeofdouble);
00088    memcpy(list,py->get_u(),sizeofdouble);
00089    memcpy(list,px->get_u_dot(),nvar*sizeofdouble);
00090    memcpy(list,py->get_u_dot(),nvar*sizeofdouble);
00091    // ***** write  record size
00092    nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00093    nlist.bptr->pf=(ADrfptr)(&ad_read_pass2_prod);
00094       ++nlist;
00095    return 0;
00096  }
00097 
00098 
00099 void read_pass2_1_prod(void);
00100 void read_pass2_2_prod(void);
00101 void read_pass2_3_prod(void);
00102 
00107 void ad_read_pass2_prod(void)
00108 {
00109   switch(df1b2variable::passnumber)
00110   {
00111   case 1:
00112     read_pass2_1_prod();
00113     break;
00114   case 2:
00115     read_pass2_2_prod();
00116     break;
00117   case 3:
00118     read_pass2_3_prod();
00119     break;
00120   default:
00121     cerr << "illegal value for df1b2variable::pass = "
00122          << df1b2variable::passnumber << endl;
00123     exit(1);
00124   }
00125 }
00126 
00131 void read_pass2_1_prod(void)
00132 {
00133   // We are going backword for bptr and nbptr
00134   // and  forward for bptr2 and nbptr2
00135   // the current entry+2 in bptr is the size of the record i.e
00136   // points to the next record
00137   //char * bptr=f1b2gradlist->bptr;
00138   //char * bptr2=f1b2gradlist2->bptr;
00139   unsigned int nvar=df1b2variable::nvar;
00140   test_smartlist& list=f1b2gradlist->list;
00141   //f1b2gradlist->nlist-=sizeof(int);
00142   int num_bytes=f1b2gradlist->nlist.bptr->numbytes;
00143   list-=num_bytes;
00144   list.saveposition(); // save pointer to beginning of record;
00145   double xu,yu;
00146   //ad_dstar xdot,ydot;
00147   //df1b2function2 * pf;
00148 
00149   // get info from tape1
00150 #if defined(SAFE_ALL)
00151   checkidentiferstring("DL",f1b2gradlist->list);
00152 #endif
00153   char * bptr=f1b2gradlist->list.bptr;
00154   df1b2_header * px=(df1b2_header *) bptr;
00155   bptr+=sizeof(df1b2_header);
00156   df1b2_header * py=(df1b2_header *) bptr;
00157   bptr+=sizeof(df1b2_header);
00158   df1b2_header * pz=(df1b2_header *) bptr;
00159   bptr+=sizeof(df1b2_header);
00160   //pf=*(df1b2function2 **) bptr;
00161   //bptr+=sizeof(char*);
00162   memcpy(&xu,bptr,sizeof(double));
00163   bptr+=sizeof(double);
00164   memcpy(&yu,bptr,sizeof(double));
00165   bptr+=sizeof(double);
00166   double * xdot=(double*)bptr;
00167   bptr+=nvar*sizeof(double);
00168   double * ydot=(double*)bptr;
00169 
00170   list.restoreposition(); // save pointer to beginning of record;
00171 
00172   // ****************************************************************
00173   // turn this off if no third derivatives are calculated
00174   // if (!no_third_derivatives)
00175   // {
00176   // save for second reverse pass
00177   // save identifier 1
00178      test_smartlist & list2 = f1b2gradlist->list2;
00179 
00180   size_t total_bytes=2*nvar*sizeof(double);
00181 // string identifier debug stuff
00182 #if defined(SAFE_ALL)
00183   char ids[]="QK";
00184   size_t slen=strlen(ids);
00185   total_bytes+=slen;
00186 #endif
00187 
00188   list2.check_buffer_size(total_bytes);
00189 
00190   void * tmpptr=list2.bptr;
00191 #if defined(SAFE_ALL)
00192   memcpy(list2,ids,slen);
00193 #endif
00194 
00195   fixed_smartlist2 & nlist2 = f1b2gradlist->nlist2;
00196   const int sizeofdouble = sizeof(double);
00197   memcpy(list2,pz->get_u_bar(),nvar*sizeofdouble);
00198   memcpy(list2,pz->get_u_dot_bar(),nvar*sizeofdouble);
00199   *nlist2.bptr=adptr_diff(list2.bptr,tmpptr);
00200   ++nlist2;
00201 
00202   // Do first reverse pass calculations
00203   for (size_t i=0;i<nvar;i++)
00204   {
00205    //px->u_bar[i]+=(pf->df1)(xu,yu)*pz->u_bar[i];
00206     px->u_bar[i]+=yu*pz->u_bar[i];
00207   }
00208   for (size_t i=0;i<nvar;i++)
00209   {
00210     //py->u_bar[i]+=(pf->df2)(xu,yu)*pz->u_bar[i];
00211     py->u_bar[i]+=xu*pz->u_bar[i];
00212   }
00213 
00214   for (size_t i=0;i<nvar;i++)
00215   {
00216     //px->u_bar[i]+=(pf->d2f11)(xu,yu)*xdot[i]*pz->u_dot_bar[i];
00217     //px->u_bar[i]+=(pf->d2f12)(xu,yu)*ydot[i]*pz->u_dot_bar[i];
00218     px->u_bar[i]+=ydot[i]*pz->u_dot_bar[i];
00219   }
00220 
00221   for (size_t i=0;i<nvar;i++)
00222   {
00223     //py->u_bar[i]+=(pf->d2f22)(xu,yu)*ydot[i]*pz->u_dot_bar[i];
00224     //py->u_bar[i]+=(pf->d2f12)(xu,yu)*xdot[i]*pz->u_dot_bar[i];
00225     py->u_bar[i]+=xdot[i]*pz->u_dot_bar[i];
00226   }
00227   for (size_t i=0;i<nvar;i++)
00228   {
00229     //px->u_dot_bar[i]+=(pf->df1)(xu,yu)*pz->u_dot_bar[i];
00230     px->u_dot_bar[i]+=yu*pz->u_dot_bar[i];
00231   }
00232   for (size_t i=0;i<nvar;i++)
00233   {
00234     //py->u_dot_bar[i]+=(pf->df2)(xu,yu)*pz->u_dot_bar[i];
00235     py->u_dot_bar[i]+=xu*pz->u_dot_bar[i];
00236   }
00237 
00238   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00239   for (size_t i=0;i<nvar;i++)
00240   {
00241     pz->u_bar[i]=0;
00242   }
00243   for (size_t i=0;i<nvar;i++)
00244   {
00245     pz->u_dot_bar[i]=0;
00246   }
00247 }
00248 
00253 void read_pass2_2_prod(void)
00254 {
00255   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00256   // We are going forward for bptr and backword for bptr2
00257   //
00258   // list 1
00259   //
00260   unsigned int nvar=df1b2variable::nvar;
00261   test_smartlist & list=f1b2gradlist->list;
00262 
00263   //int total_bytes=3*sizeof(df1b2_header)+sizeof(char*)
00264   //  +2*(nvar+1)*sizeof(double);
00265   size_t total_bytes=3*sizeof(df1b2_header)
00266     +2*(nvar+1)*sizeof(double);
00267 // string identifier debug stuff
00268 #if defined(SAFE_ALL)
00269   char ids[]="DL";
00270   size_t slen=strlen(ids);
00271   total_bytes+=slen;
00272 #endif
00273 
00274   list.check_buffer_size(total_bytes);
00275 // end of string identifier debug stuff
00276 
00277   list.saveposition(); // save pointer to beginning of record;
00278   fixed_smartlist & nlist=f1b2gradlist->nlist;
00279    // nlist-=sizeof(int);
00280   // get record size
00281   int num_bytes=nlist.bptr->numbytes;
00282   //
00283   // list 2
00284   //
00285   test_smartlist & list2=f1b2gradlist->list2;
00286   fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00287   // get record size
00288   int num_bytes2=*nlist2.bptr;
00289   --nlist2;
00290   // backup the size of the record
00291   list2-=num_bytes2;
00292   list2.saveposition(); // save pointer to beginning of record;
00293   // save the pointer to the beginning of the record
00294   // bptr and bptr2 now both point to the beginning of their records
00295 
00296   double xu,yu;
00297   double * xdot;
00298   double * ydot;
00299   //df1b2_header x,z;
00300   //df1b2function2 * pf;
00301 
00302   // get info from tape1
00303   // get info from tape1
00304 #if defined(SAFE_ALL)
00305   checkidentiferstring("DL",list);
00306   checkidentiferstring("QK",list2);
00307 #endif
00308   df1b2_header * px=(df1b2_header *) list.bptr;
00309   list.bptr+=sizeof(df1b2_header);
00310   df1b2_header * py=(df1b2_header *) list.bptr;
00311   list.bptr+=sizeof(df1b2_header);
00312   df1b2_header * pz=(df1b2_header *) list.bptr;
00313   list.bptr+=sizeof(df1b2_header);
00314   //pf=*(df1b2function2 **) list.bptr;
00315   //list.bptr+=sizeof(char*);
00316   memcpy(&xu,list.bptr,sizeof(double));
00317   list.bptr+=sizeof(double);
00318   memcpy(&yu,list.bptr,sizeof(double));
00319   list.bptr+=sizeof(double);
00320   xdot=(double*)list.bptr;
00321   list.bptr+=nvar*sizeof(double);
00322   ydot=(double*)list.bptr;
00323   list.restoreposition(num_bytes); // save pointer to beginning of record;
00324 
00325   double * zbar;
00326   double * zdotbar;
00327 
00328   zbar=(double*)list2.bptr;
00329   zdotbar=(double*)(list2.bptr+nvar*sizeof(double));
00330   list2.restoreposition(); // save pointer to beginning of record;
00331 
00332   double * x_tilde=px->get_u_tilde();
00333   double * x_dot_tilde=px->get_u_dot_tilde();
00334   double * x_bar_tilde=px->get_u_bar_tilde();
00335   double * x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00336   double * y_tilde=py->get_u_tilde();
00337   double * y_dot_tilde=py->get_u_dot_tilde();
00338   double * y_bar_tilde=py->get_u_bar_tilde();
00339   double * y_dot_bar_tilde=py->get_u_dot_bar_tilde();
00340   double * z_bar_tilde=pz->get_u_bar_tilde();
00341   double * z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00342   // Do second "reverse-reverse" pass calculations
00343 
00344   for (size_t i=0;i<nvar;i++)
00345   {
00346     z_bar_tilde[i]=0;
00347     z_dot_bar_tilde[i]=0;
00348   }
00349 
00350   // start with x and add y
00351   for (size_t i=0;i<nvar;i++)
00352   {
00353     //*x_tilde+=(pf->d2f11)(xu,yu)*zbar[i]*x_bar_tilde[i];
00354     //z_bar_tilde[i]+=(pf->df1)(xu,yu)*x_bar_tilde[i];
00355     //*y_tilde+=(pf->d2f12)(xu,yu)*zbar[i]*x_bar_tilde[i];
00356     z_bar_tilde[i]+=yu*x_bar_tilde[i];
00357     *y_tilde+=zbar[i]*x_bar_tilde[i];
00358   }
00359 
00360   for (size_t i=0;i<nvar;i++)
00361   {
00362     //*x_tilde+=(pf->d2f11)(xu,yu)*zdotbar[i]*x_dot_bar_tilde[i];
00363     //*y_tilde+=(pf->d2f12)(xu,yu)*zdotbar[i]*x_dot_bar_tilde[i];
00364     //z_dot_bar_tilde[i]+=(pf->df1)(xu,yu)*x_dot_bar_tilde[i];
00365     *y_tilde+=zdotbar[i]*x_dot_bar_tilde[i];
00366     z_dot_bar_tilde[i]+=yu*x_dot_bar_tilde[i];
00367   }
00368 
00369   /*
00370   for (i=0;i<nvar;i++)
00371   {
00372     x_dot_tilde[i]+=(pf->d2f11)(xu,yu)*zdotbar[i]*x_bar_tilde[i];
00373     z_dot_bar_tilde[i]+=(pf->d2f11)(xu,yu)*xdot[i]*x_bar_tilde[i];
00374     *x_tilde+=(pf->d3f111)(xu,yu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00375     *y_tilde+=(pf->d3f112)(xu,yu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00376   }
00377   */
00378   // start with y and add x
00379   for (size_t i=0;i<nvar;i++)
00380   {
00381     //*y_tilde+=(pf->d2f22)(xu,yu)*zbar[i]*y_bar_tilde[i];
00382     //*x_tilde+=(pf->d2f12)(xu,yu)*zbar[i]*y_bar_tilde[i];
00383     //z_bar_tilde[i]+=(pf->df2)(xu,yu)*y_bar_tilde[i];
00384     *x_tilde+=zbar[i]*y_bar_tilde[i];
00385     z_bar_tilde[i]+=xu*y_bar_tilde[i];
00386   }
00387 
00388   for (size_t i=0;i<nvar;i++)
00389   {
00390     //*y_tilde+=(pf->d2f22)(xu,yu)*zdotbar[i]*y_dot_bar_tilde[i];
00391     //*x_tilde+=(pf->d2f12)(xu,yu)*zdotbar[i]*y_dot_bar_tilde[i];
00392     //z_dot_bar_tilde[i]+=(pf->df2)(xu,yu)*y_dot_bar_tilde[i];
00393     *x_tilde+=zdotbar[i]*y_dot_bar_tilde[i];
00394     z_dot_bar_tilde[i]+=xu*y_dot_bar_tilde[i];
00395   }
00396 
00397   /*
00398   for (i=0;i<nvar;i++)
00399   {
00400     y_dot_tilde[i]+=(pf->d2f22)(xu,yu)*zdotbar[i]*y_bar_tilde[i];
00401     z_dot_bar_tilde[i]+=(pf->d2f22)(xu,yu)*ydot[i]*y_bar_tilde[i];
00402     *y_tilde+=(pf->d3f222)(xu,yu)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00403     *x_tilde+=(pf->d3f122)(xu,yu)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00404   }
00405   */
00406   for (size_t i=0;i<nvar;i++)
00407   {
00408     //*x_tilde+=(pf->d3f112)(xu,yu)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00409     //*y_tilde+=(pf->d3f122)(xu,yu)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00410     //y_dot_tilde[i]+=(pf->d2f12)(xu,yu)*zdotbar[i]*x_bar_tilde[i];
00411     //z_dot_bar_tilde[i]+=(pf->d2f12)(xu,yu)*ydot[i]*x_bar_tilde[i];
00412     y_dot_tilde[i]+=zdotbar[i]*x_bar_tilde[i];
00413     z_dot_bar_tilde[i]+=ydot[i]*x_bar_tilde[i];
00414   }
00415   for (size_t i=0;i<nvar;i++)
00416   {
00417     //*x_tilde+=(pf->d3f112)(xu,yu)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00418     //*y_tilde+=(pf->d3f122)(xu,yu)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00419     //x_dot_tilde[i]+=(pf->d2f12)(xu,yu)*zdotbar[i]*y_bar_tilde[i];
00420     //z_dot_bar_tilde[i]+=(pf->d2f12)(xu,yu)*xdot[i]*y_bar_tilde[i];
00421     x_dot_tilde[i]+=zdotbar[i]*y_bar_tilde[i];
00422     z_dot_bar_tilde[i]+=xdot[i]*y_bar_tilde[i];
00423   }
00424 }
00425 
00430 void read_pass2_3_prod(void)
00431 {
00432   // We are going backword for bptr and forward for bptr2
00433   // the current entry+2 in bptr is the size of the record i.e
00434   // points to the next record
00435   unsigned int nvar=df1b2variable::nvar;
00436   fixed_smartlist & nlist=f1b2gradlist->nlist;
00437   test_smartlist& list=f1b2gradlist->list;
00438    // nlist-=sizeof(int);
00439   // get record size
00440   int num_bytes=nlist.bptr->numbytes;
00441   // backup the size of the record
00442   list-=num_bytes;
00443   list.saveposition(); // save pointer to beginning of record;
00444   // save the pointer to the beginning of the record
00445   double xu;
00446   double yu;
00447   double * xdot;
00448   double * ydot;
00449   //df1b2_header x,z;
00450   //df1b2function2 * pf;
00451 
00452   // get info from tape1
00453   // get info from tape1
00454 #if defined(SAFE_ALL)
00455   checkidentiferstring("DL",list);
00456 #endif
00457   df1b2_header * px=(df1b2_header *) list.bptr;
00458   list.bptr+=sizeof(df1b2_header);
00459   df1b2_header * py=(df1b2_header *) list.bptr;
00460   list.bptr+=sizeof(df1b2_header);
00461   df1b2_header * pz=(df1b2_header *) list.bptr;
00462   list.bptr+=sizeof(df1b2_header);
00463   //pf=*(df1b2function2 **) list.bptr;
00464   //list.bptr+=sizeof(char*);
00465   memcpy(&xu,list.bptr,sizeof(double));
00466   list.bptr+=sizeof(double);
00467   memcpy(&yu,list.bptr,sizeof(double));
00468   list.bptr+=sizeof(double);
00469   xdot=(double*)list.bptr;
00470   list.bptr+=nvar*sizeof(double);
00471   ydot=(double*)list.bptr;
00472   list.restoreposition(); // save pointer to beginning of record;
00473 
00474   //*(px->u_tilde)+=(pf->df1)(xu,yu)* *(pz->u_tilde);
00475   //*(py->u_tilde)+=(pf->df2)(xu,yu)* *(pz->u_tilde);
00476   *(px->u_tilde)+=yu* *(pz->u_tilde);
00477   *(py->u_tilde)+=xu* *(pz->u_tilde);
00478   for (size_t i=0;i<nvar;i++)
00479   {
00480     //*(px->u_tilde)+=(pf->d2f11)(xu,yu)*xdot[i]*pz->u_dot_tilde[i];
00481     //*(py->u_tilde)+=(pf->d2f12)(xu,yu)*xdot[i]*pz->u_dot_tilde[i];
00482     *(py->u_tilde)+=xdot[i]*pz->u_dot_tilde[i];
00483     //*(py->u_tilde)+=(pf->d2f22)(xu,yu)*ydot[i]*pz->u_dot_tilde[i];
00484     //*(px->u_tilde)+=(pf->d2f12)(xu,yu)*ydot[i]*pz->u_dot_tilde[i];
00485     *(px->u_tilde)+=ydot[i]*pz->u_dot_tilde[i];
00486   }
00487   for (size_t i=0;i<nvar;i++)
00488   {
00489     //px->u_dot_tilde[i]+=(pf->df1)(xu,yu)*pz->u_dot_tilde[i];
00490     //py->u_dot_tilde[i]+=(pf->df2)(xu,yu)*pz->u_dot_tilde[i];
00491     px->u_dot_tilde[i]+=yu*pz->u_dot_tilde[i];
00492     py->u_dot_tilde[i]+=xu*pz->u_dot_tilde[i];
00493   }
00494   *(pz->u_tilde)=0;
00495   for (size_t i=0;i<nvar;i++)
00496   {
00497     pz->u_dot_tilde[i]=0;
00498   }
00499 }