ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
f1b2v10.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: f1b2v10.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier and Mollie Brooks
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 
00017 df1b2vector pow(const df1b2vector& _v,double x)
00018 {
00019   ADUNCONST(df1b2vector,v);
00020   //check_shape(x,M,"operator *");
00021   int mmin=v.indexmin();
00022   int mmax=v.indexmax();
00023   df1b2vector tmp;
00024   tmp.noallocate(mmin,mmax);
00025   for (int i=mmin;i<=mmax;i++) tmp(i)=pow(v(i),x);
00026   return tmp;
00027 }
00028 
00033 df1b2vector pow(const  df1b2vector& _x,const df1b2vector& _v)
00034 {
00035   ADUNCONST(df1b2vector,x);
00036   ADUNCONST(df1b2vector,v);
00037   //check_shape(x,M,"operator *");
00038   //check_shape(x,M,"operator *");
00039   int mmin=v.indexmin();
00040   int mmax=v.indexmax();
00041   df1b2vector tmp;
00042   tmp.noallocate(mmin,mmax);
00043   for (int i=mmin;i<=mmax;i++) tmp(i)=pow(x(i),v(i));
00044   return tmp;
00045 }
00046 
00051 df1b2vector pow(const df1b2vector& _v,const df1b2variable & _x)
00052 {
00053   ADUNCONST(df1b2vector,v);
00054   ADUNCONST(df1b2variable,x);
00055   int mmin=v.indexmin();
00056   int mmax=v.indexmax();
00057   df1b2vector tmp;
00058   tmp.noallocate(mmin,mmax);
00059   for (int i=mmin;i<=mmax;i++)
00060   {
00061     tmp(i)=pow(v(i),x);
00062   }
00063   return tmp;
00064 }
00065 
00070 df1b2vector pow(const df1b2variable & _v,const df1b2vector& _x)
00071 {
00072   ADUNCONST(df1b2variable,v);
00073   ADUNCONST(df1b2vector,x);
00074   int mmin=x.indexmin();
00075   int mmax=x.indexmax();
00076   df1b2vector tmp;
00077   tmp.noallocate(mmin,mmax);
00078   for (int i=mmin;i<=mmax;i++)
00079   {
00080     tmp(i)=pow(v,x(i));
00081   }
00082   return tmp;
00083 }
00084 
00089 df1b2vector pow(const double v,const df1b2vector& _x)
00090 {
00091   ADUNCONST(df1b2vector,x);
00092   int mmin=x.indexmin();
00093   int mmax=x.indexmax();
00094   df1b2vector tmp;
00095   tmp.noallocate(mmin,mmax);
00096   for (int i=mmin;i<=mmax;i++)
00097   {
00098     tmp(i)=pow(v,x(i));
00099   }
00100   return tmp;
00101 }
00107 df1b2vector pow(const dvector& x,  const df1b2variable& a)
00108 {
00109   RETURN_ARRAYS_INCREMENT();
00110   df1b2vector y(x.indexmin(), x.indexmax());
00111   for(int i=x.indexmin(); i<=x.indexmax(); i++)
00112   {
00113     y(i)=pow(x(i),a);
00114   }
00115 
00116   RETURN_ARRAYS_DECREMENT();
00117   return(y);
00118 }
00124 df1b2vector pow(const dvector& x,  const df1b2vector& a)
00125 {
00126   RETURN_ARRAYS_INCREMENT();
00127   df1b2vector y(x.indexmin(), x.indexmax());
00128 
00129   for(int i=x.indexmin(); i<=x.indexmax(); i++)
00130   {
00131     y(i)=pow(x(i),a(i));
00132   }
00133 
00134   RETURN_ARRAYS_DECREMENT();
00135   return(y);
00136 }
00142   dvar_vector pow(const dvar_vector& v1, const dvector& v2)
00143   {
00144     RETURN_ARRAYS_INCREMENT();
00145     dvar_vector tmp(v1.indexmin(),v1.indexmax());
00146     for (int i=v1.indexmin();i<=v1.indexmax();i++)
00147     {
00148       tmp.elem(i)=pow(v1.elem(i),v2.elem(i));
00149     }
00150     RETURN_ARRAYS_DECREMENT();
00151     return(tmp);
00152   }
00158   df1b2vector pow(df1b2vector const& _x,dvector const& v)
00159   {
00160    ADUNCONST(df1b2vector,x);
00161    int mmin=v.indexmin();
00162    int mmax=v.indexmax();
00163    df1b2vector tmp;
00164    tmp.noallocate(mmin,mmax);
00165    for (int i=mmin;i<=mmax;i++) tmp(i)=pow(x(i),v(i));
00166    return tmp;
00167   }
00168 
00173 df1b2variable operator * (const df1b2vector& _x,const df1b2vector& _y)
00174 {
00175   ADUNCONST(df1b2vector,x)
00176   ADUNCONST(df1b2vector,y)
00177   df1b2variable tmp=0.0;
00178   int mmin=x.indexmin();
00179   int mmax=x.indexmax();
00180   int i;
00181   for (i=mmin;i<=mmax;i++)
00182   {
00183     *tmp.get_u() += *x(i).get_u() *  *y(i).get_u();
00184   }
00185   double * zd=tmp.get_u_dot();
00186   for (int j=0;j<df1b2variable::nvar;j++)
00187   {
00188     *zd++ = 0;
00189   }
00190   for (i=mmin;i<=mmax;i++)
00191   {
00192     double * zd=tmp.get_u_dot();
00193     double * xd=x(i).get_u_dot();
00194     double * yd=y(i).get_u_dot();
00195     double xu= *x(i).get_u();
00196     double yu= *y(i).get_u();
00197 
00198     for (int j=0;j<df1b2variable::nvar;j++)
00199     {
00200       *zd++ += yu * *xd++ + xu * *yd++;
00201     }
00202   }
00203   // WRITE WHATEVER ON TAPE
00204   if (!df1b2_gradlist::no_derivatives)
00205     f1b2gradlist->write_pass1_prod(&x,&y,&tmp);
00206 
00207   return tmp;
00208 }
00209 
00214 void ad_read_pass2_prod_vector(void);
00215 
00216  int df1b2_gradlist::write_pass1_prod(const df1b2vector * _px,
00217    const df1b2vector * _py,df1b2variable * pz)
00218  {
00219    ADUNCONST(df1b2vector*,px)
00220    ADUNCONST(df1b2vector*,py)
00221    ncount++;
00222 #if defined(CHECK_COUNT)
00223   if (ncount >= ncount_check)
00224     cout << ncount << endl;
00225 #endif
00226    int nvar=df1b2variable::nvar;
00227 
00228    int mmin=px->indexmin();
00229    int mmax=px->indexmax();
00230 
00231    //int total_bytes=3*sizeof(df1b2_header)
00232     // +2*(nvar+1)*sizeof(double);
00233    int total_bytes= 2*sizeof(int) + 2*(mmax-mmin+1)*sizeof(df1b2_header)
00234     + sizeof(df1b2_header) + 2*(mmax-mmin+1)*sizeof(double)
00235     + 2*(mmax-mmin+1)*nvar*sizeof(double);
00236 // string identifier debug stuff
00237 #if defined(SAFE_ALL)
00238   char ids[]="DL";
00239   int slen=strlen(ids);
00240   total_bytes+=slen;
00241 #endif
00242   list.check_buffer_size(total_bytes);
00243   void * tmpptr=list.bptr;
00244 #if defined(SAFE_ALL)
00245   memcpy(list,ids,slen);
00246 #endif
00247 // end of string identifier debug stuff
00248 
00249    int i;
00250    memcpy(list,&mmin,sizeof(int));
00251    memcpy(list,&mmax,sizeof(int));
00252    for (i=mmin;i<=mmax;i++)
00253    {
00254      memcpy(list,(df1b2_header*)( &(*px)(i)  ),sizeof(df1b2_header));
00255      memcpy(list,(df1b2_header*)( &(*py)(i)  ),sizeof(df1b2_header));
00256    }
00257    memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00258 
00259 
00260    for (i=mmin;i<=mmax;i++)
00261    {
00262      memcpy(list,(*px)(i).get_u(),sizeof(double));
00263      memcpy(list,(*py)(i).get_u(),sizeof(double));
00264    }
00265    for (i=mmin;i<=mmax;i++)
00266    {
00267      memcpy(list,(*px)(i).get_u_dot(),nvar*sizeof(double));
00268      memcpy(list,(*py)(i).get_u_dot(),nvar*sizeof(double));
00269    }
00270    // ***** write  record size
00271    nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00272    //cout << "XX " << nlist.bptr->numbytes << " " << total_bytes << endl;
00273    nlist.bptr->pf=(ADrfptr)(&ad_read_pass2_prod_vector);
00274       ++nlist;
00275    return 0;
00276  }
00277 
00278 
00279 
00280 void read_pass2_1_prod_vector(void);
00281 void read_pass2_2_prod_vector(void);
00282 void read_pass2_3_prod_vector(void);
00283 
00288 void ad_read_pass2_prod_vector(void)
00289 {
00290   switch(df1b2variable::passnumber)
00291   {
00292   case 1:
00293     read_pass2_1_prod_vector();
00294     break;
00295   case 2:
00296     read_pass2_2_prod_vector();
00297     break;
00298   case 3:
00299     read_pass2_3_prod_vector();
00300     break;
00301   default:
00302     cerr << "illegal value for df1b2variable::pass = "
00303          << df1b2variable::passnumber << endl;
00304     exit(1);
00305   }
00306 }
00307 
00312 void read_pass2_1_prod_vector(void)
00313 {
00314   int nvar=df1b2variable::nvar;
00315   test_smartlist& list=f1b2gradlist->list;
00316   int num_bytes=f1b2gradlist->nlist.bptr->numbytes;
00317   list-=num_bytes;
00318   list.saveposition(); // save pointer to beginning of record;
00319 
00320 #if defined(SAFE_ALL)
00321   checkidentiferstring("DL",f1b2gradlist->list);
00322 #endif
00323   char * bptr=f1b2gradlist->list.bptr;
00324   int& mmin = *(int *) bptr;
00325   bptr+=sizeof(int);
00326   int& mmax = *(int *) bptr;
00327   bptr+=sizeof(int);
00328   df1b2_header_ptr_vector px(mmin,mmax);
00329   df1b2_header_ptr_vector py(mmin,mmax);
00330   double_ptr_vector xdot(mmin,mmax);
00331   double_ptr_vector ydot(mmin,mmax);
00332   dvector xu(mmin,mmax);
00333   dvector yu(mmin,mmax);
00334   int i;
00335   for (i=mmin;i<=mmax;i++)
00336   {
00337     // df1b2_header *
00338     px(i)=(df1b2_header *) bptr;
00339     bptr+=sizeof(df1b2_header);
00340     // df1b2_header *
00341     py(i)=(df1b2_header *) bptr;
00342     bptr+=sizeof(df1b2_header);
00343   }
00344   df1b2_header * pz=(df1b2_header *) bptr;
00345   bptr+=sizeof(df1b2_header);
00346   for (i=mmin;i<=mmax;i++)
00347   {
00348     memcpy(&(xu(i)),bptr,sizeof(double));
00349     bptr+=sizeof(double);
00350     memcpy(&(yu(i)),bptr,sizeof(double));
00351     bptr+=sizeof(double);
00352   }
00353   for (i=mmin;i<=mmax;i++)
00354   {
00355     // double *
00356     xdot(i)=(double*)bptr;
00357     bptr+=nvar*sizeof(double);
00358     // double *
00359     ydot(i)=(double*)bptr;
00360     bptr+=nvar*sizeof(double);
00361   }
00362 
00363   list.restoreposition(); // save pointer to beginning of record;
00364 
00365   // ****************************************************************
00366   // turn this off if no third derivatives are calculated
00367   // if (!no_third_derivatives)
00368   // {
00369   // save for second reverse pass
00370   // save identifier 1
00371      test_smartlist & list2 = f1b2gradlist->list2;
00372 
00373 
00374    int total_bytes=2*nvar*sizeof(double);
00375 // string identifier debug stuff
00376 #if defined(SAFE_ALL)
00377   char ids[]="QK";
00378   int slen=strlen(ids);
00379   total_bytes+=slen;
00380 #endif
00381   list2.check_buffer_size(total_bytes);
00382   void * tmpptr=list2.bptr;
00383 #if defined(SAFE_ALL)
00384   memcpy(list2,ids,slen);
00385 #endif
00386 
00387   fixed_smartlist2 & nlist2 = f1b2gradlist->nlist2;
00388   memcpy(list2,pz->get_u_bar(),nvar*sizeof(double));
00389   memcpy(list2,pz->get_u_dot_bar(),nvar*sizeof(double));
00390   *nlist2.bptr=adptr_diff(list2.bptr,tmpptr);
00391   ++nlist2;
00392 
00393   int j;
00394   // Do first reverse pass calculations
00395   for (i=mmin;i<=mmax;i++)
00396   {
00397     for (j=0;j<nvar;j++)
00398     {
00399       px(i)->u_bar[j]+=yu(i)*pz->u_bar[j];
00400     }
00401     for (j=0;j<nvar;j++)
00402     {
00403       py(i)->u_bar[j]+=xu(i)*pz->u_bar[j];
00404     }
00405   }
00406 
00407   for (i=mmin;i<=mmax;i++)
00408   {
00409     for (j=0;j<nvar;j++)
00410     {
00411       px(i)->u_bar[j]+=ydot(i)[j]*pz->u_dot_bar[j];
00412     }
00413 
00414     for (j=0;j<nvar;j++)
00415     {
00416       py(i)->u_bar[j]+=xdot(i)[j]*pz->u_dot_bar[j];
00417     }
00418     for (j=0;j<nvar;j++)
00419     {
00420       px(i)->u_dot_bar[j]+=yu(i)*pz->u_dot_bar[j];
00421     }
00422     for (j=0;j<nvar;j++)
00423     {
00424       py(i)->u_dot_bar[j]+=xu(i)*pz->u_dot_bar[j];
00425     }
00426   }
00427 
00428   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00429   for (j=0;j<nvar;j++)
00430   {
00431     pz->u_bar[j]=0;
00432   }
00433   for (j=0;j<nvar;j++)
00434   {
00435     pz->u_dot_bar[j]=0;
00436   }
00437 }
00438 
00443 void read_pass2_2_prod_vector(void)
00444 {
00445   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00446   // We are going forward for bptr and backword for bptr2
00447   //
00448   // list 1
00449   //
00450   int nvar=df1b2variable::nvar;
00451   test_smartlist & list=f1b2gradlist->list;
00452 
00453   //int total_bytes=3*sizeof(df1b2_header)+sizeof(char*)
00454   //  +2*(nvar+1)*sizeof(double);
00455   //int total_bytes=3*sizeof(df1b2_header)
00456   //  +2*(nvar+1)*sizeof(double);
00457 // string identifier debug stuff
00458 #if defined(SAFE_ALL)
00459   //char ids[]="DL";
00460   //int slen=strlen(ids);
00461   //total_bytes+=slen;
00462 #endif
00463 // end of string identifier debug stuff
00464 
00465   fixed_smartlist & nlist=f1b2gradlist->nlist;
00466   int total_bytes=nlist.bptr->numbytes;
00467   list.check_buffer_size(total_bytes);
00468   list.saveposition(); // save pointer to beginning of record;
00469   //cout << "YY " << total_bytes <<  endl;
00470   //
00471   // list 2
00472   //
00473   test_smartlist & list2=f1b2gradlist->list2;
00474   fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00475   // get record size
00476   int num_bytes2=*nlist2.bptr;
00477   --nlist2;
00478   // backup the size of the record
00479   list2-=num_bytes2;
00480   list2.saveposition(); // save pointer to beginning of record;
00481   // save the pointer to the beginning of the record
00482   // bptr and bptr2 now both point to the beginning of their records
00483 
00484   //df1b2_header x,z;
00485   //df1b2function2 * pf;
00486 
00487   // get info from tape1
00488   // get info from tape1
00489 #if defined(SAFE_ALL)
00490   checkidentiferstring("DL",list);
00491   checkidentiferstring("QK",list2);
00492 #endif
00493   int mmin = *(int *) list.bptr;
00494   list.bptr+=sizeof(int);
00495   int mmax = *(int *) list.bptr;
00496   list.bptr+=sizeof(int);
00497 
00498   df1b2_header_ptr_vector px(mmin,mmax);
00499   df1b2_header_ptr_vector py(mmin,mmax);
00500   double_ptr_vector xdot(mmin,mmax);
00501   double_ptr_vector ydot(mmin,mmax);
00502   dvector xu(mmin,mmax);
00503   dvector yu(mmin,mmax);
00504   int i;
00505   int j;
00506 
00507   for (i=mmin;i<=mmax;i++)
00508   {
00509     // df1b2_header * //
00510     px(i)=(df1b2_header *) list.bptr;
00511     list.bptr+=sizeof(df1b2_header);
00512     // df1b2_header *
00513     py(i)=(df1b2_header *) list.bptr;
00514     list.bptr+=sizeof(df1b2_header);
00515   }
00516 
00517   df1b2_header * pz=(df1b2_header *) list.bptr;
00518   list.bptr+=sizeof(df1b2_header);
00519 
00520   for (i=mmin;i<=mmax;i++)
00521   {
00522     memcpy(&(xu(i)),list.bptr,sizeof(double));
00523     list.bptr+=sizeof(double);
00524     memcpy(&(yu(i)),list.bptr,sizeof(double));
00525     list.bptr+=sizeof(double);
00526   }
00527   for (i=mmin;i<=mmax;i++)
00528   {
00529     xdot(i)=(double*)list.bptr;
00530     list.bptr+=nvar*sizeof(double);
00531     ydot(i)=(double*)list.bptr;
00532     list.bptr+=nvar*sizeof(double);
00533   }
00534   list.restoreposition(total_bytes); // save pointer to beginning of record;
00535 
00536   double * zbar;
00537   double * zdotbar;
00538 
00539 
00540   zbar=(double*)list2.bptr;
00541   zdotbar=(double*)(list2.bptr+nvar*sizeof(double));
00542   list2.restoreposition(); // save pointer to beginning of record;
00543 
00544   double * z_bar_tilde=pz->get_u_bar_tilde();
00545   double * z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00546 
00547   for (j=0;j<nvar;j++)
00548   {
00549     z_bar_tilde[j]=0;
00550     z_dot_bar_tilde[j]=0;
00551   }
00552 
00553   for (i=mmin;i<=mmax;i++)
00554   {
00555     double * x_tilde=px(i)->get_u_tilde();
00556     double * x_dot_tilde=px(i)->get_u_dot_tilde();
00557     double * x_bar_tilde=px(i)->get_u_bar_tilde();
00558     double * x_dot_bar_tilde=px(i)->get_u_dot_bar_tilde();
00559     double * y_tilde=py(i)->get_u_tilde();
00560     double * y_dot_tilde=py(i)->get_u_dot_tilde();
00561     double * y_bar_tilde=py(i)->get_u_bar_tilde();
00562     double * y_dot_bar_tilde=py(i)->get_u_dot_bar_tilde();
00563 
00564     int j;
00565 
00566     // start wjth x and add y
00567     for (j=0;j<nvar;j++)
00568     {
00569       z_bar_tilde[j]+=yu(i)*x_bar_tilde[j];
00570       *y_tilde+=zbar[j]*x_bar_tilde[j];
00571     }
00572 
00573     for (j=0;j<nvar;j++)
00574     {
00575       *y_tilde+=zdotbar[j]*x_dot_bar_tilde[j];
00576       z_dot_bar_tilde[j]+=yu(i)*x_dot_bar_tilde[j];
00577     }
00578 
00579     // start wjth y and add x
00580     for (j=0;j<nvar;j++)
00581     {
00582       *x_tilde+=zbar[j]*y_bar_tilde[j];
00583       z_bar_tilde[j]+=xu(i)*y_bar_tilde[j];
00584     }
00585 
00586     for (j=0;j<nvar;j++)
00587     {
00588       *x_tilde+=zdotbar[j]*y_dot_bar_tilde[j];
00589       z_dot_bar_tilde[j]+=xu(i)*y_dot_bar_tilde[j];
00590     }
00591 
00592     for (j=0;j<nvar;j++)
00593     {
00594       y_dot_tilde[j]+=zdotbar[j]*x_bar_tilde[j];
00595       z_dot_bar_tilde[j]+=ydot(i)[j]*x_bar_tilde[j];
00596     }
00597     for (j=0;j<nvar;j++)
00598     {
00599       x_dot_tilde[j]+=zdotbar[j]*y_bar_tilde[j];
00600       z_dot_bar_tilde[j]+=xdot(i)[j]*y_bar_tilde[j];
00601     }
00602   }
00603 }
00604 
00609 void read_pass2_3_prod_vector(void)
00610 {
00611   // We are going backword for bptr and forward for bptr2
00612   // the current entry+2 in bptr is the size of the record i.e
00613   // points to the next record
00614   int nvar=df1b2variable::nvar;
00615   fixed_smartlist & nlist=f1b2gradlist->nlist;
00616   test_smartlist& list=f1b2gradlist->list;
00617    // nlist-=sizeof(int);
00618   // get record size
00619   int num_bytes=nlist.bptr->numbytes;
00620   // backup the size of the record
00621   list-=num_bytes;
00622   list.saveposition(); // save pointer to beginning of record;
00623   // save the pointer to the beginning of the record
00624 
00625   //df1b2_header x,z;
00626   //df1b2function2 * pf;
00627 
00628   // get info from tape1
00629   // get info from tape1
00630 #if defined(SAFE_ALL)
00631   checkidentiferstring("DL",list);
00632 #endif
00633   int& mmin = *(int *) list.bptr;
00634   list.bptr+=sizeof(int);
00635   int& mmax = *(int *) list.bptr;
00636   list.bptr+=sizeof(int);
00637   df1b2_header_ptr_vector px(mmin,mmax);
00638   df1b2_header_ptr_vector py(mmin,mmax);
00639   double_ptr_vector xdot(mmin,mmax);
00640   double_ptr_vector ydot(mmin,mmax);
00641   dvector xu(mmin,mmax);
00642   dvector yu(mmin,mmax);
00643   int i;
00644   int j;
00645   for (i=mmin;i<=mmax;i++)
00646   {
00647     // df1b2_header *
00648     px(i)=(df1b2_header *) list.bptr;
00649     list.bptr+=sizeof(df1b2_header);
00650     // df1b2_header *
00651     py(i)=(df1b2_header *) list.bptr;
00652     list.bptr+=sizeof(df1b2_header);
00653   }
00654   df1b2_header * pz=(df1b2_header *) list.bptr;
00655   list.bptr+=sizeof(df1b2_header);
00656   for (i=mmin;i<=mmax;i++)
00657   {
00658     memcpy(&(xu(i)),list.bptr,sizeof(double));
00659     list.bptr+=sizeof(double);
00660     memcpy(&(yu(i)),list.bptr,sizeof(double));
00661     list.bptr+=sizeof(double);
00662   }
00663   for (i=mmin;i<=mmax;i++)
00664   {
00665     // double *
00666     xdot(i)=(double*)list.bptr;
00667     list.bptr+=nvar*sizeof(double);
00668     // double *
00669     ydot(i)=(double*)list.bptr;
00670     list.bptr+=nvar*sizeof(double);
00671   }
00672 
00673   list.restoreposition(); // save pointer to beginning of record;
00674 
00675   for (i=mmin;i<=mmax;i++)
00676   {
00677     *(px(i)->u_tilde)+=yu(i) * *(pz->u_tilde);
00678     *(py(i)->u_tilde)+=xu(i) * *(pz->u_tilde);
00679     for (j=0;j<nvar;j++)
00680     {
00681       *(py(i)->u_tilde)+=xdot(i)[j]*pz->u_dot_tilde[j];
00682       *(px(i)->u_tilde)+=ydot(i)[j]*pz->u_dot_tilde[j];
00683     }
00684     for (j=0;j<nvar;j++)
00685     {
00686       px(i)->u_dot_tilde[j]+=yu(i)*pz->u_dot_tilde[j];
00687       py(i)->u_dot_tilde[j]+=xu(i)*pz->u_dot_tilde[j];
00688     }
00689   }
00690   *(pz->u_tilde)=0;
00691   for (j=0;j<nvar;j++)
00692   {
00693     pz->u_dot_tilde[j]=0;
00694   }
00695 }