ADMB Documentation  11.1x.2730
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2f12.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2f12.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 
00017 void ad_read_pass2(void);
00018 
00023  int df1b2_gradlist::write_pass1(const df1b2variable * _px,
00024    const df1b2variable * _py,df1b2variable * pz,df1b2function2 * pf)
00025  {
00026    ADUNCONST(df1b2variable*,px)
00027    ADUNCONST(df1b2variable*,py)
00028    ncount++;
00029 #if defined(CHECK_COUNT)
00030   if (ncount >= ncount_check)
00031     ncount_checker(ncount,ncount_check);
00032 #endif
00033 
00034   unsigned int nvar=df1b2variable::nvar;
00035 
00036    size_t total_bytes=3*sizeof(df1b2_header)+sizeof(char*)
00037      +2*(nvar+1)*sizeof(double);
00038 // string identifier debug stuff
00039 #if defined(SAFE_ALL)
00040   char ids[]="BY";
00041   int slen=strlen(ids);
00042   total_bytes+=slen;
00043 #endif
00044   list.check_buffer_size(total_bytes);
00045   void * tmpptr=list.bptr;
00046 #if defined(SAFE_ALL)
00047   memcpy(list,ids,slen);
00048 #endif
00049 // end of string identifier debug stuff
00050 
00051    memcpy(list,(df1b2_header*)(px),sizeof(df1b2_header));
00052    memcpy(list,(df1b2_header*)(py),sizeof(df1b2_header));
00053    memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00054    memcpy(list,&pf,sizeof(char *));
00055    //*(char**)(list.bptr)=(char*)pf;
00056    memcpy(list,px->get_u(),sizeof(double));
00057    memcpy(list,py->get_u(),sizeof(double));
00058    memcpy(list,px->get_u_dot(),nvar*sizeof(double));
00059    memcpy(list,py->get_u_dot(),nvar*sizeof(double));
00060    // ***** write  record size
00061    nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00062    nlist.bptr->pf=(ADrfptr)(&ad_read_pass2);
00063    ++nlist;
00064    return 0;
00065  }
00066 
00067 
00068 void read_pass2_1(void);
00069 void read_pass2_2(void);
00070 void read_pass2_3(void);
00071 
00076 void ad_read_pass2(void)
00077 {
00078   switch(df1b2variable::passnumber)
00079   {
00080   case 1:
00081     read_pass2_1();
00082     break;
00083   case 2:
00084     read_pass2_2();
00085     break;
00086   case 3:
00087     read_pass2_3();
00088     break;
00089   default:
00090     cerr << "illegal value for df1b2variable::pass = "
00091          << df1b2variable::passnumber << endl;
00092     exit(1);
00093   }
00094 }
00095 
00100 void read_pass2_1(void)
00101 {
00102   // We are going backword for bptr and nbptr
00103   // and  forward for bptr2 and nbptr2
00104   // the current entry+2 in bptr is the size of the record i.e
00105   // points to the next record
00106   //char * bptr=f1b2gradlist->bptr;
00107   //char * bptr2=f1b2gradlist2->bptr;
00108 
00109   unsigned int nvar=df1b2variable::nvar;
00110 
00111   test_smartlist& list=f1b2gradlist->list;
00112   //f1b2gradlist->nlist-=sizeof(int);
00113   int num_bytes=f1b2gradlist->nlist.bptr->numbytes;
00114   list-=num_bytes;
00115   list.saveposition(); // save pointer to beginning of record;
00116   double xu,yu;
00117   //ad_dstar xdot,ydot;
00118   df1b2function2 * pf;
00119 
00120   // get info from tape1
00121 #if defined(SAFE_ARRAYS)
00122   checkidentiferstring("BY",f1b2gradlist->list);
00123 #endif
00124   char * bptr=f1b2gradlist->list.bptr;
00125   df1b2_header * px=(df1b2_header *) bptr;
00126   bptr+=sizeof(df1b2_header);
00127   df1b2_header * py=(df1b2_header *) bptr;
00128   bptr+=sizeof(df1b2_header);
00129   df1b2_header * pz=(df1b2_header *) bptr;
00130   bptr+=sizeof(df1b2_header);
00131   pf=*(df1b2function2 **) bptr;
00132   bptr+=sizeof(char*);
00133   memcpy(&xu,bptr,sizeof(double));
00134   bptr+=sizeof(double);
00135   memcpy(&yu,bptr,sizeof(double));
00136   bptr+=sizeof(double);
00137   double * xdot=(double*)bptr;
00138   bptr+=nvar*sizeof(double);
00139   double * ydot=(double*)bptr;
00140 
00141   list.restoreposition(); // save pointer to beginning of record;
00142 
00143   // ****************************************************************
00144   // turn this off if no third derivatives are calculated
00145   // if (!no_third_derivatives)
00146   // {
00147   // save for second reverse pass
00148   // save identifier 1
00149      test_smartlist & list2 = f1b2gradlist->list2;
00150 
00151 
00152    size_t total_bytes=2*nvar*sizeof(double);
00153 // string identifier debug stuff
00154 #if defined(SAFE_ALL)
00155   char ids[]="FW";
00156   int slen=strlen(ids);
00157   total_bytes+=slen;
00158 #endif
00159   list2.check_buffer_size(total_bytes);
00160   void * tmpptr=list2.bptr;
00161 #if defined(SAFE_ALL)
00162   memcpy(list2,ids,slen);
00163 #endif
00164 
00165      fixed_smartlist2 & nlist2 = f1b2gradlist->nlist2;
00166      memcpy(list2,pz->get_u_bar(),nvar*sizeof(double));
00167      memcpy(list2,pz->get_u_dot_bar(),nvar*sizeof(double));
00168      *nlist2.bptr=adptr_diff(list2.bptr,tmpptr);
00169      ++nlist2;
00170   // }
00171   //
00172   // ****************************************************************
00173 #if defined(PRINT_DERS)
00174  print_derivatives(pf->funname,(pf->f)(xu,yu),(pf->df1)(xu,yu),
00175   (pf->df2)(xu,yu),(pf->d2f11)(xu,yu),(pf->d2f12)(xu,yu),(pf->d2f22)(xu,yu),
00176   (pf->d3f111)(xu,yu),(pf->d3f112)(xu,yu),(pf->d3f122)(xu,yu),
00177   (pf->d3f222)(xu,yu),1);
00178  print_derivatives(pz,"z");
00179  print_derivatives(px,"x");
00180  print_derivatives(py,"y");
00181 #endif
00182 #if defined(__DERCHECK__)
00183   if (derchecker)
00184   if (derchecker->node_number)
00185   {
00186     if (derchecker->counter == derchecker->node_number)
00187     {
00188       switch (derchecker->pass_number) // increment the variable of interest
00189       {
00190       case 2:
00191         switch(derchecker->vartype)
00192         {
00193         case 1:
00194           if (!derchecker->dotflag)
00195             px->u_bar[derchecker->index-1]+=derchecker->delta;
00196           else
00197             px->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00198           break;
00199         case 2:
00200           if (!derchecker->dotflag)
00201             py->u_bar[derchecker->index-1]+=derchecker->delta;
00202           else
00203             py->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00204           break;
00205         case 3:
00206           if (!derchecker->dotflag)
00207             pz->u_bar[derchecker->index-1]+=derchecker->delta;
00208           else
00209             pz->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00210           break;
00211         default:
00212           cerr << "Invalid index value for dercheck_index was "
00213                << derchecker->index << endl;
00214           break;
00215         }
00216         break;
00217       case 3:
00218         switch(derchecker->vartype)
00219         {
00220         case 1:
00221           if (!derchecker->dotflag)
00222             px->u_bar[derchecker->index-1]-=derchecker->delta;
00223           else
00224             px->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00225           break;
00226         case 2:
00227           if (!derchecker->dotflag)
00228             py->u_bar[derchecker->index-1]-=derchecker->delta;
00229           else
00230             py->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00231           break;
00232         case 3:
00233           if (!derchecker->dotflag)
00234             pz->u_bar[derchecker->index-1]-=derchecker->delta;
00235           else
00236             pz->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00237           break;
00238         default:
00239           cerr << "Invalid index value for dercheck_index was "
00240                << derchecker->index << endl;
00241           break;
00242         }
00243         break;
00244       }
00245     }
00246   }
00247 #endif
00248 
00249   // Do first reverse pass calculations
00250   for (size_t i=0;i<nvar;i++)
00251   {
00252     //double df1=(pf->df1)(xu,yu);
00253     px->u_bar[i]+=(pf->df1)(xu,yu)*pz->u_bar[i];
00254   }
00255   for (size_t i=0;i<nvar;i++)
00256   {
00257     py->u_bar[i]+=(pf->df2)(xu,yu)*pz->u_bar[i];
00258   }
00259   for (size_t i=0;i<nvar;i++)
00260   {
00261     px->u_bar[i]+=(pf->d2f11)(xu,yu)*xdot[i]*pz->u_dot_bar[i];
00262     px->u_bar[i]+=(pf->d2f12)(xu,yu)*ydot[i]*pz->u_dot_bar[i];
00263 #if defined(ADDEBUG_PRINT)
00264     cout << px->u_bar[i] << " " << pz->u_dot_bar[i] << " " << addebug_count
00265          << endl;
00266 #endif
00267   }
00268   for (size_t i=0;i<nvar;i++)
00269   {
00270     //py->u_bar[i]+=(pf->d2f22)(*(px->u),*(py->u))*ydot[i]*pz->u_dot_bar[i];
00271     //py->u_bar[i]+=(pf->d2f12)(*(px->u),*(py->u))*xdot[i]*pz->u_dot_bar[i];
00272     py->u_bar[i]+=(pf->d2f22)(xu,yu)*ydot[i]*pz->u_dot_bar[i];
00273     py->u_bar[i]+=(pf->d2f12)(xu,yu)*xdot[i]*pz->u_dot_bar[i];
00274 #if defined(ADDEBUG_PRINT)
00275     cout << py->u_bar[i] << " " << pz->u_dot_bar[i] << " " << addebug_count
00276          << endl;
00277 #endif
00278   }
00279   for (size_t i=0;i<nvar;i++)
00280   {
00281     //px->u_dot_bar[i]+=(pf->df1)(*(px->u),*(py->u))*pz->u_dot_bar[i];
00282     px->u_dot_bar[i]+=(pf->df1)(xu,yu)*pz->u_dot_bar[i];
00283 #if defined(ADDEBUG_PRINT)
00284     cout << px->u_dot_bar[i] << " " << addebug_count << endl;
00285     cout << pz->u_dot_bar[i] << " " << addebug_count << endl;
00286 #endif
00287   }
00288   for (size_t i=0;i<nvar;i++)
00289   {
00290     py->u_dot_bar[i]+=(pf->df2)(xu,yu)*pz->u_dot_bar[i];
00291 #if defined(ADDEBUG_PRINT)
00292     cout << py->u_dot_bar[i] << " " << addebug_count << endl;
00293     cout << pz->u_dot_bar[i] << " " << addebug_count << endl;
00294 #endif
00295   }
00296 
00297   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00298   for (size_t i=0;i<nvar;i++)
00299   {
00300     pz->u_bar[i]=0;
00301   }
00302   for (size_t i=0;i<nvar;i++)
00303   {
00304     pz->u_dot_bar[i]=0;
00305   }
00306 
00307 #if defined(PRINT_DERS)
00308  print_derivatives(pz,"z");
00309  print_derivatives(px,"x");
00310  print_derivatives(py,"y");
00311 #endif
00312 }
00313 
00318 void read_pass2_2(void)
00319 {
00320   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00321   // We are going forward for bptr and backword for bptr2
00322   //
00323   // list 1
00324   //
00325 
00326   unsigned int nvar=df1b2variable::nvar;
00327 
00328   test_smartlist & list=f1b2gradlist->list;
00329 
00330   size_t total_bytes=3*sizeof(df1b2_header)+sizeof(char*)
00331     +2*(nvar+1)*sizeof(double);
00332 // string identifier debug stuff
00333 #if defined(SAFE_ALL)
00334   char ids[]="BY";
00335   int slen=strlen(ids);
00336   total_bytes+=slen;
00337 #endif
00338   list.check_buffer_size(total_bytes);
00339 // end of string identifier debug stuff
00340 
00341   list.saveposition(); // save pointer to beginning of record;
00342   fixed_smartlist & nlist=f1b2gradlist->nlist;
00343    // nlist-=sizeof(int);
00344   // get record size
00345   int num_bytes=nlist.bptr->numbytes;
00346     // nlist+=nlist_record_size;
00347   //
00348   // list 2
00349   //
00350   test_smartlist & list2=f1b2gradlist->list2;
00351   fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00352   // get record size
00353   int num_bytes2=*nlist2.bptr;
00354   --nlist2;
00355   // backup the size of the record
00356   list2-=num_bytes2;
00357   list2.saveposition(); // save pointer to beginning of record;
00358   // save the pointer to the beginning of the record
00359   // bptr and bptr2 now both point to the beginning of their records
00360 
00361   double xu,yu;
00362   double * xdot;
00363   double * ydot;
00364   //df1b2_header x,z;
00365   df1b2function2 * pf;
00366 
00367   // get info from tape1
00368   // get info from tape1
00369 #if defined(SAFE_ARRAYS)
00370   checkidentiferstring("BY",list);
00371   checkidentiferstring("FW",list2);
00372 #endif
00373   df1b2_header * px=(df1b2_header *) list.bptr;
00374   list.bptr+=sizeof(df1b2_header);
00375   df1b2_header * py=(df1b2_header *) list.bptr;
00376   list.bptr+=sizeof(df1b2_header);
00377   df1b2_header * pz=(df1b2_header *) list.bptr;
00378   list.bptr+=sizeof(df1b2_header);
00379   pf=*(df1b2function2 **) list.bptr;
00380   list.bptr+=sizeof(char*);
00381   memcpy(&xu,list.bptr,sizeof(double));
00382   list.bptr+=sizeof(double);
00383   memcpy(&yu,list.bptr,sizeof(double));
00384   list.bptr+=sizeof(double);
00385   xdot=(double*)list.bptr;
00386   list.bptr+=nvar*sizeof(double);
00387   ydot=(double*)list.bptr;
00388   list.restoreposition(num_bytes); // save pointer to beginning of record;
00389 
00390   double * zbar;
00391   double * zdotbar;
00392 
00393   zbar=(double*)list2.bptr;
00394   zdotbar=(double*)(list2.bptr+nvar*sizeof(double));
00395   list2.restoreposition(); // save pointer to beginning of record;
00396 
00397   double * x_tilde=px->get_u_tilde();
00398   double * x_dot_tilde=px->get_u_dot_tilde();
00399   double * x_bar_tilde=px->get_u_bar_tilde();
00400   double * x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00401   double * y_tilde=py->get_u_tilde();
00402   double * y_dot_tilde=py->get_u_dot_tilde();
00403   double * y_bar_tilde=py->get_u_bar_tilde();
00404   double * y_dot_bar_tilde=py->get_u_dot_bar_tilde();
00405   double * z_bar_tilde=pz->get_u_bar_tilde();
00406   double * z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00407   // Do second "reverse-reverse" pass calculations
00408 #if defined(PRINT_DERS)
00409  print_derivatives(pf->funname,(pf->f)(xu,yu),(pf->df1)(xu,yu),
00410   (pf->df2)(xu,yu),(pf->d2f11)(xu,yu),(pf->d2f12)(xu,yu),(pf->d2f22)(xu,yu),
00411   (pf->d3f111)(xu,yu),(pf->d3f112)(xu,yu),(pf->d3f122)(xu,yu),
00412   (pf->d3f222)(xu,yu),1);
00413  print_derivatives(pz,"z");
00414  print_derivatives(px,"x");
00415  print_derivatives(py,"y");
00416 #endif
00417 
00418   for (size_t i=0;i<nvar;i++)
00419   {
00420     z_bar_tilde[i]=0;
00421     z_dot_bar_tilde[i]=0;
00422   }
00423 
00424   // start with x and add y
00425   for (size_t i=0;i<nvar;i++)
00426   {
00427     *x_tilde+=(pf->d2f11)(xu,yu)*zbar[i]*x_bar_tilde[i];
00428     z_bar_tilde[i]+=(pf->df1)(xu,yu)*x_bar_tilde[i];
00429     *y_tilde+=(pf->d2f12)(xu,yu)*zbar[i]*x_bar_tilde[i];
00430   }
00431 
00432   for (size_t i=0;i<nvar;i++)
00433   {
00434     *x_tilde+=(pf->d2f11)(xu,yu)*zdotbar[i]*x_dot_bar_tilde[i];
00435     *y_tilde+=(pf->d2f12)(xu,yu)*zdotbar[i]*x_dot_bar_tilde[i];
00436     z_dot_bar_tilde[i]+=(pf->df1)(xu,yu)*x_dot_bar_tilde[i];
00437   }
00438 
00439   for (size_t i=0;i<nvar;i++)
00440   {
00441     x_dot_tilde[i]+=(pf->d2f11)(xu,yu)*zdotbar[i]*x_bar_tilde[i];
00442     z_dot_bar_tilde[i]+=(pf->d2f11)(xu,yu)*xdot[i]*x_bar_tilde[i];
00443     *x_tilde+=(pf->d3f111)(xu,yu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00444     *y_tilde+=(pf->d3f112)(xu,yu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00445   }
00446   // start with y and add x
00447   for (size_t i=0;i<nvar;i++)
00448   {
00449     *y_tilde+=(pf->d2f22)(xu,yu)*zbar[i]*y_bar_tilde[i];
00450     *x_tilde+=(pf->d2f12)(xu,yu)*zbar[i]*y_bar_tilde[i];
00451     z_bar_tilde[i]+=(pf->df2)(xu,yu)*y_bar_tilde[i];
00452   }
00453 
00454   for (size_t i=0;i<nvar;i++)
00455   {
00456     *y_tilde+=(pf->d2f22)(xu,yu)*zdotbar[i]*y_dot_bar_tilde[i];
00457     *x_tilde+=(pf->d2f12)(xu,yu)*zdotbar[i]*y_dot_bar_tilde[i];
00458     z_dot_bar_tilde[i]+=(pf->df2)(xu,yu)*y_dot_bar_tilde[i];
00459   }
00460 
00461   for (size_t i=0;i<nvar;i++)
00462   {
00463     y_dot_tilde[i]+=(pf->d2f22)(xu,yu)*zdotbar[i]*y_bar_tilde[i];
00464     z_dot_bar_tilde[i]+=(pf->d2f22)(xu,yu)*ydot[i]*y_bar_tilde[i];
00465     *y_tilde+=(pf->d3f222)(xu,yu)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00466     *x_tilde+=(pf->d3f122)(xu,yu)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00467   }
00468   for (size_t i=0;i<nvar;i++)
00469   {
00470     *x_tilde+=(pf->d3f112)(xu,yu)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00471     *y_tilde+=(pf->d3f122)(xu,yu)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00472     y_dot_tilde[i]+=(pf->d2f12)(xu,yu)*zdotbar[i]*x_bar_tilde[i];
00473     z_dot_bar_tilde[i]+=(pf->d2f12)(xu,yu)*ydot[i]*x_bar_tilde[i];
00474   }
00475   for (size_t i=0;i<nvar;i++)
00476   {
00477     *x_tilde+=(pf->d3f112)(xu,yu)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00478     *y_tilde+=(pf->d3f122)(xu,yu)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00479     x_dot_tilde[i]+=(pf->d2f12)(xu,yu)*zdotbar[i]*y_bar_tilde[i];
00480     z_dot_bar_tilde[i]+=(pf->d2f12)(xu,yu)*xdot[i]*y_bar_tilde[i];
00481   }
00482 #if defined(__DERCHECK__)
00483   if (derchecker->node_number)
00484   {
00485     if (derchecker->counter == derchecker->node_number)
00486     {
00487       if (derchecker->pass_number==1) // increment the variable of interest
00488       {
00489         switch(derchecker->vartype)
00490         {
00491         case 1:
00492           if (!derchecker->dotflag)
00493             derchecker->der_value=
00494               px->u_bar_tilde[derchecker->index-1];
00495           else
00496             derchecker->der_value=
00497               px->u_dot_bar_tilde[derchecker->index-1];
00498           break;
00499         case 2:
00500           if (!derchecker->dotflag)
00501             derchecker->der_value=
00502               py->u_bar_tilde[derchecker->index-1];
00503           else
00504             derchecker->der_value=
00505               py->u_dot_bar_tilde[derchecker->index-1];
00506           break;
00507         case 3:
00508           if (!derchecker->dotflag)
00509             derchecker->der_value=
00510               pz->u_bar_tilde[derchecker->index-1];
00511           else
00512             derchecker->der_value=
00513               pz->u_dot_bar_tilde[derchecker->index-1];
00514           break;
00515         default:
00516           cerr << "Invalid index value for dercheck_index was "
00517                << derchecker->index << endl;
00518         }
00519       }
00520     }
00521   }
00522 #endif
00523 #if defined(PRINT_DERS)
00524  print_derivatives(pz,"z");
00525  print_derivatives(px,"x");
00526  print_derivatives(py,"y");
00527 #endif
00528 }
00529 
00534 void read_pass2_3(void)
00535 {
00536   // We are going backword for bptr and forward for bptr2
00537   // the current entry+2 in bptr is the size of the record i.e
00538   // points to the next record
00539 
00540   unsigned int nvar=df1b2variable::nvar;
00541 
00542   fixed_smartlist & nlist=f1b2gradlist->nlist;
00543   test_smartlist& list=f1b2gradlist->list;
00544    // nlist-=sizeof(int);
00545   // get record size
00546   int num_bytes=nlist.bptr->numbytes;
00547   // backup the size of the record
00548   list-=num_bytes;
00549   list.saveposition(); // save pointer to beginning of record;
00550   // save the pointer to the beginning of the record
00551   double xu;
00552   double yu;
00553   double * xdot;
00554   double * ydot;
00555   //df1b2_header x,z;
00556   df1b2function2 * pf;
00557 
00558   // get info from tape1
00559   // get info from tape1
00560 #if defined(SAFE_ARRAYS)
00561   checkidentiferstring("BY",list);
00562 #endif
00563   df1b2_header * px=(df1b2_header *) list.bptr;
00564   list.bptr+=sizeof(df1b2_header);
00565   df1b2_header * py=(df1b2_header *) list.bptr;
00566   list.bptr+=sizeof(df1b2_header);
00567   df1b2_header * pz=(df1b2_header *) list.bptr;
00568   list.bptr+=sizeof(df1b2_header);
00569   pf=*(df1b2function2 **) list.bptr;
00570   list.bptr+=sizeof(char*);
00571   memcpy(&xu,list.bptr,sizeof(double));
00572   list.bptr+=sizeof(double);
00573   memcpy(&yu,list.bptr,sizeof(double));
00574   list.bptr+=sizeof(double);
00575   xdot=(double*)list.bptr;
00576   list.bptr+=nvar*sizeof(double);
00577   ydot=(double*)list.bptr;
00578   list.restoreposition(); // save pointer to beginning of record;
00579 #if defined(PRINT_DERS)
00580  print_derivatives(pf->funname,(pf->f)(xu,yu),(pf->df1)(xu,yu),
00581   (pf->df2)(xu,yu),(pf->d2f11)(xu,yu),(pf->d2f12)(xu,yu),(pf->d2f22)(xu,yu),
00582   (pf->d3f111)(xu,yu),(pf->d3f112)(xu,yu),(pf->d3f122)(xu,yu),
00583   (pf->d3f222)(xu,yu),1);
00584  print_derivatives(pz,"z");
00585  print_derivatives(px,"x");
00586  print_derivatives(py,"y");
00587 #endif
00588 
00589   *(px->u_tilde)+=(pf->df1)(xu,yu)* *(pz->u_tilde);
00590   *(py->u_tilde)+=(pf->df2)(xu,yu)* *(pz->u_tilde);
00591   for (size_t i=0;i<nvar;i++)
00592   {
00593     *(px->u_tilde)+=(pf->d2f11)(xu,yu)*xdot[i]*pz->u_dot_tilde[i];
00594     *(py->u_tilde)+=(pf->d2f12)(xu,yu)*xdot[i]*pz->u_dot_tilde[i];
00595     *(py->u_tilde)+=(pf->d2f22)(xu,yu)*ydot[i]*pz->u_dot_tilde[i];
00596     *(px->u_tilde)+=(pf->d2f12)(xu,yu)*ydot[i]*pz->u_dot_tilde[i];
00597   }
00598   for (size_t i=0;i<nvar;i++)
00599   {
00600     px->u_dot_tilde[i]+=(pf->df1)(xu,yu)*pz->u_dot_tilde[i];
00601     py->u_dot_tilde[i]+=(pf->df2)(xu,yu)*pz->u_dot_tilde[i];
00602   }
00603   *(pz->u_tilde)=0;
00604   for (size_t i=0;i<nvar;i++)
00605   {
00606     pz->u_dot_tilde[i]=0;
00607   }
00608 #if defined(PRINT_DERS)
00609  print_derivatives(pz,"z");
00610  print_derivatives(px,"x");
00611  print_derivatives(py,"y");
00612 #endif
00613 }