ADMB Documentation  11.5.3247
 All Classes Files Functions Variables Typedefs Friends Defines
df32fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 df1b2variable* df3_two_variable::ind_var[2];
00013 int df3_two_variable::num_ind_var = 0;
00014 
00018 df3_two_variable::df3_two_variable(const df3_two_variable& x)
00019 {
00020   v[0] = x.v[0];
00021   v[1] = x.v[1];
00022   v[2] = x.v[2];
00023   v[3] = x.v[3];
00024   v[4] = x.v[4];
00025   v[5] = x.v[5];
00026   v[6] = x.v[6];
00027   v[7] = x.v[7];
00028   v[8] = x.v[8];
00029   v[9] = x.v[9];
00030 }
00031 
00036  df3_two_vector::df3_two_vector(const df3_two_vector& m2)
00037  {
00038    index_min=m2.index_min;
00039    index_max=m2.index_max;
00040    shape=m2.shape;
00041    if (shape)
00042    {
00043      (shape->ncopies)++;
00044    }
00045    v = m2.v;
00046  }
00047 
00052  df3_two_vector::~df3_two_vector()
00053  {
00054    deallocate();
00055  }
00056 
00061  void df3_two_vector::deallocate(void)
00062  {
00063    if(shape)
00064    {
00065      if (shape->ncopies)
00066      {
00067        (shape->ncopies)--;
00068      }
00069      else
00070      {
00071        v = (df3_two_variable*) (shape->trueptr);
00072        delete [] v;
00073        v = NULL;
00074        delete shape;
00075        shape=0;
00076      }
00077    }
00078  }
00079 
00084  dvector value(const df3_two_vector& v)
00085  {
00086    int mmin=v.indexmin();
00087    int mmax=v.indexmax();
00088    dvector cv(mmin,mmax);
00089    for (int i=mmin;i<=mmax;i++)
00090    {
00091      cv(i)=value(v(i));
00092    }
00093    return cv;
00094  }
00095 
00100   void df3_two_vector::initialize(void)
00101   {
00102     int mmin=indexmin();
00103     int mmax=indexmax();
00104     for (int i=mmin;i<=mmax;i++)
00105     {
00106       (*this)(i)=0.0;
00107     }
00108   }
00109 
00114   df3_two_vector::df3_two_vector(void)
00115   {
00116     allocate();
00117   }
00118 
00123   df3_two_vector::df3_two_vector(int min,int max)
00124   {
00125     allocate(min,max);
00126   }
00127 
00132   void df3_two_vector::allocate(int min,int max)
00133   {
00134     index_min=min;
00135     index_max=max;
00136     v=new df3_two_variable[max-min+1];
00137     if (v==0)
00138     {
00139       cerr << "error allocating memory in df3_two_vector" << endl;
00140       ad_exit(1);
00141     }
00142     if ( (shape=new vector_shapex(min,max,v)) == NULL)
00143     {
00144       cerr << "Error trying to allocate memory for df3_two_vector"
00145            << endl;;
00146       ad_exit(1);
00147     }
00148     v-=min;
00149   }
00150 
00155   void df3_two_vector::allocate(void)
00156   {
00157     index_min=0;
00158     index_max=-1;
00159     v=0;
00160     shape=0;
00161   }
00162 
00167  dmatrix value(const df3_two_matrix& v)
00168  {
00169    int rmin=v.indexmin();
00170    int rmax=v.indexmax();
00171    dmatrix cm(rmin,rmax);
00172    for (int i=rmin;i<=rmax;i++)
00173    {
00174      int cmin=v(i).indexmin();
00175      int cmax=v(i).indexmax();
00176      cm(i).allocate(cmin,cmax);
00177      for (int j=cmin;j<=cmax;j++)
00178      {
00179        cm(i,j)=value(v(i,j));
00180      }
00181    }
00182    return cm;
00183  }
00184 
00189  df3_two_matrix::df3_two_matrix(const df3_two_matrix& m2)
00190  {
00191    index_min=m2.index_min;
00192    index_max=m2.index_max;
00193    shape=m2.shape;
00194    if (shape)
00195    {
00196      (shape->ncopies)++;
00197    }
00198    v = m2.v;
00199  }
00200 
00205  df3_two_matrix::~df3_two_matrix()
00206  {
00207    deallocate();
00208  }
00209 
00214  void df3_two_matrix::deallocate(void)
00215  {
00216    if (shape)
00217    {
00218      if (shape->ncopies)
00219      {
00220        (shape->ncopies)--;
00221      }
00222      else
00223      {
00224        v = (df3_two_vector*) (shape->get_pointer());
00225        delete [] v;
00226        v=0;
00227        delete shape;
00228        shape=0;
00229      }
00230    }
00231  }
00232 
00237   void df3_two_matrix::initialize(void)
00238   {
00239     int mmin=indexmin();
00240     int mmax=indexmax();
00241     for (int i=mmin;i<=mmax;i++)
00242     {
00243       (*this)(i).initialize();
00244     }
00245   }
00246 
00251   df3_two_matrix::df3_two_matrix(int rmin,int rmax,int cmin,int cmax)
00252   {
00253     index_min=rmin;
00254     index_max=rmax;
00255     v=new df3_two_vector[rmax-rmin+1];
00256     if (v==0)
00257     {
00258       cerr << "error allocating memory in df3_two_matrix" << endl;
00259       ad_exit(1);
00260     }
00261     if ( (shape=new mat_shapex(v)) == NULL)
00262     {
00263       cerr << "Error trying to allocate memory for df3_two_vector"
00264            << endl;;
00265     }
00266     v-=rmin;
00267 
00268     for (int i=rmin;i<=rmax;i++)
00269     {
00270       v[i].allocate(cmin,cmax);
00271     }
00272   }
00276 df3_two_variable& df3_two_variable::operator-=(const df3_two_variable& _v)
00277 {
00278   *get_u() -= *_v.get_u();
00279   *get_u_x() -= *_v.get_u_x();
00280   *get_u_y() -= *_v.get_u_y();
00281   *get_u_xx() -= *_v.get_u_xx();
00282   *get_u_xy() -= *_v.get_u_xy();
00283   *get_u_yy() -= *_v.get_u_yy();
00284   *get_u_xxx() -= *_v.get_u_xxx();
00285   *get_u_xxy() -= *_v.get_u_xxy();
00286   *get_u_xyy() -= *_v.get_u_xyy();
00287   *get_u_yyy() -= *_v.get_u_yyy();
00288   return *this;
00289 }
00290 
00295 df3_two_variable operator-(const df3_two_variable& v)
00296 {
00297   df3_two_variable z;
00298 
00299   *z.get_u() = - *v.get_u();
00300 
00301   *z.get_u_x() = -(*v.get_u_x());
00302   *z.get_u_y() = -(*v.get_u_y());
00303   *z.get_u_xx() = -(*v.get_u_xx());
00304   *z.get_u_xy() = -(*v.get_u_xy());
00305   *z.get_u_yy() = -(*v.get_u_yy());
00306   *z.get_u_xxx() = -(*v.get_u_xxx());
00307   *z.get_u_xxy() = -(*v.get_u_xxy());
00308   *z.get_u_xyy() = -(*v.get_u_xyy());
00309   *z.get_u_yyy() = -(*v.get_u_yyy());
00310 
00311   return z;
00312 }
00316 df3_two_variable& df3_two_variable::operator+=(const df3_two_variable& _v)
00317 {
00318   *get_u() += *_v.get_u();
00319   *get_u_x() += *_v.get_u_x();
00320   *get_u_y() += *_v.get_u_y();
00321   *get_u_xx() += *_v.get_u_xx();
00322   *get_u_xy() += *_v.get_u_xy();
00323   *get_u_yy() += *_v.get_u_yy();
00324   *get_u_xxx() += *_v.get_u_xxx();
00325   *get_u_xxy() += *_v.get_u_xxy();
00326   *get_u_xyy() += *_v.get_u_xyy();
00327   *get_u_yyy() += *_v.get_u_yyy();
00328 
00329   return *this;
00330 }
00334 df3_two_variable& df3_two_variable::operator+=(double _v)
00335 {
00336   *get_u() += _v;
00337 
00338   return *this;
00339 }
00343 df3_two_variable& df3_two_variable::operator-=(double _v)
00344 {
00345   *get_u() -= _v;
00346   return *this;
00347 }
00351 df3_two_variable& df3_two_variable::operator*=(const df3_two_variable& _v)
00352 {
00353   df3_two_variable x = *this * _v;
00354   *this = x;
00355   return *this;
00356 }
00360 df3_two_variable& df3_two_variable::operator*=(double _v)
00361 {
00362   *get_u() *= _v;
00363   *get_u_x() *= _v;
00364   *get_u_y() *= _v;
00365   *get_u_xx() *= _v;
00366   *get_u_xy() *= _v;
00367   *get_u_yy() *= _v;
00368   *get_u_xxx() *= _v;
00369   *get_u_xxy() *= _v;
00370   *get_u_xyy() *= _v;
00371   *get_u_yyy() *= _v;
00372 
00373   return *this;
00374 }
00375 
00380   df3_two_variable& df3_two_variable::operator /= (const df3_two_variable& y)
00381   {
00382     df3_two_variable x=*this / y;
00383     *this=x;
00384     return *this;
00385   }
00386 
00391 int operator <(const df3_two_variable & x, double n)
00392 {
00393    return value(x) < n;
00394 }
00395 
00400 int operator >(const df3_two_variable & x, double n)
00401 {
00402    return value(x) > n;
00403 }
00404 
00409 int operator >=(const df3_two_variable & x, double n)
00410 {
00411    return value(x) >= n;
00412 }
00413 
00418 int operator ==(const df3_two_variable & x, const df3_two_variable & n)
00419 {
00420    return value(x) == value(n);
00421 }
00422 
00427 int operator ==(const df3_two_variable & x, double n)
00428 {
00429    return value(x) == n;
00430 }
00431 
00436 int operator ==(double x, const df3_two_variable & n)
00437 {
00438    return x == value(n);
00439 }
00440 
00445 int operator <(const df3_two_variable & x, const df3_two_variable & n)
00446 {
00447    return value(x) < value(n);
00448 }
00449 
00454 int operator >(const df3_two_variable & x, const df3_two_variable & n)
00455 {
00456    return value(x) > value(n);
00457 }
00458 
00463 void set_derivatives( df3_two_variable& z,const df3_two_variable& x,double u,
00464   double zp,double zp2,double zp3)
00465 {
00466     //*z.get_u() = u;
00467 
00468     *z.get_u_x() = zp* *x.get_u_x();
00469 
00470     *z.get_u_y() = zp* *x.get_u_y();
00471 
00472     *z.get_u_xx() = zp2 * square(*x.get_u_x())
00473                    + zp * *x.get_u_xx();
00474 
00475     *z.get_u_xy() = zp2 * *x.get_u_x() * *x.get_u_y()
00476                    + zp * *x.get_u_xy();
00477 
00478     *z.get_u_yy() = zp2 * square(*x.get_u_y())
00479                    + zp * *x.get_u_yy();
00480 
00481     *z.get_u_xxx() = zp3 * cube(*x.get_u_x())
00482                    + 3.0 * zp2 * *x.get_u_x() * *x.get_u_xx()
00483                    + zp * *x.get_u_xxx();
00484 
00485     *z.get_u_xxy() = zp3 * square(*x.get_u_x())* *x.get_u_y()
00486                    + 2.0* zp2 * *x.get_u_x() * *x.get_u_xy()
00487                    + zp2 * *x.get_u_y() * *x.get_u_xx()
00488                    + zp * *x.get_u_xxy();
00489 
00490     *z.get_u_xyy() = zp3 * square(*x.get_u_y())* *x.get_u_x()
00491                    + 2.0* zp2 * *x.get_u_y() * *x.get_u_xy()
00492                    + zp2 * *x.get_u_x() * *x.get_u_yy()
00493                    + zp * *x.get_u_xyy();
00494 
00495 
00496     *z.get_u_yyy() = zp3 * cube(*x.get_u_y())
00497                    + 3.0 * zp2 * *x.get_u_y() * *x.get_u_yy()
00498                    + zp * *x.get_u_yyy();
00499 }
00500 
00505 void set_derivatives( df3_two_variable& z, const df3_two_variable& x,
00506   const df3_two_variable& y, double u,
00507   double f_u,double f_v,double f_uu,double f_uv,double f_vv,double f_uuu,
00508   double f_uuv,double f_uvv,double f_vvv)
00509 {
00510     *z.get_u() = u;
00511 
00512     *z.get_u_x() = f_u* *x.get_u_x()
00513                  + f_v* *y.get_u_x();
00514 
00515     *z.get_u_y() = f_u* *x.get_u_y()
00516                  + f_v* *y.get_u_y();
00517 
00518     *z.get_u_xx() = f_uu * square(*x.get_u_x())
00519                   + f_u  * *x.get_u_xx()
00520                   + f_vv * square(*y.get_u_x())
00521                   + f_v  * *y.get_u_xx()
00522             + 2.0 * f_uv * *x.get_u_x() * *y.get_u_x();
00523 
00524     *z.get_u_xy() = f_uu * *x.get_u_x() * *x.get_u_y()
00525                   + f_u  * *x.get_u_xy()
00526                   + f_vv *  *y.get_u_x() * *y.get_u_y()
00527                   + f_v  * *y.get_u_xy()
00528                   + f_uv * (*x.get_u_x() * *y.get_u_y()
00529                          +  *x.get_u_y() * *y.get_u_x());
00530 
00531     *z.get_u_yy() = f_uu * square(*x.get_u_y())
00532                   + f_u  * *x.get_u_yy()
00533                   + f_vv * square(*y.get_u_y())
00534                   + f_v  * *y.get_u_yy()
00535             + 2.0 * f_uv * *x.get_u_y() * *y.get_u_y();
00536 
00537 
00538     *z.get_u_xxx() = f_uuu * cube(*x.get_u_x())
00539                    + 3.0 * f_uu * *x.get_u_x() * *x.get_u_xx()
00540                    + f_u * *x.get_u_xxx()
00541                    + f_vvv * cube(*y.get_u_x())
00542                    + 3.0 * f_vv * *y.get_u_x() * *y.get_u_xx()
00543                    + f_v * *y.get_u_xxx()
00544                    + 3.0 * f_uuv * square(*x.get_u_x()) * *y.get_u_x()
00545                    + 3.0 * f_uvv * *x.get_u_x()* square(*y.get_u_x())
00546                    + 3.0* f_uv * *x.get_u_xx() * *y.get_u_x()
00547                    + 3.0* f_uv * *x.get_u_x() * *y.get_u_xx();
00548 
00549     *z.get_u_xxy() = f_uuu * square(*x.get_u_x())* *x.get_u_y()
00550                    + 2.0 * f_uu * *x.get_u_x() * *x.get_u_xy()
00551                    + f_uu * *x.get_u_y() * *x.get_u_xx()
00552                    + f_u * *x.get_u_xxy()
00553                    + f_vvv * square(*y.get_u_x())* *y.get_u_y()
00554                    + 2.0 * f_vv * *y.get_u_x() * *y.get_u_xy()
00555                    + f_vv * *y.get_u_xx() * *y.get_u_y()
00556                    + f_v * *y.get_u_xxy()
00557                    + f_uuv * square(*x.get_u_x()) * *y.get_u_y()
00558                    + 2.0* f_uuv * *x.get_u_x() * *x.get_u_y() * *y.get_u_x()
00559                    + f_uvv * *x.get_u_y()* square(*y.get_u_x())
00560                    + 2.0 * f_uvv * *x.get_u_x()* *y.get_u_x() * *y.get_u_y()
00561                    + f_uv * *x.get_u_xx() * *y.get_u_y()
00562                    + f_uv * *x.get_u_y() * *y.get_u_xx()
00563                    + 2.0* f_uv * *x.get_u_xy() * *y.get_u_x()
00564                    + 2.0* f_uv * *x.get_u_x() * *y.get_u_xy();
00565 
00566     *z.get_u_xyy() = f_uuu * square(*x.get_u_y())* *x.get_u_x()
00567                    + 2.0 * f_uu * *x.get_u_y() * *x.get_u_xy()
00568                    + f_uu * *x.get_u_x() * *x.get_u_yy()
00569                    + f_u * *x.get_u_xyy()
00570                    + f_vvv * square(*y.get_u_y())* *y.get_u_x()
00571                    + 2.0 * f_vv * *y.get_u_y() * *y.get_u_xy()
00572                    + f_vv * *y.get_u_yy() * *y.get_u_x()
00573                    + f_v * *y.get_u_xyy()
00574                    + f_uuv * square(*x.get_u_y()) * *y.get_u_x()
00575                    + 2.0* f_uuv * *x.get_u_y() * *x.get_u_x() * *y.get_u_y()
00576                    + f_uvv * *x.get_u_x()* square(*y.get_u_y())
00577                    + 2.0 * f_uvv * *x.get_u_y()* *y.get_u_y() * *y.get_u_x()
00578                    + f_uv * *x.get_u_yy() * *y.get_u_x()
00579                    + f_uv * *x.get_u_x() * *y.get_u_yy()
00580                    + 2.0* f_uv * *x.get_u_xy() * *y.get_u_y()
00581                    + 2.0* f_uv * *x.get_u_y() * *y.get_u_xy();
00582 
00583 
00584     *z.get_u_yyy() = f_uuu * cube(*x.get_u_y())
00585                    + 3.0 * f_uu * *x.get_u_y() * *x.get_u_yy()
00586                    + f_u * *x.get_u_yyy()
00587                    + f_vvv * cube(*y.get_u_y())
00588                    + 3.0 * f_vv * *y.get_u_y() * *y.get_u_yy()
00589                    + f_v * *y.get_u_yyy()
00590                    + 3.0 * f_uuv * square(*x.get_u_y()) * *y.get_u_y()
00591                    + 3.0 * f_uvv * *x.get_u_y()* square(*y.get_u_y())
00592                    + 3.0 * f_uv * *x.get_u_yy() * *y.get_u_y()
00593                    + 3.0 * f_uv * *x.get_u_y() * *y.get_u_yy();
00594 }
00595 
00600   df3_two_variable sqrt(const df3_two_variable& x)
00601   {
00602     df3_two_variable z;
00603     double u=::sqrt(*x.get_u());
00604     *z.get_u()=u;
00605     double xinv=1.0/(*x.get_u());
00606     double zp=0.5/u;
00607     double zp2=-0.5*zp*xinv;
00608     double zp3=-1.5*zp2*xinv;
00609 
00610 
00611     set_derivatives(z,x,u,zp,zp2,zp3);
00612 
00613     return z;
00614   }
00615 
00620   df3_two_variable atan(const df3_two_variable& x)
00621   {
00622     df3_two_variable z;
00623     double cx=value(x);
00624     double d=1.0/(1+square(cx));
00625     double d2=square(d);
00626     double u=::atan(cx);
00627     *z.get_u()=u;
00628     double zp=d;
00629     double zp2=-2.0*cx*d2;
00630     double zp3=-2.0*d2+8*cx*cx*d*d2;
00631 
00632     set_derivatives(z,x,u,zp,zp2,zp3);
00633     return z;
00634   }
00635 
00640   df3_two_variable square(const df3_two_variable& x)
00641   {
00642     df3_two_variable z;
00643     double u=value(x);
00644     *z.get_u()=u*u;
00645     double zp=2.0*u;
00646     double zp2=2.0;
00647     double zp3=0.0;
00648 
00649     set_derivatives(z,x,u,zp,zp2,zp3);
00650     return z;
00651   }
00652 
00657   df3_two_variable tan(const df3_two_variable& x)
00658   {
00659     df3_two_variable z;
00660     double u=::tan(*x.get_u());
00661     *z.get_u()=u;
00662     double v=1.0/::cos(*x.get_u());
00663     double w=::sin(*x.get_u());
00664     double v2=v*v;
00665     double zp=v2;
00666     double zp2=2.0*w*v2*v;
00667     double zp3=(4.0*w*w+2.0)*v2*v2;
00668 
00669     set_derivatives(z,x,u,zp,zp2,zp3);
00670     return z;
00671   }
00672 
00677   df3_two_variable sin(const df3_two_variable& x)
00678   {
00679     df3_two_variable z;
00680     double u=::sin(*x.get_u());
00681     *z.get_u()=u;
00682     double zp=::cos(*x.get_u());
00683     double zp2=-u;
00684     double zp3=-zp;
00685 
00686     set_derivatives(z,x,u,zp,zp2,zp3);
00687     return z;
00688   }
00689 
00694   df3_two_variable fabs(const df3_two_variable& v)
00695   {
00696     df3_two_variable z;
00697     if (value(v)>=0)
00698     {
00699       *z.get_u() = *v.get_u();
00700       *z.get_u_x() = *v.get_u_x();
00701       *z.get_u_y() = *v.get_u_y();
00702       *z.get_u_xx() = *v.get_u_xx();
00703       *z.get_u_xy() = *v.get_u_xy();
00704       *z.get_u_yy() = *v.get_u_yy();
00705       *z.get_u_xxx() = *v.get_u_xxx();
00706       *z.get_u_xxy() = *v.get_u_xxy();
00707       *z.get_u_xyy() = *v.get_u_xyy();
00708       *z.get_u_yyy() = *v.get_u_yyy();
00709     }
00710     else
00711     {
00712       *z.get_u() = -*v.get_u();
00713       *z.get_u_x() = -*v.get_u_x();
00714       *z.get_u_y() = -*v.get_u_y();
00715       *z.get_u_xx() = -*v.get_u_xx();
00716       *z.get_u_xy() = -*v.get_u_xy();
00717       *z.get_u_yy() = -*v.get_u_yy();
00718       *z.get_u_xxx() = -*v.get_u_xxx();
00719       *z.get_u_xxy() = -*v.get_u_xxy();
00720       *z.get_u_xyy() = -*v.get_u_xyy();
00721       *z.get_u_yyy() = -*v.get_u_yyy();
00722     }
00723 
00724     return z;
00725   }
00726 
00731   df3_two_variable log(const df3_two_variable& x)
00732   {
00733     df3_two_variable z;
00734     double u=::log(*x.get_u());
00735     *z.get_u()=u;
00736     double zp=1/(*x.get_u());
00737     double zp2=-zp*zp;
00738     double zp3=-2.0*zp*zp2;
00739 
00740     set_derivatives(z,x,u,zp,zp2,zp3);
00741     return z;
00742   }
00743 
00748   df3_two_variable exp(const df3_two_variable& x)
00749   {
00750     df3_two_variable z;
00751     double u=::exp(*x.get_u());
00752     *z.get_u()=u;
00753     double zp=u;
00754     double zp2=u;
00755     double zp3=u;
00756 
00757     set_derivatives(z,x,u,zp,zp2,zp3);
00758     return z;
00759   }
00760 
00765   df3_two_variable pow(const df3_two_variable& x,
00766                        const df3_two_variable& y)
00767   {
00768     df3_two_variable z;
00769     z=exp(y*log(x));
00770     return z;
00771   }
00772 
00777   df3_two_variable inv(const df3_two_variable& x)
00778   {
00779     df3_two_variable z;
00780     double xinv=1.0/(*x.get_u());
00781     *z.get_u()=xinv;
00782     double zp=-xinv*xinv;
00783     double zp2=-2.0*zp*xinv;
00784     double zp3=-3.0*zp2*xinv;
00785 
00786     set_derivatives(z,x,xinv,zp,zp2,zp3);
00787 
00788     return z;
00789   }
00790 
00795   df3_two_variable& df3_two_variable::operator = (const df3_two_variable& x)
00796   {
00797     *get_u() = *x.get_u();
00798     *get_u_x() = *x.get_u_x();
00799     *get_u_y() = *x.get_u_y();
00800     *get_u_xx() = *x.get_u_xx();
00801     *get_u_xy() = *x.get_u_xy();
00802     *get_u_yy() = *x.get_u_yy();
00803     *get_u_xxx() = *x.get_u_xxx();
00804     *get_u_xxy() = *x.get_u_xxy();
00805     *get_u_xyy() = *x.get_u_xyy();
00806     *get_u_yyy() = *x.get_u_yyy();
00807     return *this;
00808   }
00809 
00814   df3_two_variable& df3_two_variable::operator = (double x)
00815   {
00816     *get_u() = x;
00817     *get_u_x() =0.0;
00818     *get_u_y() =0.0;
00819     *get_u_xx() =0.0;
00820     *get_u_xy() =0.0;
00821     *get_u_yy() =0.0;
00822     *get_u_xxx() =0.0;
00823     *get_u_xxy() =0.0;
00824     *get_u_xyy() =0.0;
00825     *get_u_yyy() =0.0;
00826     return *this;
00827   }
00828 
00833   df3_two_variable operator * (const df3_two_variable& x,
00834     const df3_two_variable& y)
00835   {
00836     df3_two_variable z;
00837     double u= *x.get_u() * *y.get_u();
00838     *z.get_u() = u;
00839     double f_u=*y.get_u();
00840     double f_v=*x.get_u();
00841     double f_uu=0.0;
00842     double f_uv=1.0;
00843     double f_vv=0.0;
00844     double f_uuu=0.0;
00845     double f_uuv=0.0;
00846     double f_uvv=0.0;
00847     double f_vvv=0.0;
00848     set_derivatives(z,x,y,u,
00849       f_u, f_v,
00850       f_uu, f_uv, f_vv,
00851       f_uuu, f_uuv, f_uvv, f_vvv);
00852     return z;
00853   }
00854 
00859   df3_two_variable operator * (double x,
00860     const df3_two_variable& y)
00861   {
00862     df3_two_variable z;
00863     *z.get_u() = x *  *y.get_u();
00864     *z.get_u_x() = x * *y.get_u_x();
00865     *z.get_u_y() = x * *y.get_u_y();
00866     *z.get_u_xx() = x * *y.get_u_xx();
00867     *z.get_u_xy() = x * *y.get_u_xy();
00868     *z.get_u_yy() = x * *y.get_u_yy();
00869     *z.get_u_xxx() = x * *y.get_u_xxx();
00870     *z.get_u_xxy() = x * *y.get_u_xxy();
00871     *z.get_u_xyy() = x * *y.get_u_xyy();
00872     *z.get_u_yyy() = x * *y.get_u_yyy();
00873 
00874     return z;
00875   }
00876 
00881   df3_two_variable operator * (const df3_two_variable& y,
00882     double x)
00883   {
00884     df3_two_variable z;
00885     *z.get_u() = x *  *y.get_u();
00886     *z.get_u_x() = x * *y.get_u_x();
00887     *z.get_u_y() = x * *y.get_u_y();
00888     *z.get_u_xx() = x * *y.get_u_xx();
00889     *z.get_u_xy() = x * *y.get_u_xy();
00890     *z.get_u_yy() = x * *y.get_u_yy();
00891     *z.get_u_xxx() = x * *y.get_u_xxx();
00892     *z.get_u_xxy() = x * *y.get_u_xxy();
00893     *z.get_u_xyy() = x * *y.get_u_xyy();
00894     *z.get_u_yyy() = x * *y.get_u_yyy();
00895 
00896     return z;
00897   }
00898 
00903   df3_two_variable operator / (const df3_two_variable& x,
00904     double y)
00905   {
00906     double u=1/y;
00907     return x*u;
00908   }
00909 
00914   df3_two_variable operator / (const df3_two_variable& x,
00915     const df3_two_variable& y)
00916   {
00917     df3_two_variable u=inv(y);
00918     return x*u;
00919   }
00920 
00925   df3_two_variable operator / (const double x,
00926     const df3_two_variable& y)
00927   {
00928     df3_two_variable u=inv(y);
00929     return x*u;
00930   }
00931 
00936   df3_two_variable operator + (const double x,const df3_two_variable& y)
00937   {
00938     df3_two_variable z;
00939     *z.get_u() =  x + *y.get_u();
00940     *z.get_u_x() = *y.get_u_x();
00941     *z.get_u_y() = *y.get_u_y();
00942     *z.get_u_xx() = *y.get_u_xx();
00943     *z.get_u_xy() = *y.get_u_xy();
00944     *z.get_u_yy() = *y.get_u_yy();
00945     *z.get_u_xxx() = *y.get_u_xxx();
00946     *z.get_u_xxy() = *y.get_u_xxy();
00947     *z.get_u_xyy() = *y.get_u_xyy();
00948     *z.get_u_yyy() = *y.get_u_yyy();
00949     return z;
00950   }
00951 
00956   df3_two_variable operator + (const df3_two_variable& x,const double y)
00957   {
00958     df3_two_variable z;
00959     *z.get_u() =  *x.get_u() + y;
00960     *z.get_u_x() = *x.get_u_x();
00961     *z.get_u_y() = *x.get_u_y();
00962     *z.get_u_xx() = *x.get_u_xx();
00963     *z.get_u_xy() = *x.get_u_xy();
00964     *z.get_u_yy() = *x.get_u_yy();
00965     *z.get_u_xxx() = *x.get_u_xxx();
00966     *z.get_u_xxy() = *x.get_u_xxy();
00967     *z.get_u_xyy() = *x.get_u_xyy();
00968     *z.get_u_yyy() = *x.get_u_yyy();
00969     return z;
00970   }
00971 
00976   df3_two_variable operator + (const df3_two_variable& x,
00977     const df3_two_variable& y)
00978   {
00979     df3_two_variable z;
00980     *z.get_u() = *x.get_u() + *y.get_u();
00981     *z.get_u_x() = *x.get_u_x() + *y.get_u_x();
00982     *z.get_u_y() = *x.get_u_y()+*y.get_u_y();
00983     *z.get_u_xx() = *x.get_u_xx()+ *y.get_u_xx();
00984     *z.get_u_xy() = *x.get_u_xy()+ *y.get_u_xy();
00985     *z.get_u_yy() = *x.get_u_yy()+ *y.get_u_yy();
00986     *z.get_u_xxx() = *x.get_u_xxx()+ *y.get_u_xxx();
00987     *z.get_u_xxy() = *x.get_u_xxy()+ *y.get_u_xxy();
00988     *z.get_u_xyy() = *x.get_u_xyy()+ *y.get_u_xyy();
00989     *z.get_u_yyy() = *x.get_u_yyy()+ *y.get_u_yyy();
00990     return z;
00991   }
00992 
00997   df3_two_variable operator - (const df3_two_variable& x,
00998     const df3_two_variable& y)
00999   {
01000     df3_two_variable z;
01001     *z.get_u() = *x.get_u() - *y.get_u();
01002     *z.get_u_x() = *x.get_u_x()  - *y.get_u_x();
01003     *z.get_u_y() = *x.get_u_y() - *y.get_u_y();
01004     *z.get_u_xx() = *x.get_u_xx() - *y.get_u_xx();
01005     *z.get_u_xy() = *x.get_u_xy() - *y.get_u_xy();
01006     *z.get_u_yy() = *x.get_u_yy() - *y.get_u_yy();
01007     *z.get_u_xxx() = *x.get_u_xxx() - *y.get_u_xxx();
01008     *z.get_u_xxy() = *x.get_u_xxy() - *y.get_u_xxy();
01009     *z.get_u_xyy() = *x.get_u_xyy() - *y.get_u_xyy();
01010     *z.get_u_yyy() = *x.get_u_yyy() - *y.get_u_yyy();
01011     return z;
01012   }
01013 
01018   df3_two_variable operator - (double x,
01019     const df3_two_variable& y)
01020   {
01021     df3_two_variable z;
01022     *z.get_u() = x - *y.get_u();
01023     *z.get_u_x() = - *y.get_u_x();
01024     *z.get_u_y() = - *y.get_u_y();
01025     *z.get_u_xx() = - *y.get_u_xx();
01026     *z.get_u_xy() = - *y.get_u_xy();
01027     *z.get_u_yy() = - *y.get_u_yy();
01028     *z.get_u_xxx() = - *y.get_u_xxx();
01029     *z.get_u_xxy() = - *y.get_u_xxy();
01030     *z.get_u_xyy() = - *y.get_u_xyy();
01031     *z.get_u_yyy() = - *y.get_u_yyy();
01032     return z;
01033   }
01034 
01039   df3_two_variable operator - (const df3_two_variable& x,
01040     double y)
01041   {
01042     df3_two_variable z;
01043     *z.get_u() = *x.get_u()-y;
01044     *z.get_u_x() = *x.get_u_x();
01045     *z.get_u_y() = *x.get_u_y();
01046     *z.get_u_xx() = *x.get_u_xx();
01047     *z.get_u_xy() = *x.get_u_xy();
01048     *z.get_u_yy() = *x.get_u_yy();
01049     *z.get_u_xxx() = *x.get_u_xxx();
01050     *z.get_u_xxy() = *x.get_u_xxy();
01051     *z.get_u_xyy() = *x.get_u_xyy();
01052     *z.get_u_yyy() = *x.get_u_yyy();
01053     return z;
01054   }
01055 
01060 init_df3_two_variable::init_df3_two_variable(const df1b2variable& _v)
01061 {
01062   /*
01063   if (ind_var != 0)
01064   {
01065     cerr << " can only have 1 independent_variable in a reverse funnel"
01066            << endl;
01067     ad_exit(1);
01068   }
01069   */
01070   if (num_ind_var>1)
01071   {
01072     cerr << "can only have 2 independent_variables in df3_two_variable"
01073        " function" << endl;
01074     ad_exit(1);
01075   }
01076   else
01077   {
01078     ADUNCONST(df1b2variable,v)
01079 
01080     ind_var[num_ind_var++]=&v;
01081     *get_u() =  *v.get_u();
01082     switch(num_ind_var)
01083     {
01084     case 1:
01085       *get_u_x() = 1.0;
01086       *get_u_y() = 0.0;
01087       break;
01088     case 2:
01089       *get_u_x() = 0.0;
01090       *get_u_y() = 1.0;
01091       break;
01092     default:
01093       cerr << "illegal num_ind_var value of " << num_ind_var
01094            << " in  df3_two_variable function" << endl;
01095       ad_exit(1);
01096     }
01097     *get_u_xx() = 0.0;
01098     *get_u_xy() = 0.0;
01099     *get_u_yy() = 0.0;
01100     *get_u_xxx() = 0.0;
01101     *get_u_xxy() = 0.0;
01102     *get_u_xyy() = 0.0;
01103     *get_u_yyy() = 0.0;
01104   }
01105 }
01106 
01111   init_df3_two_variable::init_df3_two_variable(double v)
01112   {
01113     *get_u() =  v;
01114     *get_u_x() = 0.0;
01115     *get_u_y() = 0.0;
01116     *get_u_xx() = 0.0;
01117     *get_u_xy() = 0.0;
01118     *get_u_yy() = 0.0;
01119     *get_u_xxx() = 0.0;
01120     *get_u_xxy() = 0.0;
01121     *get_u_xyy() = 0.0;
01122     *get_u_yyy() = 0.0;
01123   }
01124 
01125   df3_two_variable::df3_two_variable(void)
01126   {
01127   }
01128 
01133 df3_two_matrix choleski_decomp(const df3_two_matrix& MM)
01134 {
01135   // kludge to deal with constantness
01136   df3_two_matrix & M= (df3_two_matrix &) MM;
01137   int rmin=M.indexmin();
01138   int cmin=M(rmin).indexmin();
01139   int rmax=M.indexmax();
01140   int cmax=M(rmin).indexmax();
01141   if (rmin !=1 || cmin !=1)
01142   {
01143     cerr << "minimum row and column inidices must equal 1 in "
01144       "df1b2matrix choleski_decomp(const df3_two_atrix& MM)"
01145          << endl;
01146     ad_exit(1);
01147   }
01148   if (rmax !=cmax)
01149   {
01150     cerr << "Error in df1b2matrix choleski_decomp(const df3_two_matrix& MM)"
01151       " Matrix not square" << endl;
01152     ad_exit(1);
01153   }
01154 
01155   int n=rmax-rmin+1;
01156   df3_two_matrix L(1,n,1,n);
01157 #ifndef SAFE_INITIALIZE
01158     L.initialize();
01159 #endif
01160 
01161   int i,j,k;
01162   df3_two_variable tmp;
01163 
01164     if (value(M(1,1))<=0)
01165     {
01166       cerr << "Error matrix not positive definite in choleski_decomp"
01167         <<endl;
01168       ad_exit(1);
01169     }
01170 
01171   L(1,1)=sqrt(M(1,1));
01172   for (i=2;i<=n;i++)
01173   {
01174     L(i,1)=M(i,1)/L(1,1);
01175   }
01176 
01177   for (i=2;i<=n;i++)
01178   {
01179     for (j=2;j<=i-1;j++)
01180     {
01181       tmp=M(i,j);
01182       for (k=1;k<=j-1;k++)
01183       {
01184         tmp-=L(i,k)*L(j,k);
01185       }
01186       L(i,j)=tmp/L(j,j);
01187     }
01188     tmp=M(i,i);
01189     for (k=1;k<=i-1;k++)
01190     {
01191       tmp-=L(i,k)*L(i,k);
01192     }
01193 
01194     if (value(tmp)<=0)
01195     {
01196       cerr << "Error matrix not positive definite in choleski_decomp"
01197         <<endl;
01198       ad_exit(1);
01199     }
01200 
01201     L(i,i)=sqrt(tmp);
01202   }
01203 
01204   return L;
01205 }
01206 
01211 df1b2matrix& df1b2matrix::operator = (const df3_two_matrix& M)
01212 {
01213   int rmin=M.indexmin();
01214   int rmax=M.indexmax();
01215   if (rmin != indexmin() || rmax != indexmax())
01216   {
01217     cerr << "unequal shape in "
01218      "df1b2matrix& df1b2matrix::operator = (const df3_two_matrix& M)"
01219       << endl;
01220     ad_exit(1);
01221   }
01222 
01223   for (int i=rmin;i<=rmax;i++)
01224   {
01225     (*this)(i)=M(i);
01226   }
01227   return *this;
01228 }
01229 
01234 df1b2vector& df1b2vector::operator = (const df3_two_vector& M)
01235 {
01236   int rmin=M.indexmin();
01237   int rmax=M.indexmax();
01238   if (rmin != indexmin() || rmax != indexmax())
01239   {
01240     cerr << "unequal shape in "
01241      "df1b2vector& df1b2vector::operator = (const df3_two_vector& M)"
01242       << endl;
01243     ad_exit(1);
01244   }
01245 
01246   for (int i=rmin;i<=rmax;i++)
01247   {
01248     (*this)(i)=M(i);
01249   }
01250   return *this;
01251 }
01252 
01257 df1b2variable& df1b2variable::operator = (const df3_two_variable& v)
01258 {
01259   df1b2variable * px=df3_two_variable::ind_var[0];
01260   df1b2variable * py=df3_two_variable::ind_var[1];
01261   df3_two_variable::num_ind_var=0;
01262   df3_two_variable::ind_var[0]=0;
01263   df3_two_variable::ind_var[1]=0;
01264   //df1b2variable * px=0;
01265   double  dfx= *v.get_u_x();
01266   double  dfy= *v.get_u_y();
01267   double * xd=px->get_u_dot();
01268   double * yd=py->get_u_dot();
01269   double * zd=get_u_dot();
01270   *get_u()=*v.get_u();
01271   for (unsigned int i=0;i<df1b2variable::nvar;i++)
01272   {
01273     *zd++ = dfx * *xd++ + dfy * *yd++;
01274   }
01275 
01276   f1b2gradlist->write_pass1(px,py,this,
01277     *v.get_u_x(),
01278     *v.get_u_y(),
01279     *v.get_u_xx(),
01280     *v.get_u_xy(),
01281     *v.get_u_yy(),
01282     *v.get_u_xxx(),
01283     *v.get_u_xxy(),
01284     *v.get_u_xyy(),
01285     *v.get_u_yyy());
01286   return *this;
01287 }
01288 
01293 df1b2variable div(const df1b2variable& x,const df1b2variable& y)
01294 {
01295   df1b2variable z;
01296   double xu=*x.get_u();
01297   double yu=*y.get_u();
01298   double yinv=1.0/yu;
01299   *z.get_u()=xu*yinv;
01300 
01301   double dfx= yinv;
01302   double dfy= -xu*yinv*yinv;
01303   double dfxx= 0.0;
01304   double dfxy=-yinv*yinv;
01305   double dfyy=2.0*xu*yinv*yinv*yinv;
01306   double dfxxx= 0.0;
01307   double dfxxy= 0.0;
01308   double dfxyy=2.0*yinv*yinv*yinv;
01309   double dfyyy=-6.0*xu*yinv*yinv*yinv*yinv;
01310 
01311   double * xd=x.get_u_dot();
01312   double * yd=y.get_u_dot();
01313   double * zd=z.get_u_dot();
01314 
01315   for (unsigned int i=0;i<df1b2variable::nvar;i++)
01316   {
01317     *zd++ = dfx * *xd++ + dfy * *yd++;
01318   }
01319 
01320   f1b2gradlist->write_pass1(&x,&y,&z,
01321     dfx,
01322     dfy,
01323     dfxx,dfxy,dfyy,
01324     dfxxx,dfxxy,dfxyy,dfyyy);
01325   return z;
01326 }
01327 
01332 df1b2variable mypow(const df1b2variable& x,double y)
01333 {
01334   df1b2variable z;
01335   double xu=*x.get_u();
01336   *z.get_u()=::pow(xu,y);
01337 
01338   double dfx= y*::pow(xu,y-1.0);
01339   double dfxx= y*(y-1.0)*::pow(xu,y-2.0);
01340   double dfxxx= y*(y-1.0)*(y-2.0)*::pow(xu,y-3.0);
01341   double * xd=x.get_u_dot();
01342   double * zd=z.get_u_dot();
01343 
01344   for (unsigned int i=0;i<df1b2variable::nvar;i++)
01345   {
01346     *zd++ = dfx * *xd++ ;
01347   }
01348 
01349   f1b2gradlist->write_pass1(&x,&z,
01350     dfx,
01351     dfxx,
01352     dfxxx);
01353   return z;
01354 }