ADMB Documentation  11.1.2397
 All Classes Files Functions Variables Typedefs Friends Defines
df33fun1.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df33fun1.cpp 2370 2014-09-17 21:27:31Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 #include "df33fun.h"
00013 
00014 #ifndef OPT_LIB
00015   #include <cassert>
00016   #include <climits>
00017 #endif
00018 
00019 void ad_read_pass2_dvdvdv(void);
00020 
00025  int df1b2_gradlist::write_pass1(const df1b2variable * _px,
00026    const df1b2variable * _py,const df1b2variable * pw,
00027    const df1b2variable * pz,
00028    double df_x, double df_y, double df_z,
00029    double df_xx, double df_xy,double df_xz,double df_yy,
00030    double df_yz,double df_zz,
00031    double df_xxx,
00032    double df_xxy,
00033    double df_xxz,
00034    double df_xyy,
00035    double df_xyz,
00036    double df_xzz,
00037    double df_yyy,
00038    double df_yyz,
00039    double df_yzz,
00040    double df_zzz)
00041  {
00042    ADUNCONST(df1b2variable*,px)
00043    ADUNCONST(df1b2variable*,py)
00044    ncount++;
00045 #if defined(CHECK_COUNT)
00046   if (ncount >= ncount_check)
00047     ncount_checker(ncount,ncount_check);
00048 #endif
00049    int nvar=df1b2variable::nvar;
00050 
00051    size_t total_bytes=4*sizeof(df1b2_header)+sizeof(char*)
00052      +(3*nvar+22)*sizeof(double);
00053 
00054 // string identifier debug stuff
00055 #if defined(SAFE_ALL)
00056   char ids[]="U8";
00057   int slen=strlen(ids);
00058   total_bytes+=slen;
00059 #endif
00060 
00061 #ifndef OPT_LIB
00062   assert(total_bytes <= INT_MAX);
00063 #endif
00064 
00065   list.check_buffer_size((int)total_bytes);
00066 
00067   void * tmpptr=list.bptr;
00068 #if defined(SAFE_ALL)
00069   memcpy(list,ids,slen);
00070 #endif
00071 // end of string identifier debug stuff
00072 
00073    memcpy(list,(df1b2_header*)(px),sizeof(df1b2_header));
00074    memcpy(list,(df1b2_header*)(py),sizeof(df1b2_header));
00075    memcpy(list,(df1b2_header*)(pw),sizeof(df1b2_header));
00076    memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00077    //memcpy(list,&pf,sizeof(char *));
00078    //*(char**)(list.bptr)=(char*)pf;
00079 
00080   /*
00081    memcpy(list,&df_x,sizeof(double));
00082    memcpy(list,&df_y,sizeof(double));
00083    memcpy(list,&df_xx,sizeof(double));
00084    memcpy(list,&df_xy,sizeof(double));
00085    memcpy(list,&df_yy,sizeof(double));
00086    memcpy(list,&df_xxx,sizeof(double));
00087    memcpy(list,&df_xxy,sizeof(double));
00088    memcpy(list,&df_xyy,sizeof(double));
00089    memcpy(list,&df_yyy,sizeof(double));
00090   */
00091    const int sizeofdouble = sizeof(double);
00092    memcpy(list,&df_x,sizeofdouble);
00093    memcpy(list,&df_y,sizeofdouble);
00094    memcpy(list,&df_z,sizeofdouble);
00095    memcpy(list,&df_xx,sizeofdouble);
00096    memcpy(list,&df_xy,sizeofdouble);
00097    memcpy(list,&df_xz,sizeofdouble);
00098    memcpy(list,&df_yy,sizeofdouble);
00099    memcpy(list,&df_yz,sizeofdouble);
00100    memcpy(list,&df_zz,sizeofdouble);
00101    memcpy(list,&df_xxx,sizeofdouble);
00102    memcpy(list,&df_xxy,sizeofdouble);
00103    memcpy(list,&df_xxz,sizeofdouble);
00104    memcpy(list,&df_xyy,sizeofdouble);
00105    memcpy(list,&df_xyz,sizeofdouble);
00106    memcpy(list,&df_xzz,sizeofdouble);
00107    memcpy(list,&df_yyy,sizeofdouble);
00108    memcpy(list,&df_yyz,sizeofdouble);
00109    memcpy(list,&df_yzz,sizeofdouble);
00110    memcpy(list,&df_zzz,sizeofdouble);
00111 
00112    memcpy(list,px->get_u(),sizeofdouble);
00113    memcpy(list,py->get_u(),sizeofdouble);
00114    memcpy(list,pw->get_u(),sizeofdouble);
00115    memcpy(list,px->get_u_dot(),nvar*sizeofdouble);
00116    memcpy(list,py->get_u_dot(),nvar*sizeofdouble);
00117    memcpy(list,pw->get_u_dot(),nvar*sizeofdouble);
00118    // ***** write  record size
00119    nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00120    nlist.bptr->pf=(ADrfptr)(&ad_read_pass2_dvdvdv);
00121    ++nlist;
00122    return 0;
00123  }
00124 
00125 
00126 void read_pass2_1_dvdvdv(void);
00127 void read_pass2_2_dvdvdv(void);
00128 void read_pass2_3_dvdvdv(void);
00129 
00134 void ad_read_pass2_dvdvdv(void)
00135 {
00136   switch(df1b2variable::passnumber)
00137   {
00138   case 1:
00139     read_pass2_1_dvdvdv();
00140     break;
00141   case 2:
00142     read_pass2_2_dvdvdv();
00143     break;
00144   case 3:
00145     read_pass2_3_dvdvdv();
00146     break;
00147   default:
00148     cerr << "illegal value for df1b2variable::pass = "
00149          << df1b2variable::passnumber << endl;
00150     exit(1);
00151   }
00152 }
00153 
00158 void read_pass2_1_dvdvdv(void)
00159 {
00160   // We are going backword for bptr and nbptr
00161   // and  forward for bptr2 and nbptr2
00162   // the current entry+2 in bptr is the size of the record i.e
00163   // points to the next record
00164   //char * bptr=f1b2gradlist->bptr;
00165   //char * bptr2=f1b2gradlist2->bptr;
00166   int nvar=df1b2variable::nvar;
00167   test_smartlist& list=f1b2gradlist->list;
00168   //f1b2gradlist->nlist-=sizeof(int);
00169   int num_bytes=f1b2gradlist->nlist.bptr->numbytes;
00170   list-=num_bytes;
00171   list.saveposition(); // save pointer to beginning of record;
00172   double xu,yu,wu;
00173   //ad_dstar xdot,ydot;
00174   //df1b2function2 * pf;
00175 
00176   // get info from tape1
00177 #if defined(SAFE_ARRAYS)
00178   checkidentiferstring("U8",f1b2gradlist->list);
00179 #endif
00180   char * bptr=f1b2gradlist->list.bptr;
00181   df1b2_header * px=(df1b2_header *) bptr;
00182   bptr+=sizeof(df1b2_header);
00183   df1b2_header * py=(df1b2_header *) bptr;
00184   bptr+=sizeof(df1b2_header);
00185   df1b2_header * pw=(df1b2_header *) bptr;
00186   bptr+=sizeof(df1b2_header);
00187   df1b2_header * pz=(df1b2_header *) bptr;
00188   bptr+=sizeof(df1b2_header);
00189   //pf=*(df1b2function2 **) bptr;
00190   //bptr+=sizeof(char*);
00191 
00192   double df1=*(double*) bptr;
00193   bptr+=sizeof(double);
00194 
00195   double df2=*(double*) bptr;
00196   bptr+=sizeof(double);
00197 
00198   double df3=*(double*) bptr;
00199   bptr+=sizeof(double);
00200 
00201   double d2f11=*(double*) bptr;
00202   bptr+=sizeof(double);
00203 
00204   double d2f12=*(double*) bptr;
00205   bptr+=sizeof(double);
00206 
00207   double d2f13=*(double*) bptr;
00208   bptr+=sizeof(double);
00209 
00210   double d2f22=*(double*) bptr;
00211   bptr+=sizeof(double);
00212 
00213   double d2f23=*(double*) bptr;
00214   bptr+=sizeof(double);
00215 
00216   double d2f33=*(double*) bptr;
00217   bptr+=sizeof(double);
00218 
00219 #if defined(PRINT_DERS)
00220   double d3f111=*(double*) bptr;
00221 #endif
00222   bptr+=sizeof(double);
00223 
00224 #if defined(PRINT_DERS)
00225   double d3f112=*(double*) bptr;
00226 #endif
00227   bptr+=sizeof(double);
00228 
00229 #if defined(PRINT_DERS)
00230   double d3f113=*(double*) bptr;
00231 #endif
00232   bptr+=sizeof(double);
00233 
00234 #if defined(PRINT_DERS)
00235   double d3f122=*(double*) bptr;
00236 #endif
00237   bptr+=sizeof(double);
00238 
00239 #if defined(PRINT_DERS)
00240   double d3f123=*(double*) bptr;
00241 #endif
00242   bptr+=sizeof(double);
00243 
00244 #if defined(PRINT_DERS)
00245   double d3f133=*(double*) bptr;
00246 #endif
00247   bptr+=sizeof(double);
00248 
00249 #if defined(PRINT_DERS)
00250   double d3f222=*(double*) bptr;
00251 #endif
00252   bptr+=sizeof(double);
00253 
00254 #if defined(PRINT_DERS)
00255   double d3f223=*(double*) bptr;
00256 #endif
00257   bptr+=sizeof(double);
00258 
00259 #if defined(PRINT_DERS)
00260   double d3f233=*(double*) bptr;
00261 #endif
00262   bptr+=sizeof(double);
00263 
00264 #if defined(PRINT_DERS)
00265   double d3f333=*(double*) bptr;
00266 #endif
00267   bptr+=sizeof(double);
00268 
00269   memcpy(&xu,bptr,sizeof(double));
00270   bptr+=sizeof(double);
00271   memcpy(&yu,bptr,sizeof(double));
00272   bptr+=sizeof(double);
00273   memcpy(&wu,bptr,sizeof(double));
00274   bptr+=sizeof(double);
00275   double * xdot=(double*)bptr;
00276   bptr+=nvar*sizeof(double);
00277   double * ydot=(double*)bptr;
00278   bptr+=nvar*sizeof(double);
00279   double * wdot=(double*)bptr;
00280 
00281   list.restoreposition(); // save pointer to beginning of record;
00282 
00283   // ****************************************************************
00284   // turn this off if no third derivatives are calculated
00285   // if (!no_third_derivatives)
00286   // {
00287   // save for second reverse pass
00288   // save identifier 1
00289      test_smartlist & list2 = f1b2gradlist->list2;
00290 
00291 
00292    size_t total_bytes=2*nvar*sizeof(double);
00293 // string identifier debug stuff
00294 #if defined(SAFE_ALL)
00295   char ids[]="FX";
00296   int slen=strlen(ids);
00297   total_bytes+=slen;
00298 #endif
00299 
00300 #ifndef OPT_LIB
00301   assert(total_bytes <= INT_MAX);
00302 #endif
00303 
00304   list2.check_buffer_size((int)total_bytes);
00305   void * tmpptr=list2.bptr;
00306 #if defined(SAFE_ALL)
00307   memcpy(list2,ids,slen);
00308 #endif
00309 
00310      fixed_smartlist2 & nlist2 = f1b2gradlist->nlist2;
00311   const int sizeofdouble = sizeof(double);
00312      memcpy(list2,pz->get_u_bar(),nvar*sizeofdouble);
00313      memcpy(list2,pz->get_u_dot_bar(),nvar*sizeofdouble);
00314      *nlist2.bptr=adptr_diff(list2.bptr,tmpptr);
00315      ++nlist2;
00316   // }
00317   //
00318   // ****************************************************************
00319 #if defined(PRINT_DERS)
00320  print_derivatives(funname,(f),(df1),
00321   (df2),(d2f11),(d2f12),(d2f22),
00322   (d3f111),(d3f112),(d3f122),
00323   (d3f222),1);
00324  print_derivatives(pz,"z");
00325  print_derivatives(px,"x");
00326  print_derivatives(py,"y");
00327 #endif
00328 #if defined(__DERCHECK__)
00329   if (derchecker)
00330   if (derchecker->node_number)
00331   {
00332     if (derchecker->counter == derchecker->node_number)
00333     {
00334       switch (derchecker->pass_number) // increment the variable of interest
00335       {
00336       case 2:
00337         switch(derchecker->vartype)
00338         {
00339         case 1:
00340           if (!derchecker->dotflag)
00341             px->u_bar[derchecker->index-1]+=derchecker->delta;
00342           else
00343             px->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00344           break;
00345         case 2:
00346           if (!derchecker->dotflag)
00347             py->u_bar[derchecker->index-1]+=derchecker->delta;
00348           else
00349             py->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00350           break;
00351         case 3:
00352           if (!derchecker->dotflag)
00353             pz->u_bar[derchecker->index-1]+=derchecker->delta;
00354           else
00355             pz->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00356           break;
00357         default:
00358           cerr << "Invalid index value for dercheck_index was "
00359                << derchecker->index << endl;
00360           break;
00361         }
00362         break;
00363       case 3:
00364         switch(derchecker->vartype)
00365         {
00366         case 1:
00367           if (!derchecker->dotflag)
00368             px->u_bar[derchecker->index-1]-=derchecker->delta;
00369           else
00370             px->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00371           break;
00372         case 2:
00373           if (!derchecker->dotflag)
00374             py->u_bar[derchecker->index-1]-=derchecker->delta;
00375           else
00376             py->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00377           break;
00378         case 3:
00379           if (!derchecker->dotflag)
00380             pz->u_bar[derchecker->index-1]-=derchecker->delta;
00381           else
00382             pz->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00383           break;
00384         default:
00385           cerr << "Invalid index value for dercheck_index was "
00386                << derchecker->index << endl;
00387           break;
00388         }
00389         break;
00390       }
00391     }
00392   }
00393 #endif
00394 
00395   // Do first reverse pass calculations
00396   int i;
00397   for (i=0;i<nvar;i++)
00398   {
00399     px->u_bar[i]+=(df1)*pz->u_bar[i];
00400   }
00401   for (i=0;i<nvar;i++)
00402   {
00403     py->u_bar[i]+=(df2)*pz->u_bar[i];
00404   }
00405   for (i=0;i<nvar;i++)
00406   {
00407     pw->u_bar[i]+=(df3)*pz->u_bar[i];
00408   }
00409   for (i=0;i<nvar;i++)
00410   {
00411     px->u_bar[i]+=(d2f11)*xdot[i]*pz->u_dot_bar[i];
00412     px->u_bar[i]+=(d2f12)*ydot[i]*pz->u_dot_bar[i];
00413     px->u_bar[i]+=(d2f13)*wdot[i]*pz->u_dot_bar[i];
00414   }
00415   for (i=0;i<nvar;i++)
00416   {
00417     py->u_bar[i]+=(d2f12)*xdot[i]*pz->u_dot_bar[i];
00418     py->u_bar[i]+=(d2f22)*ydot[i]*pz->u_dot_bar[i];
00419     py->u_bar[i]+=(d2f13)*wdot[i]*pz->u_dot_bar[i];
00420   }
00421   for (i=0;i<nvar;i++)
00422   {
00423     pw->u_bar[i]+=(d2f13)*xdot[i]*pz->u_dot_bar[i];
00424     pw->u_bar[i]+=(d2f23)*ydot[i]*pz->u_dot_bar[i];
00425     pw->u_bar[i]+=(d2f33)*wdot[i]*pz->u_dot_bar[i];
00426   }
00427   for (i=0;i<nvar;i++)
00428   {
00429     px->u_dot_bar[i]+=(df1)*pz->u_dot_bar[i];
00430   }
00431   for (i=0;i<nvar;i++)
00432   {
00433     py->u_dot_bar[i]+=(df2)*pz->u_dot_bar[i];
00434   }
00435   for (i=0;i<nvar;i++)
00436   {
00437     pw->u_dot_bar[i]+=(df3)*pz->u_dot_bar[i];
00438   }
00439 
00440   for (i=0;i<nvar;i++)
00441   {
00442     pz->u_bar[i]=0;
00443   }
00444   for (i=0;i<nvar;i++)
00445   {
00446     pz->u_dot_bar[i]=0;
00447   }
00448 }
00449 
00454 void read_pass2_2_dvdvdv(void)
00455 {
00456   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00457   // We are going forward for bptr and backword for bptr2
00458   //
00459   // list 1
00460   //
00461   int nvar=df1b2variable::nvar;
00462   test_smartlist & list=f1b2gradlist->list;
00463 
00464   size_t total_bytes=3*sizeof(df1b2_header)+sizeof(char*)
00465     +2*(nvar+1)*sizeof(double);
00466 // string identifier debug stuff
00467 #if defined(SAFE_ALL)
00468   char ids[]="U8";
00469   int slen=strlen(ids);
00470   total_bytes+=slen;
00471 #endif
00472 
00473 #ifndef OPT_LIB
00474   assert(total_bytes <= INT_MAX);
00475 #endif
00476 
00477   list.check_buffer_size((int)total_bytes);
00478 // end of string identifier debug stuff
00479 
00480   list.saveposition(); // save pointer to beginning of record;
00481   fixed_smartlist & nlist=f1b2gradlist->nlist;
00482    // nlist-=sizeof(int);
00483   // get record size
00484   int num_bytes=nlist.bptr->numbytes;
00485     // nlist+=nlist_record_size;
00486   //
00487   // list 2
00488   //
00489   test_smartlist & list2=f1b2gradlist->list2;
00490   fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00491   // get record size
00492   int num_bytes2=*nlist2.bptr;
00493   --nlist2;
00494   // backup the size of the record
00495   list2-=num_bytes2;
00496   list2.saveposition(); // save pointer to beginning of record;
00497   // save the pointer to the beginning of the record
00498   // bptr and bptr2 now both point to the beginning of their records
00499 
00500   double xu,yu,wu;
00501   //df1b2_header x,z;
00502   //df1b2function2 * pf;
00503 
00504   // get info from tape1
00505   // get info from tape1
00506 #if defined(SAFE_ARRAYS)
00507   checkidentiferstring("U8",list);
00508   checkidentiferstring("FX",list2);
00509 #endif
00510   /*
00511   df1b2_header * px=(df1b2_header *) list.bptr;
00512   list.bptr+=sizeof(df1b2_header);
00513   df1b2_header * py=(df1b2_header *) list.bptr;
00514   list.bptr+=sizeof(df1b2_header);
00515   df1b2_header * pz=(df1b2_header *) list.bptr;
00516   list.bptr+=sizeof(df1b2_header);
00517   pf=*(df1b2function2 **) list.bptr;
00518   list.bptr+=sizeof(char*);
00519   memcpy(&xu,list.bptr,sizeof(double));
00520   list.bptr+=sizeof(double);
00521   memcpy(&yu,list.bptr,sizeof(double));
00522   list.bptr+=sizeof(double);
00523   xdot=(double*)list.bptr;
00524   list.bptr+=nvar*sizeof(double);
00525   ydot=(double*)list.bptr;
00526   */
00527   //char * bptr=f1b2gradlist->list.bptr;
00528   df1b2_header * px=(df1b2_header *) list.bptr;
00529   list.bptr+=sizeof(df1b2_header);
00530   df1b2_header * py=(df1b2_header *) list.bptr;
00531   list.bptr+=sizeof(df1b2_header);
00532   df1b2_header * pw=(df1b2_header *) list.bptr;
00533   list.bptr+=sizeof(df1b2_header);
00534   df1b2_header * pz=(df1b2_header *) list.bptr;
00535   list.bptr+=sizeof(df1b2_header);
00536 
00537   double df1=*(double*) list.bptr;
00538   list.bptr+=sizeof(double);
00539 
00540   double df2=*(double*) list.bptr;
00541   list.bptr+=sizeof(double);
00542 
00543   double df3=*(double*) list.bptr;
00544   list.bptr+=sizeof(double);
00545 
00546   double d2f11=*(double*) list.bptr;
00547   list.bptr+=sizeof(double);
00548 
00549   double d2f12=*(double*) list.bptr;
00550   list.bptr+=sizeof(double);
00551 
00552   double d2f13=*(double*) list.bptr;
00553   list.bptr+=sizeof(double);
00554 
00555   double d2f22=*(double*) list.bptr;
00556   list.bptr+=sizeof(double);
00557 
00558   double d2f23=*(double*) list.bptr;
00559   list.bptr+=sizeof(double);
00560 
00561   double d2f33=*(double*) list.bptr;
00562   list.bptr+=sizeof(double);
00563 
00564   double d3f111=*(double*) list.bptr;
00565   list.bptr+=sizeof(double);
00566 
00567   double d3f112=*(double*) list.bptr;
00568   list.bptr+=sizeof(double);
00569 
00570   double d3f113=*(double*) list.bptr;
00571   list.bptr+=sizeof(double);
00572 
00573   double d3f122=*(double*) list.bptr;
00574   list.bptr+=sizeof(double);
00575 
00576   double d3f123=*(double*) list.bptr;
00577   list.bptr+=sizeof(double);
00578 
00579   double d3f133=*(double*) list.bptr;
00580   list.bptr+=sizeof(double);
00581 
00582   double d3f222=*(double*) list.bptr;
00583   list.bptr+=sizeof(double);
00584 
00585   double d3f223=*(double*) list.bptr;
00586   list.bptr+=sizeof(double);
00587 
00588   double d3f233=*(double*) list.bptr;
00589   list.bptr+=sizeof(double);
00590 
00591   double d3f333=*(double*) list.bptr;
00592   list.bptr+=sizeof(double);
00593 
00594   memcpy(&xu,list.bptr,sizeof(double));
00595   list.bptr+=sizeof(double);
00596   memcpy(&yu,list.bptr,sizeof(double));
00597   list.bptr+=sizeof(double);
00598   memcpy(&wu,list.bptr,sizeof(double));
00599   list.bptr+=sizeof(double);
00600   double * xdot=(double*)list.bptr;
00601   list.bptr+=nvar*sizeof(double);
00602   double * ydot=(double*)list.bptr;
00603   list.bptr+=nvar*sizeof(double);
00604   double * wdot=(double*)list.bptr;
00605 
00606   list.restoreposition(num_bytes); // save pointer to beginning of record;
00607 
00608   double * zbar;
00609   double * zdotbar;
00610 
00611 
00612   zbar=(double*)list2.bptr;
00613   zdotbar=(double*)(list2.bptr+nvar*sizeof(double));
00614   list2.restoreposition(); // save pointer to beginning of record;
00615 
00616   double * x_tilde=px->get_u_tilde();
00617   double * x_dot_tilde=px->get_u_dot_tilde();
00618   double * x_bar_tilde=px->get_u_bar_tilde();
00619   double * x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00620   double * y_tilde=py->get_u_tilde();
00621   double * y_dot_tilde=py->get_u_dot_tilde();
00622   double * y_bar_tilde=py->get_u_bar_tilde();
00623   double * y_dot_bar_tilde=py->get_u_dot_bar_tilde();
00624   double * w_tilde=pw->get_u_tilde();
00625   double * w_dot_tilde=pw->get_u_dot_tilde();
00626   double * w_bar_tilde=pw->get_u_bar_tilde();
00627   double * w_dot_bar_tilde=pw->get_u_dot_bar_tilde();
00628   double * z_bar_tilde=pz->get_u_bar_tilde();
00629   double * z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00630   // Do second "reverse-reverse" pass calculations
00631 #if defined(PRINT_DERS)
00632  print_derivatives(funname,(f),(df1),
00633   (df2),(d2f11),(d2f12),(d2f22),
00634   (d3f111),(d3f112),(d3f122),
00635   (d3f222),1);
00636  print_derivatives(pz,"z");
00637  print_derivatives(px,"x");
00638  print_derivatives(pw,"x");
00639  print_derivatives(py,"y");
00640 #endif
00641 
00642   int i;
00643 
00644   for (i=0;i<nvar;i++)
00645   {
00646     z_bar_tilde[i]=0;
00647     z_dot_bar_tilde[i]=0;
00648   }
00649 
00650   for (i=0;i<nvar;i++)
00651   {
00652     *x_tilde+=(d2f11)*zbar[i]*x_bar_tilde[i];
00653     *y_tilde+=(d2f12)*zbar[i]*x_bar_tilde[i];
00654     *w_tilde+=(d2f13)*zbar[i]*x_bar_tilde[i];
00655 
00656     *x_tilde+=(d2f12)*zbar[i]*y_bar_tilde[i];
00657     *y_tilde+=(d2f22)*zbar[i]*y_bar_tilde[i];
00658     *w_tilde+=(d2f23)*zbar[i]*y_bar_tilde[i];
00659 
00660     *x_tilde+=(d2f13)*zbar[i]*w_bar_tilde[i];
00661     *y_tilde+=(d2f23)*zbar[i]*w_bar_tilde[i];
00662     *w_tilde+=(d2f33)*zbar[i]*w_bar_tilde[i];
00663   }
00664 
00665   for (i=0;i<nvar;i++)
00666   {
00667     *x_tilde+=(d2f11)*zdotbar[i]*x_dot_bar_tilde[i];
00668     *y_tilde+=(d2f12)*zdotbar[i]*x_dot_bar_tilde[i];
00669     *w_tilde+=(d2f13)*zdotbar[i]*x_dot_bar_tilde[i];
00670 
00671     *x_tilde+=(d2f12)*zdotbar[i]*y_dot_bar_tilde[i];
00672     *y_tilde+=(d2f22)*zdotbar[i]*y_dot_bar_tilde[i];
00673     *w_tilde+=(d2f23)*zdotbar[i]*y_dot_bar_tilde[i];
00674 
00675     *x_tilde+=(d2f13)*zdotbar[i]*w_dot_bar_tilde[i];
00676     *y_tilde+=(d2f23)*zdotbar[i]*w_dot_bar_tilde[i];
00677     *w_tilde+=(d2f33)*zdotbar[i]*w_dot_bar_tilde[i];
00678   }
00679 
00680   for (i=0;i<nvar;i++)
00681   {
00682     *x_tilde+=(d3f111)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00683     *y_tilde+=(d3f112)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00684     *w_tilde+=(d3f113)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00685 
00686     *x_tilde+=(d3f112)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00687     *y_tilde+=(d3f122)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00688     *w_tilde+=(d3f123)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00689 
00690     *x_tilde+=(d3f113)*xdot[i]*zdotbar[i]*w_bar_tilde[i];
00691     *y_tilde+=(d3f123)*xdot[i]*zdotbar[i]*w_bar_tilde[i];
00692     *w_tilde+=(d3f133)*xdot[i]*zdotbar[i]*w_bar_tilde[i];
00693 
00694     *x_tilde+=(d3f112)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00695     *y_tilde+=(d3f122)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00696     *w_tilde+=(d3f123)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00697 
00698     *x_tilde+=(d3f122)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00699     *y_tilde+=(d3f222)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00700     *w_tilde+=(d3f223)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00701 
00702     *x_tilde+=(d3f123)*ydot[i]*zdotbar[i]*w_bar_tilde[i];
00703     *y_tilde+=(d3f223)*ydot[i]*zdotbar[i]*w_bar_tilde[i];
00704     *w_tilde+=(d3f233)*ydot[i]*zdotbar[i]*w_bar_tilde[i];
00705 
00706     *x_tilde+=(d3f113)*wdot[i]*zdotbar[i]*x_bar_tilde[i];
00707     *y_tilde+=(d3f123)*wdot[i]*zdotbar[i]*x_bar_tilde[i];
00708     *w_tilde+=(d3f133)*wdot[i]*zdotbar[i]*x_bar_tilde[i];
00709 
00710     *x_tilde+=(d3f123)*wdot[i]*zdotbar[i]*y_bar_tilde[i];
00711     *y_tilde+=(d3f223)*wdot[i]*zdotbar[i]*y_bar_tilde[i];
00712     *w_tilde+=(d3f233)*wdot[i]*zdotbar[i]*y_bar_tilde[i];
00713 
00714     *x_tilde+=(d3f133)*wdot[i]*zdotbar[i]*w_bar_tilde[i];
00715     *y_tilde+=(d3f233)*wdot[i]*zdotbar[i]*w_bar_tilde[i];
00716     *w_tilde+=(d3f333)*wdot[i]*zdotbar[i]*w_bar_tilde[i];
00717   }
00718 
00719   for (i=0;i<nvar;i++)
00720   {
00721     z_bar_tilde[i]+=(df1)*x_bar_tilde[i];
00722     z_bar_tilde[i]+=(df2)*y_bar_tilde[i];
00723     z_bar_tilde[i]+=(df3)*w_bar_tilde[i];
00724   }
00725 
00726   for (i=0;i<nvar;i++)
00727   {
00728     z_dot_bar_tilde[i]+=(df1)*x_dot_bar_tilde[i];
00729     z_dot_bar_tilde[i]+=(df2)*y_dot_bar_tilde[i];
00730     z_dot_bar_tilde[i]+=(df3)*w_dot_bar_tilde[i];
00731   }
00732 
00733   for (i=0;i<nvar;i++)
00734   {
00735     x_dot_tilde[i]+=(d2f11)*zdotbar[i]*x_bar_tilde[i];
00736     y_dot_tilde[i]+=(d2f12)*zdotbar[i]*x_bar_tilde[i];
00737     w_dot_tilde[i]+=(d2f13)*zdotbar[i]*x_bar_tilde[i];
00738 
00739     x_dot_tilde[i]+=(d2f12)*zdotbar[i]*y_bar_tilde[i];
00740     y_dot_tilde[i]+=(d2f22)*zdotbar[i]*y_bar_tilde[i];
00741     w_dot_tilde[i]+=(d2f23)*zdotbar[i]*y_bar_tilde[i];
00742 
00743     x_dot_tilde[i]+=(d2f13)*zdotbar[i]*w_bar_tilde[i];
00744     y_dot_tilde[i]+=(d2f23)*zdotbar[i]*w_bar_tilde[i];
00745     w_dot_tilde[i]+=(d2f33)*zdotbar[i]*w_bar_tilde[i];
00746   }
00747 
00748   for (i=0;i<nvar;i++)
00749   {
00750     z_dot_bar_tilde[i]+=(d2f11)*xdot[i]*x_bar_tilde[i];
00751     z_dot_bar_tilde[i]+=(d2f12)*xdot[i]*y_bar_tilde[i];
00752     z_dot_bar_tilde[i]+=(d2f13)*xdot[i]*w_bar_tilde[i];
00753 
00754     z_dot_bar_tilde[i]+=(d2f12)*ydot[i]*x_bar_tilde[i];
00755     z_dot_bar_tilde[i]+=(d2f22)*ydot[i]*y_bar_tilde[i];
00756     z_dot_bar_tilde[i]+=(d2f23)*ydot[i]*w_bar_tilde[i];
00757 
00758     z_dot_bar_tilde[i]+=(d2f13)*wdot[i]*x_bar_tilde[i];
00759     z_dot_bar_tilde[i]+=(d2f23)*wdot[i]*y_bar_tilde[i];
00760     z_dot_bar_tilde[i]+=(d2f33)*wdot[i]*w_bar_tilde[i];
00761   }
00762 
00763 
00764 #if defined(__DERCHECK__)
00765   if (derchecker->node_number)
00766   {
00767     if (derchecker->counter == derchecker->node_number)
00768     {
00769       if (derchecker->pass_number==1) // increment the variable of interest
00770       {
00771         switch(derchecker->vartype)
00772         {
00773         case 1:
00774           if (!derchecker->dotflag)
00775             derchecker->der_value=
00776               px->u_bar_tilde[derchecker->index-1];
00777           else
00778             derchecker->der_value=
00779               px->u_dot_bar_tilde[derchecker->index-1];
00780           break;
00781         case 2:
00782           if (!derchecker->dotflag)
00783             derchecker->der_value=
00784               py->u_bar_tilde[derchecker->index-1];
00785           else
00786             derchecker->der_value=
00787               py->u_dot_bar_tilde[derchecker->index-1];
00788           break;
00789         case 3:
00790           if (!derchecker->dotflag)
00791             derchecker->der_value=
00792               pz->u_bar_tilde[derchecker->index-1];
00793           else
00794             derchecker->der_value=
00795               pz->u_dot_bar_tilde[derchecker->index-1];
00796           break;
00797         default:
00798           cerr << "Invalid index value for dercheck_index was "
00799                << derchecker->index << endl;
00800         }
00801       }
00802     }
00803   }
00804 #endif
00805 #if defined(PRINT_DERS)
00806  print_derivatives(pz,"z");
00807  print_derivatives(px,"x");
00808  print_derivatives(py,"y");
00809 #endif
00810 }
00811 
00816 void read_pass2_3_dvdvdv(void)
00817 {
00818   // We are going backword for bptr and forward for bptr2
00819   // the current entry+2 in bptr is the size of the record i.e
00820   // points to the next record
00821   int nvar=df1b2variable::nvar;
00822   fixed_smartlist & nlist=f1b2gradlist->nlist;
00823   test_smartlist& list=f1b2gradlist->list;
00824    // nlist-=sizeof(int);
00825   // get record size
00826   int num_bytes=nlist.bptr->numbytes;
00827   // backup the size of the record
00828   list-=num_bytes;
00829   list.saveposition(); // save pointer to beginning of record;
00830   // save the pointer to the beginning of the record
00831   double xu;
00832   double yu;
00833   double wu;
00834   //df1b2_header x,z;
00835   //df1b2function2 * pf;
00836 
00837   // get info from tape1
00838   // get info from tape1
00839 #if defined(SAFE_ARRAYS)
00840   checkidentiferstring("U8",list);
00841 #endif
00842 
00843   df1b2_header * px=(df1b2_header *) list.bptr;
00844   list.bptr+=sizeof(df1b2_header);
00845   df1b2_header * py=(df1b2_header *) list.bptr;
00846   list.bptr+=sizeof(df1b2_header);
00847   df1b2_header * pw=(df1b2_header *) list.bptr;
00848   list.bptr+=sizeof(df1b2_header);
00849   df1b2_header * pz=(df1b2_header *) list.bptr;
00850   list.bptr+=sizeof(df1b2_header);
00851 
00852   double df1=*(double*) list.bptr;
00853   list.bptr+=sizeof(double);
00854 
00855   double df2=*(double*) list.bptr;
00856   list.bptr+=sizeof(double);
00857 
00858   double df3=*(double*) list.bptr;
00859   list.bptr+=sizeof(double);
00860 
00861   double d2f11=*(double*) list.bptr;
00862   list.bptr+=sizeof(double);
00863 
00864   double d2f12=*(double*) list.bptr;
00865   list.bptr+=sizeof(double);
00866 
00867   double d2f13=*(double*) list.bptr;
00868   list.bptr+=sizeof(double);
00869 
00870   double d2f22=*(double*) list.bptr;
00871   list.bptr+=sizeof(double);
00872 
00873   double d2f23=*(double*) list.bptr;
00874   list.bptr+=sizeof(double);
00875 
00876   double d2f33=*(double*) list.bptr;
00877   list.bptr+=sizeof(double);
00878 
00879 #if defined(PRINT_DERS)
00880   double d3f111=*(double*) list.bptr;
00881 #endif
00882   list.bptr+=sizeof(double);
00883 
00884 #if defined(PRINT_DERS)
00885   double d3f112=*(double*) list.bptr;
00886 #endif
00887   list.bptr+=sizeof(double);
00888 
00889 #if defined(PRINT_DERS)
00890   double d3f113=*(double*) list.bptr;
00891 #endif
00892   list.bptr+=sizeof(double);
00893 
00894 #if defined(PRINT_DERS)
00895   double d3f122=*(double*) list.bptr;
00896 #endif
00897   list.bptr+=sizeof(double);
00898 
00899 #if defined(PRINT_DERS)
00900   double d3f123=*(double*) list.bptr;
00901 #endif
00902   list.bptr+=sizeof(double);
00903 
00904 #if defined(PRINT_DERS)
00905   double d3f133=*(double*) list.bptr;
00906 #endif
00907   list.bptr+=sizeof(double);
00908 
00909 #if defined(PRINT_DERS)
00910   double d3f222=*(double*) list.bptr;
00911 #endif
00912   list.bptr+=sizeof(double);
00913 
00914 #if defined(PRINT_DERS)
00915   double d3f223=*(double*) list.bptr;
00916 #endif
00917   list.bptr+=sizeof(double);
00918 
00919 #if defined(PRINT_DERS)
00920   double d3f233=*(double*) list.bptr;
00921 #endif
00922   list.bptr+=sizeof(double);
00923 
00924 #if defined(PRINT_DERS)
00925   double d3f333=*(double*) list.bptr;
00926 #endif
00927   list.bptr+=sizeof(double);
00928 
00929   memcpy(&xu,list.bptr,sizeof(double));
00930   list.bptr+=sizeof(double);
00931   memcpy(&yu,list.bptr,sizeof(double));
00932   list.bptr+=sizeof(double);
00933   memcpy(&wu,list.bptr,sizeof(double));
00934   list.bptr+=sizeof(double);
00935   double * xdot=(double*)list.bptr;
00936   list.bptr+=nvar*sizeof(double);
00937   double * ydot=(double*)list.bptr;
00938   list.bptr+=nvar*sizeof(double);
00939   double * wdot=(double*)list.bptr;
00940 
00941   list.restoreposition(); // save pointer to beginning of record;
00942   int i;
00943 #if defined(PRINT_DERS)
00944  print_derivatives(funname,(f),(df1),
00945   (df2),(d2f11),(d2f12),(d2f22),
00946   (d3f111),(d3f112),(d3f122),
00947   (d3f222),1);
00948  print_derivatives(pz,"z");
00949  print_derivatives(px,"x");
00950  print_derivatives(py,"y");
00951  print_derivatives(pw,"y");
00952 #endif
00953 
00954   *(px->u_tilde)+=(df1)* *(pz->u_tilde);
00955   *(py->u_tilde)+=(df2)* *(pz->u_tilde);
00956   *(pw->u_tilde)+=(df3)* *(pz->u_tilde);
00957 
00958   for (i=0;i<nvar;i++)
00959   {
00960     *(px->u_tilde)+=(d2f11)*xdot[i]*pz->u_dot_tilde[i];
00961     *(py->u_tilde)+=(d2f12)*xdot[i]*pz->u_dot_tilde[i];
00962     *(pw->u_tilde)+=(d2f13)*xdot[i]*pz->u_dot_tilde[i];
00963 
00964     *(px->u_tilde)+=(d2f12)*ydot[i]*pz->u_dot_tilde[i];
00965     *(py->u_tilde)+=(d2f22)*ydot[i]*pz->u_dot_tilde[i];
00966     *(pw->u_tilde)+=(d2f23)*ydot[i]*pz->u_dot_tilde[i];
00967 
00968     *(px->u_tilde)+=(d2f13)*wdot[i]*pz->u_dot_tilde[i];
00969     *(py->u_tilde)+=(d2f23)*wdot[i]*pz->u_dot_tilde[i];
00970     *(pw->u_tilde)+=(d2f33)*wdot[i]*pz->u_dot_tilde[i];
00971   }
00972   for (i=0;i<nvar;i++)
00973   {
00974     px->u_dot_tilde[i]+=(df1)*pz->u_dot_tilde[i];
00975     py->u_dot_tilde[i]+=(df2)*pz->u_dot_tilde[i];
00976     pw->u_dot_tilde[i]+=(df3)*pz->u_dot_tilde[i];
00977   }
00978   *(pz->u_tilde)=0;
00979   for (i=0;i<nvar;i++)
00980   {
00981     pz->u_dot_tilde[i]=0;
00982   }
00983 
00984 
00985 #if defined(PRINT_DERS)
00986  print_derivatives(pz,"z");
00987  print_derivatives(px,"x");
00988  print_derivatives(py,"y");
00989  print_derivatives(pw,"y");
00990 #endif
00991 }