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