ADMB Documentation  11.1.2382
 All Classes Files Functions Variables Typedefs Friends Defines
df32fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df32fun.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012   df1b2variable * df3_two_variable::ind_var[2];
00013   int df3_two_variable::num_ind_var=0;
00014 
00019   df3_two_variable::df3_two_variable(const df3_two_variable& x)
00020   {
00021     v[0]=x.v[0];
00022     v[1]=x.v[1];
00023     v[2]=x.v[2];
00024     v[3]=x.v[3];
00025     v[4]=x.v[4];
00026     v[5]=x.v[5];
00027     v[6]=x.v[6];
00028     v[7]=x.v[7];
00029     v[8]=x.v[8];
00030     v[9]=x.v[9];
00031   }
00032 
00037  df3_two_vector::df3_two_vector(const df3_two_vector& m2)
00038  {
00039    index_min=m2.index_min;
00040    index_max=m2.index_max;
00041    shape=m2.shape;
00042    if (shape)
00043    {
00044      (shape->ncopies)++;
00045    }
00046    v = m2.v;
00047  }
00048 
00053  df3_two_vector::~df3_two_vector()
00054  {
00055    deallocate();
00056  }
00057 
00062  void df3_two_vector::deallocate(void)
00063  {
00064    if(shape)
00065    {
00066      if (shape->ncopies)
00067      {
00068        (shape->ncopies)--;
00069      }
00070      else
00071      {
00072        v = (df3_two_variable*) (shape->trueptr);
00073        delete [] v;
00074        v = NULL;
00075        delete shape;
00076        shape=0;
00077      }
00078    }
00079  }
00080 
00085  dvector value(const df3_two_vector& v)
00086  {
00087    int mmin=v.indexmin();
00088    int mmax=v.indexmax();
00089    dvector cv(mmin,mmax);
00090    for (int i=mmin;i<=mmax;i++)
00091    {
00092      cv(i)=value(v(i));
00093    }
00094    return cv;
00095  }
00096 
00101   void df3_two_vector::initialize(void)
00102   {
00103     int mmin=indexmin();
00104     int mmax=indexmax();
00105     for (int i=mmin;i<=mmax;i++)
00106     {
00107       (*this)(i)=0.0;
00108     }
00109   }
00110 
00115   df3_two_vector::df3_two_vector(void)
00116   {
00117     allocate();
00118   }
00119 
00124   df3_two_vector::df3_two_vector(int min,int max)
00125   {
00126     allocate(min,max);
00127   }
00128 
00133   void df3_two_vector::allocate(int min,int max)
00134   {
00135     index_min=min;
00136     index_max=max;
00137     v=new df3_two_variable[max-min+1];
00138     if (v==0)
00139     {
00140       cerr << "error allocating memory in df3_two_vector" << endl;
00141       ad_exit(1);
00142     }
00143     if ( (shape=new vector_shapex(min,max,v)) == NULL)
00144     {
00145       cerr << "Error trying to allocate memory for df3_two_vector"
00146            << endl;;
00147       ad_exit(1);
00148     }
00149     v-=min;
00150   }
00151 
00156   void df3_two_vector::allocate(void)
00157   {
00158     index_min=0;
00159     index_max=-1;
00160     v=0;
00161     shape=0;
00162   }
00163 
00168  dmatrix value(const df3_two_matrix& v)
00169  {
00170    int rmin=v.indexmin();
00171    int rmax=v.indexmax();
00172    dmatrix cm(rmin,rmax);
00173    for (int i=rmin;i<=rmax;i++)
00174    {
00175      int cmin=v(i).indexmin();
00176      int cmax=v(i).indexmax();
00177      cm(i).allocate(cmin,cmax);
00178      for (int j=cmin;j<=cmax;j++)
00179      {
00180        cm(i,j)=value(v(i,j));
00181      }
00182    }
00183    return cm;
00184  }
00185 
00190  df3_two_matrix::df3_two_matrix(const df3_two_matrix& m2)
00191  {
00192    index_min=m2.index_min;
00193    index_max=m2.index_max;
00194    shape=m2.shape;
00195    if (shape)
00196    {
00197      (shape->ncopies)++;
00198    }
00199    v = m2.v;
00200  }
00201 
00206  df3_two_matrix::~df3_two_matrix()
00207  {
00208    deallocate();
00209  }
00210 
00215  void df3_two_matrix::deallocate(void)
00216  {
00217    if (shape)
00218    {
00219      if (shape->ncopies)
00220      {
00221        (shape->ncopies)--;
00222      }
00223      else
00224      {
00225        v = (df3_two_vector*) (shape->get_pointer());
00226        delete [] v;
00227        v=0;
00228        delete shape;
00229        shape=0;
00230      }
00231    }
00232  }
00233 
00238   void df3_two_matrix::initialize(void)
00239   {
00240     int mmin=indexmin();
00241     int mmax=indexmax();
00242     for (int i=mmin;i<=mmax;i++)
00243     {
00244       (*this)(i).initialize();
00245     }
00246   }
00247 
00252   df3_two_matrix::df3_two_matrix(int rmin,int rmax,int cmin,int cmax)
00253   {
00254     index_min=rmin;
00255     index_max=rmax;
00256     v=new df3_two_vector[rmax-rmin+1];
00257     if (v==0)
00258     {
00259       cerr << "error allocating memory in df3_two_matrix" << endl;
00260       ad_exit(1);
00261     }
00262     if ( (shape=new mat_shapex(v)) == NULL)
00263     {
00264       cerr << "Error trying to allocate memory for df3_two_vector"
00265            << endl;;
00266     }
00267     v-=rmin;
00268 
00269     for (int i=rmin;i<=rmax;i++)
00270     {
00271       v[i].allocate(cmin,cmax);
00272     }
00273   }
00274 
00279   df3_two_variable& df3_two_variable::operator -= (const df3_two_variable& v)
00280   {
00281     *get_u() -= *v.get_u();
00282 
00283     *get_u_x() -= *v.get_u_x();
00284     *get_u_y() -= *v.get_u_y();
00285     *get_u_xx() -= *v.get_u_xx();
00286     *get_u_xy() -= *v.get_u_xy();
00287     *get_u_yy() -= *v.get_u_yy();
00288     *get_u_xxx() -= *v.get_u_xxx();
00289     *get_u_xxy() -= *v.get_u_xxy();
00290     *get_u_xyy() -= *v.get_u_xyy();
00291     *get_u_yyy() -= *v.get_u_yyy();
00292     return *this;
00293   }
00294 
00299 df3_two_variable operator-(const df3_two_variable& v)
00300 {
00301   df3_two_variable z;
00302 
00303   *z.get_u() = - *v.get_u();
00304 
00305   *z.get_u_x() = -(*v.get_u_x());
00306   *z.get_u_y() = -(*v.get_u_y());
00307   *z.get_u_xx() = -(*v.get_u_xx());
00308   *z.get_u_xy() = -(*v.get_u_xy());
00309   *z.get_u_yy() = -(*v.get_u_yy());
00310   *z.get_u_xxx() = -(*v.get_u_xxx());
00311   *z.get_u_xxy() = -(*v.get_u_xxy());
00312   *z.get_u_xyy() = -(*v.get_u_xyy());
00313   *z.get_u_yyy() = -(*v.get_u_yyy());
00314 
00315   return z;
00316 }
00317 
00322   df3_two_variable& df3_two_variable::operator += (const df3_two_variable& v)
00323   {
00324     *get_u() += *v.get_u();
00325     *get_u_x() += *v.get_u_x();
00326     *get_u_y() += *v.get_u_y();
00327     *get_u_xx() += *v.get_u_xx();
00328     *get_u_xy() += *v.get_u_xy();
00329     *get_u_yy() += *v.get_u_yy();
00330     *get_u_xxx() += *v.get_u_xxx();
00331     *get_u_xxy() += *v.get_u_xxy();
00332     *get_u_xyy() += *v.get_u_xyy();
00333     *get_u_yyy() += *v.get_u_yyy();
00334 
00335     return *this;
00336   }
00337 
00342   df3_two_variable& df3_two_variable::operator += (double v)
00343   {
00344     *get_u() += v;
00345     return *this;
00346   }
00347 
00352   df3_two_variable& df3_two_variable::operator -= (double v)
00353   {
00354     *get_u() -= v;
00355     return *this;
00356   }
00357 
00362   df3_two_variable& df3_two_variable::operator *= (const df3_two_variable& v)
00363   {
00364     df3_two_variable x=*this * v;
00365     *this=x;
00366     return *this;
00367   }
00368 
00373   df3_two_variable& df3_two_variable::operator *= (double v)
00374   {
00375     *get_u() *= v;
00376     *get_u_x() *= v;
00377     *get_u_y() *= v;
00378     *get_u_xx() *= v;
00379     *get_u_xy() *= v;
00380     *get_u_yy() *= v;
00381     *get_u_xxx() *= v;
00382     *get_u_xxy() *= v;
00383     *get_u_xyy() *= v;
00384     *get_u_yyy() *= v;
00385     return *this;
00386   }
00387 
00392   df3_two_variable& df3_two_variable::operator /= (const df3_two_variable& y)
00393   {
00394     df3_two_variable x=*this / y;
00395     *this=x;
00396     return *this;
00397   }
00398 
00403 int operator <(const df3_two_variable & x, double n)
00404 {
00405    return value(x) < n;
00406 }
00407 
00412 int operator >(const df3_two_variable & x, double n)
00413 {
00414    return value(x) > n;
00415 }
00416 
00421 int operator >=(const df3_two_variable & x, double n)
00422 {
00423    return value(x) >= n;
00424 }
00425 
00430 int operator ==(const df3_two_variable & x, const df3_two_variable & n)
00431 {
00432    return value(x) == value(n);
00433 }
00434 
00439 int operator ==(const df3_two_variable & x, double n)
00440 {
00441    return value(x) == n;
00442 }
00443 
00448 int operator ==(double x, const df3_two_variable & n)
00449 {
00450    return x == value(n);
00451 }
00452 
00457 int operator <(const df3_two_variable & x, const df3_two_variable & n)
00458 {
00459    return value(x) < value(n);
00460 }
00461 
00466 int operator >(const df3_two_variable & x, const df3_two_variable & n)
00467 {
00468    return value(x) > value(n);
00469 }
00470 
00475 void set_derivatives( df3_two_variable& z,const df3_two_variable& x,double u,
00476   double zp,double zp2,double zp3)
00477 {
00478     //*z.get_u() = u;
00479 
00480     *z.get_u_x() = zp* *x.get_u_x();
00481 
00482     *z.get_u_y() = zp* *x.get_u_y();
00483 
00484     *z.get_u_xx() = zp2 * square(*x.get_u_x())
00485                    + zp * *x.get_u_xx();
00486 
00487     *z.get_u_xy() = zp2 * *x.get_u_x() * *x.get_u_y()
00488                    + zp * *x.get_u_xy();
00489 
00490     *z.get_u_yy() = zp2 * square(*x.get_u_y())
00491                    + zp * *x.get_u_yy();
00492 
00493     *z.get_u_xxx() = zp3 * cube(*x.get_u_x())
00494                    + 3.0 * zp2 * *x.get_u_x() * *x.get_u_xx()
00495                    + zp * *x.get_u_xxx();
00496 
00497     *z.get_u_xxy() = zp3 * square(*x.get_u_x())* *x.get_u_y()
00498                    + 2.0* zp2 * *x.get_u_x() * *x.get_u_xy()
00499                    + zp2 * *x.get_u_y() * *x.get_u_xx()
00500                    + zp * *x.get_u_xxy();
00501 
00502     *z.get_u_xyy() = zp3 * square(*x.get_u_y())* *x.get_u_x()
00503                    + 2.0* zp2 * *x.get_u_y() * *x.get_u_xy()
00504                    + zp2 * *x.get_u_x() * *x.get_u_yy()
00505                    + zp * *x.get_u_xyy();
00506 
00507 
00508     *z.get_u_yyy() = zp3 * cube(*x.get_u_y())
00509                    + 3.0 * zp2 * *x.get_u_y() * *x.get_u_yy()
00510                    + zp * *x.get_u_yyy();
00511 }
00512 
00517 void set_derivatives( df3_two_variable& z, const df3_two_variable& x,
00518   const df3_two_variable& y, double u,
00519   double f_u,double f_v,double f_uu,double f_uv,double f_vv,double f_uuu,
00520   double f_uuv,double f_uvv,double f_vvv)
00521 {
00522     *z.get_u() = u;
00523 
00524     *z.get_u_x() = f_u* *x.get_u_x()
00525                  + f_v* *y.get_u_x();
00526 
00527     *z.get_u_y() = f_u* *x.get_u_y()
00528                  + f_v* *y.get_u_y();
00529 
00530     *z.get_u_xx() = f_uu * square(*x.get_u_x())
00531                   + f_u  * *x.get_u_xx()
00532                   + f_vv * square(*y.get_u_x())
00533                   + f_v  * *y.get_u_xx()
00534             + 2.0 * f_uv * *x.get_u_x() * *y.get_u_x();
00535 
00536     *z.get_u_xy() = f_uu * *x.get_u_x() * *x.get_u_y()
00537                   + f_u  * *x.get_u_xy()
00538                   + f_vv *  *y.get_u_x() * *y.get_u_y()
00539                   + f_v  * *y.get_u_xy()
00540                   + f_uv * (*x.get_u_x() * *y.get_u_y()
00541                          +  *x.get_u_y() * *y.get_u_x());
00542 
00543     *z.get_u_yy() = f_uu * square(*x.get_u_y())
00544                   + f_u  * *x.get_u_yy()
00545                   + f_vv * square(*y.get_u_y())
00546                   + f_v  * *y.get_u_yy()
00547             + 2.0 * f_uv * *x.get_u_y() * *y.get_u_y();
00548 
00549 
00550     *z.get_u_xxx() = f_uuu * cube(*x.get_u_x())
00551                    + 3.0 * f_uu * *x.get_u_x() * *x.get_u_xx()
00552                    + f_u * *x.get_u_xxx()
00553                    + f_vvv * cube(*y.get_u_x())
00554                    + 3.0 * f_vv * *y.get_u_x() * *y.get_u_xx()
00555                    + f_v * *y.get_u_xxx()
00556                    + 3.0 * f_uuv * square(*x.get_u_x()) * *y.get_u_x()
00557                    + 3.0 * f_uvv * *x.get_u_x()* square(*y.get_u_x())
00558                    + 3.0* f_uv * *x.get_u_xx() * *y.get_u_x()
00559                    + 3.0* f_uv * *x.get_u_x() * *y.get_u_xx();
00560 
00561     *z.get_u_xxy() = f_uuu * square(*x.get_u_x())* *x.get_u_y()
00562                    + 2.0 * f_uu * *x.get_u_x() * *x.get_u_xy()
00563                    + f_uu * *x.get_u_y() * *x.get_u_xx()
00564                    + f_u * *x.get_u_xxy()
00565                    + f_vvv * square(*y.get_u_x())* *y.get_u_y()
00566                    + 2.0 * f_vv * *y.get_u_x() * *y.get_u_xy()
00567                    + f_vv * *y.get_u_xx() * *y.get_u_y()
00568                    + f_v * *y.get_u_xxy()
00569                    + f_uuv * square(*x.get_u_x()) * *y.get_u_y()
00570                    + 2.0* f_uuv * *x.get_u_x() * *x.get_u_y() * *y.get_u_x()
00571                    + f_uvv * *x.get_u_y()* square(*y.get_u_x())
00572                    + 2.0 * f_uvv * *x.get_u_x()* *y.get_u_x() * *y.get_u_y()
00573                    + f_uv * *x.get_u_xx() * *y.get_u_y()
00574                    + f_uv * *x.get_u_y() * *y.get_u_xx()
00575                    + 2.0* f_uv * *x.get_u_xy() * *y.get_u_x()
00576                    + 2.0* f_uv * *x.get_u_x() * *y.get_u_xy();
00577 
00578     *z.get_u_xyy() = f_uuu * square(*x.get_u_y())* *x.get_u_x()
00579                    + 2.0 * f_uu * *x.get_u_y() * *x.get_u_xy()
00580                    + f_uu * *x.get_u_x() * *x.get_u_yy()
00581                    + f_u * *x.get_u_xyy()
00582                    + f_vvv * square(*y.get_u_y())* *y.get_u_x()
00583                    + 2.0 * f_vv * *y.get_u_y() * *y.get_u_xy()
00584                    + f_vv * *y.get_u_yy() * *y.get_u_x()
00585                    + f_v * *y.get_u_xyy()
00586                    + f_uuv * square(*x.get_u_y()) * *y.get_u_x()
00587                    + 2.0* f_uuv * *x.get_u_y() * *x.get_u_x() * *y.get_u_y()
00588                    + f_uvv * *x.get_u_x()* square(*y.get_u_y())
00589                    + 2.0 * f_uvv * *x.get_u_y()* *y.get_u_y() * *y.get_u_x()
00590                    + f_uv * *x.get_u_yy() * *y.get_u_x()
00591                    + f_uv * *x.get_u_x() * *y.get_u_yy()
00592                    + 2.0* f_uv * *x.get_u_xy() * *y.get_u_y()
00593                    + 2.0* f_uv * *x.get_u_y() * *y.get_u_xy();
00594 
00595 
00596     *z.get_u_yyy() = f_uuu * cube(*x.get_u_y())
00597                    + 3.0 * f_uu * *x.get_u_y() * *x.get_u_yy()
00598                    + f_u * *x.get_u_yyy()
00599                    + f_vvv * cube(*y.get_u_y())
00600                    + 3.0 * f_vv * *y.get_u_y() * *y.get_u_yy()
00601                    + f_v * *y.get_u_yyy()
00602                    + 3.0 * f_uuv * square(*x.get_u_y()) * *y.get_u_y()
00603                    + 3.0 * f_uvv * *x.get_u_y()* square(*y.get_u_y())
00604                    + 3.0 * f_uv * *x.get_u_yy() * *y.get_u_y()
00605                    + 3.0 * f_uv * *x.get_u_y() * *y.get_u_yy();
00606 }
00607 
00612   df3_two_variable sqrt(const df3_two_variable& x)
00613   {
00614     df3_two_variable z;
00615     double u=::sqrt(*x.get_u());
00616     *z.get_u()=u;
00617     double xinv=1.0/(*x.get_u());
00618     double zp=0.5/u;
00619     double zp2=-0.5*zp*xinv;
00620     double zp3=-1.5*zp2*xinv;
00621 
00622 
00623     set_derivatives(z,x,u,zp,zp2,zp3);
00624 
00625     return z;
00626   }
00627 
00632   df3_two_variable atan(const df3_two_variable& x)
00633   {
00634     df3_two_variable z;
00635     double cx=value(x);
00636     double d=1.0/(1+square(cx));
00637     double d2=square(d);
00638     double u=::atan(cx);
00639     *z.get_u()=u;
00640     double zp=d;
00641     double zp2=-2.0*cx*d2;
00642     double zp3=-2.0*d2+8*cx*cx*d*d2;
00643 
00644     set_derivatives(z,x,u,zp,zp2,zp3);
00645     return z;
00646   }
00647 
00652   df3_two_variable square(const df3_two_variable& x)
00653   {
00654     df3_two_variable z;
00655     double u=value(x);
00656     *z.get_u()=u*u;
00657     double zp=2.0*u;
00658     double zp2=2.0;
00659     double zp3=0.0;
00660 
00661     set_derivatives(z,x,u,zp,zp2,zp3);
00662     return z;
00663   }
00664 
00669   df3_two_variable tan(const df3_two_variable& x)
00670   {
00671     df3_two_variable z;
00672     double u=::tan(*x.get_u());
00673     *z.get_u()=u;
00674     double v=1.0/::cos(*x.get_u());
00675     double w=::sin(*x.get_u());
00676     double v2=v*v;
00677     double zp=v2;
00678     double zp2=2.0*w*v2*v;
00679     double zp3=(4.0*w*w+2.0)*v2*v2;
00680 
00681     set_derivatives(z,x,u,zp,zp2,zp3);
00682     return z;
00683   }
00684 
00689   df3_two_variable sin(const df3_two_variable& x)
00690   {
00691     df3_two_variable z;
00692     double u=::sin(*x.get_u());
00693     *z.get_u()=u;
00694     double zp=::cos(*x.get_u());
00695     double zp2=-u;
00696     double zp3=-zp;
00697 
00698     set_derivatives(z,x,u,zp,zp2,zp3);
00699     return z;
00700   }
00701 
00706   df3_two_variable fabs(const df3_two_variable& v)
00707   {
00708     df3_two_variable z;
00709     if (value(v)>=0)
00710     {
00711       *z.get_u() = *v.get_u();
00712       *z.get_u_x() = *v.get_u_x();
00713       *z.get_u_y() = *v.get_u_y();
00714       *z.get_u_xx() = *v.get_u_xx();
00715       *z.get_u_xy() = *v.get_u_xy();
00716       *z.get_u_yy() = *v.get_u_yy();
00717       *z.get_u_xxx() = *v.get_u_xxx();
00718       *z.get_u_xxy() = *v.get_u_xxy();
00719       *z.get_u_xyy() = *v.get_u_xyy();
00720       *z.get_u_yyy() = *v.get_u_yyy();
00721     }
00722     else
00723     {
00724       *z.get_u() = -*v.get_u();
00725       *z.get_u_x() = -*v.get_u_x();
00726       *z.get_u_y() = -*v.get_u_y();
00727       *z.get_u_xx() = -*v.get_u_xx();
00728       *z.get_u_xy() = -*v.get_u_xy();
00729       *z.get_u_yy() = -*v.get_u_yy();
00730       *z.get_u_xxx() = -*v.get_u_xxx();
00731       *z.get_u_xxy() = -*v.get_u_xxy();
00732       *z.get_u_xyy() = -*v.get_u_xyy();
00733       *z.get_u_yyy() = -*v.get_u_yyy();
00734     }
00735 
00736     return z;
00737   }
00738 
00743   df3_two_variable log(const df3_two_variable& x)
00744   {
00745     df3_two_variable z;
00746     double u=::log(*x.get_u());
00747     *z.get_u()=u;
00748     double zp=1/(*x.get_u());
00749     double zp2=-zp*zp;
00750     double zp3=-2.0*zp*zp2;
00751 
00752     set_derivatives(z,x,u,zp,zp2,zp3);
00753     return z;
00754   }
00755 
00760   df3_two_variable exp(const df3_two_variable& x)
00761   {
00762     df3_two_variable z;
00763     double u=::exp(*x.get_u());
00764     *z.get_u()=u;
00765     double zp=u;
00766     double zp2=u;
00767     double zp3=u;
00768 
00769     set_derivatives(z,x,u,zp,zp2,zp3);
00770     return z;
00771   }
00772 
00777   df3_two_variable pow(const df3_two_variable& x,
00778                        const df3_two_variable& y)
00779   {
00780     df3_two_variable z;
00781     z=exp(y*log(x));
00782     return z;
00783   }
00784 
00789   df3_two_variable inv(const df3_two_variable& x)
00790   {
00791     df3_two_variable z;
00792     double xinv=1.0/(*x.get_u());
00793     *z.get_u()=xinv;
00794     double zp=-xinv*xinv;
00795     double zp2=-2.0*zp*xinv;
00796     double zp3=-3.0*zp2*xinv;
00797 
00798     set_derivatives(z,x,xinv,zp,zp2,zp3);
00799 
00800     return z;
00801   }
00802 
00807   df3_two_variable& df3_two_variable::operator = (const df3_two_variable& x)
00808   {
00809     *get_u() = *x.get_u();
00810     *get_u_x() = *x.get_u_x();
00811     *get_u_y() = *x.get_u_y();
00812     *get_u_xx() = *x.get_u_xx();
00813     *get_u_xy() = *x.get_u_xy();
00814     *get_u_yy() = *x.get_u_yy();
00815     *get_u_xxx() = *x.get_u_xxx();
00816     *get_u_xxy() = *x.get_u_xxy();
00817     *get_u_xyy() = *x.get_u_xyy();
00818     *get_u_yyy() = *x.get_u_yyy();
00819     return *this;
00820   }
00821 
00826   df3_two_variable& df3_two_variable::operator = (double x)
00827   {
00828     *get_u() = x;
00829     *get_u_x() =0.0;
00830     *get_u_y() =0.0;
00831     *get_u_xx() =0.0;
00832     *get_u_xy() =0.0;
00833     *get_u_yy() =0.0;
00834     *get_u_xxx() =0.0;
00835     *get_u_xxy() =0.0;
00836     *get_u_xyy() =0.0;
00837     *get_u_yyy() =0.0;
00838     return *this;
00839   }
00840 
00845   df3_two_variable operator * (const df3_two_variable& x,
00846     const df3_two_variable& y)
00847   {
00848     df3_two_variable z;
00849     double u= *x.get_u() * *y.get_u();
00850     *z.get_u() = u;
00851     double f_u=*y.get_u();
00852     double f_v=*x.get_u();
00853     double f_uu=0.0;
00854     double f_uv=1.0;
00855     double f_vv=0.0;
00856     double f_uuu=0.0;
00857     double f_uuv=0.0;
00858     double f_uvv=0.0;
00859     double f_vvv=0.0;
00860     set_derivatives(z,x,y,u,
00861       f_u, f_v,
00862       f_uu, f_uv, f_vv,
00863       f_uuu, f_uuv, f_uvv, f_vvv);
00864     return z;
00865   }
00866 
00871   df3_two_variable operator * (double x,
00872     const df3_two_variable& y)
00873   {
00874     df3_two_variable z;
00875     *z.get_u() = x *  *y.get_u();
00876     *z.get_u_x() = x * *y.get_u_x();
00877     *z.get_u_y() = x * *y.get_u_y();
00878     *z.get_u_xx() = x * *y.get_u_xx();
00879     *z.get_u_xy() = x * *y.get_u_xy();
00880     *z.get_u_yy() = x * *y.get_u_yy();
00881     *z.get_u_xxx() = x * *y.get_u_xxx();
00882     *z.get_u_xxy() = x * *y.get_u_xxy();
00883     *z.get_u_xyy() = x * *y.get_u_xyy();
00884     *z.get_u_yyy() = x * *y.get_u_yyy();
00885 
00886     return z;
00887   }
00888 
00893   df3_two_variable operator * (const df3_two_variable& y,
00894     double x)
00895   {
00896     df3_two_variable z;
00897     *z.get_u() = x *  *y.get_u();
00898     *z.get_u_x() = x * *y.get_u_x();
00899     *z.get_u_y() = x * *y.get_u_y();
00900     *z.get_u_xx() = x * *y.get_u_xx();
00901     *z.get_u_xy() = x * *y.get_u_xy();
00902     *z.get_u_yy() = x * *y.get_u_yy();
00903     *z.get_u_xxx() = x * *y.get_u_xxx();
00904     *z.get_u_xxy() = x * *y.get_u_xxy();
00905     *z.get_u_xyy() = x * *y.get_u_xyy();
00906     *z.get_u_yyy() = x * *y.get_u_yyy();
00907 
00908     return z;
00909   }
00910 
00915   df3_two_variable operator / (const df3_two_variable& x,
00916     double y)
00917   {
00918     double u=1/y;
00919     return x*u;
00920   }
00921 
00926   df3_two_variable operator / (const df3_two_variable& x,
00927     const df3_two_variable& y)
00928   {
00929     df3_two_variable u=inv(y);
00930     return x*u;
00931   }
00932 
00937   df3_two_variable operator / (const double x,
00938     const df3_two_variable& y)
00939   {
00940     df3_two_variable u=inv(y);
00941     return x*u;
00942   }
00943 
00948   df3_two_variable operator + (const double x,const df3_two_variable& y)
00949   {
00950     df3_two_variable z;
00951     *z.get_u() =  x + *y.get_u();
00952     *z.get_u_x() = *y.get_u_x();
00953     *z.get_u_y() = *y.get_u_y();
00954     *z.get_u_xx() = *y.get_u_xx();
00955     *z.get_u_xy() = *y.get_u_xy();
00956     *z.get_u_yy() = *y.get_u_yy();
00957     *z.get_u_xxx() = *y.get_u_xxx();
00958     *z.get_u_xxy() = *y.get_u_xxy();
00959     *z.get_u_xyy() = *y.get_u_xyy();
00960     *z.get_u_yyy() = *y.get_u_yyy();
00961     return z;
00962   }
00963 
00968   df3_two_variable operator + (const df3_two_variable& x,const double y)
00969   {
00970     df3_two_variable z;
00971     *z.get_u() =  *x.get_u() + y;
00972     *z.get_u_x() = *x.get_u_x();
00973     *z.get_u_y() = *x.get_u_y();
00974     *z.get_u_xx() = *x.get_u_xx();
00975     *z.get_u_xy() = *x.get_u_xy();
00976     *z.get_u_yy() = *x.get_u_yy();
00977     *z.get_u_xxx() = *x.get_u_xxx();
00978     *z.get_u_xxy() = *x.get_u_xxy();
00979     *z.get_u_xyy() = *x.get_u_xyy();
00980     *z.get_u_yyy() = *x.get_u_yyy();
00981     return z;
00982   }
00983 
00988   df3_two_variable operator + (const df3_two_variable& x,
00989     const df3_two_variable& y)
00990   {
00991     df3_two_variable z;
00992     *z.get_u() = *x.get_u() + *y.get_u();
00993     *z.get_u_x() = *x.get_u_x() + *y.get_u_x();
00994     *z.get_u_y() = *x.get_u_y()+*y.get_u_y();
00995     *z.get_u_xx() = *x.get_u_xx()+ *y.get_u_xx();
00996     *z.get_u_xy() = *x.get_u_xy()+ *y.get_u_xy();
00997     *z.get_u_yy() = *x.get_u_yy()+ *y.get_u_yy();
00998     *z.get_u_xxx() = *x.get_u_xxx()+ *y.get_u_xxx();
00999     *z.get_u_xxy() = *x.get_u_xxy()+ *y.get_u_xxy();
01000     *z.get_u_xyy() = *x.get_u_xyy()+ *y.get_u_xyy();
01001     *z.get_u_yyy() = *x.get_u_yyy()+ *y.get_u_yyy();
01002     return z;
01003   }
01004 
01009   df3_two_variable operator - (const df3_two_variable& x,
01010     const df3_two_variable& y)
01011   {
01012     df3_two_variable z;
01013     *z.get_u() = *x.get_u() - *y.get_u();
01014     *z.get_u_x() = *x.get_u_x()  - *y.get_u_x();
01015     *z.get_u_y() = *x.get_u_y() - *y.get_u_y();
01016     *z.get_u_xx() = *x.get_u_xx() - *y.get_u_xx();
01017     *z.get_u_xy() = *x.get_u_xy() - *y.get_u_xy();
01018     *z.get_u_yy() = *x.get_u_yy() - *y.get_u_yy();
01019     *z.get_u_xxx() = *x.get_u_xxx() - *y.get_u_xxx();
01020     *z.get_u_xxy() = *x.get_u_xxy() - *y.get_u_xxy();
01021     *z.get_u_xyy() = *x.get_u_xyy() - *y.get_u_xyy();
01022     *z.get_u_yyy() = *x.get_u_yyy() - *y.get_u_yyy();
01023     return z;
01024   }
01025 
01030   df3_two_variable operator - (double x,
01031     const df3_two_variable& y)
01032   {
01033     df3_two_variable z;
01034     *z.get_u() = x - *y.get_u();
01035     *z.get_u_x() = - *y.get_u_x();
01036     *z.get_u_y() = - *y.get_u_y();
01037     *z.get_u_xx() = - *y.get_u_xx();
01038     *z.get_u_xy() = - *y.get_u_xy();
01039     *z.get_u_yy() = - *y.get_u_yy();
01040     *z.get_u_xxx() = - *y.get_u_xxx();
01041     *z.get_u_xxy() = - *y.get_u_xxy();
01042     *z.get_u_xyy() = - *y.get_u_xyy();
01043     *z.get_u_yyy() = - *y.get_u_yyy();
01044     return z;
01045   }
01046 
01051   df3_two_variable operator - (const df3_two_variable& x,
01052     double y)
01053   {
01054     df3_two_variable z;
01055     *z.get_u() = *x.get_u()-y;
01056     *z.get_u_x() = *x.get_u_x();
01057     *z.get_u_y() = *x.get_u_y();
01058     *z.get_u_xx() = *x.get_u_xx();
01059     *z.get_u_xy() = *x.get_u_xy();
01060     *z.get_u_yy() = *x.get_u_yy();
01061     *z.get_u_xxx() = *x.get_u_xxx();
01062     *z.get_u_xxy() = *x.get_u_xxy();
01063     *z.get_u_xyy() = *x.get_u_xyy();
01064     *z.get_u_yyy() = *x.get_u_yyy();
01065     return z;
01066   }
01067 
01072   init_df3_two_variable::init_df3_two_variable(const df1b2variable& _v)
01073   {
01074     ADUNCONST(df1b2variable,v)
01075     /*
01076     if (ind_var != 0)
01077     {
01078       cerr << " can only have 1 independent_variable in a reverse funnel"
01079            << endl;
01080       ad_exit(1);
01081     }
01082     */
01083     if (num_ind_var>1)
01084     {
01085       cerr << "can only have 2 independent_variables in df3_two_variable"
01086        " function" << endl;
01087       ad_exit(1);
01088     }
01089     ind_var[num_ind_var++]=&v;
01090     *get_u() =  *v.get_u();
01091     switch(num_ind_var)
01092     {
01093     case 1:
01094       *get_u_x() = 1.0;
01095       *get_u_y() = 0.0;
01096       break;
01097     case 2:
01098       *get_u_x() = 0.0;
01099       *get_u_y() = 1.0;
01100       break;
01101     default:
01102       cerr << "illegal num_ind_var value of " << num_ind_var
01103            << " in  df3_two_variable function" << endl;
01104       ad_exit(1);
01105     }
01106     *get_u_xx() = 0.0;
01107     *get_u_xy() = 0.0;
01108     *get_u_yy() = 0.0;
01109     *get_u_xxx() = 0.0;
01110     *get_u_xxy() = 0.0;
01111     *get_u_xyy() = 0.0;
01112     *get_u_yyy() = 0.0;
01113   }
01114 
01119   init_df3_two_variable::init_df3_two_variable(double v)
01120   {
01121     *get_u() =  v;
01122     *get_u_x() = 0.0;
01123     *get_u_y() = 0.0;
01124     *get_u_xx() = 0.0;
01125     *get_u_xy() = 0.0;
01126     *get_u_yy() = 0.0;
01127     *get_u_xxx() = 0.0;
01128     *get_u_xxy() = 0.0;
01129     *get_u_xyy() = 0.0;
01130     *get_u_yyy() = 0.0;
01131   }
01132 
01133   df3_two_variable::df3_two_variable(void)
01134   {
01135   }
01136 
01141 df3_two_matrix choleski_decomp(const df3_two_matrix& MM)
01142 {
01143   // kludge to deal with constantness
01144   df3_two_matrix & M= (df3_two_matrix &) MM;
01145   int rmin=M.indexmin();
01146   int cmin=M(rmin).indexmin();
01147   int rmax=M.indexmax();
01148   int cmax=M(rmin).indexmax();
01149   if (rmin !=1 || cmin !=1)
01150   {
01151     cerr << "minimum row and column inidices must equal 1 in "
01152       "df1b2matrix choleski_decomp(const df3_two_atrix& MM)"
01153          << endl;
01154     ad_exit(1);
01155   }
01156   if (rmax !=cmax)
01157   {
01158     cerr << "Error in df1b2matrix choleski_decomp(const df3_two_matrix& MM)"
01159       " Matrix not square" << endl;
01160     ad_exit(1);
01161   }
01162 
01163   int n=rmax-rmin+1;
01164   df3_two_matrix L(1,n,1,n);
01165 #ifndef SAFE_INITIALIZE
01166     L.initialize();
01167 #endif
01168 
01169   int i,j,k;
01170   df3_two_variable tmp;
01171 
01172     if (value(M(1,1))<=0)
01173     {
01174       cerr << "Error matrix not positive definite in choleski_decomp"
01175         <<endl;
01176       ad_exit(1);
01177     }
01178 
01179   L(1,1)=sqrt(M(1,1));
01180   for (i=2;i<=n;i++)
01181   {
01182     L(i,1)=M(i,1)/L(1,1);
01183   }
01184 
01185   for (i=2;i<=n;i++)
01186   {
01187     for (j=2;j<=i-1;j++)
01188     {
01189       tmp=M(i,j);
01190       for (k=1;k<=j-1;k++)
01191       {
01192         tmp-=L(i,k)*L(j,k);
01193       }
01194       L(i,j)=tmp/L(j,j);
01195     }
01196     tmp=M(i,i);
01197     for (k=1;k<=i-1;k++)
01198     {
01199       tmp-=L(i,k)*L(i,k);
01200     }
01201 
01202     if (value(tmp)<=0)
01203     {
01204       cerr << "Error matrix not positive definite in choleski_decomp"
01205         <<endl;
01206       ad_exit(1);
01207     }
01208 
01209     L(i,i)=sqrt(tmp);
01210   }
01211 
01212   return L;
01213 }
01214 
01219 df1b2matrix& df1b2matrix::operator = (const df3_two_matrix& M)
01220 {
01221   int rmin=M.indexmin();
01222   int rmax=M.indexmax();
01223   if (rmin != indexmin() || rmax != indexmax())
01224   {
01225     cerr << "unequal shape in "
01226      "df1b2matrix& df1b2matrix::operator = (const df3_two_matrix& M)"
01227       << endl;
01228     ad_exit(1);
01229   }
01230 
01231   for (int i=rmin;i<=rmax;i++)
01232   {
01233     (*this)(i)=M(i);
01234   }
01235   return *this;
01236 }
01237 
01242 df1b2vector& df1b2vector::operator = (const df3_two_vector& M)
01243 {
01244   int rmin=M.indexmin();
01245   int rmax=M.indexmax();
01246   if (rmin != indexmin() || rmax != indexmax())
01247   {
01248     cerr << "unequal shape in "
01249      "df1b2vector& df1b2vector::operator = (const df3_two_vector& M)"
01250       << endl;
01251     ad_exit(1);
01252   }
01253 
01254   for (int i=rmin;i<=rmax;i++)
01255   {
01256     (*this)(i)=M(i);
01257   }
01258   return *this;
01259 }
01260 
01265 df1b2variable& df1b2variable::operator = (const df3_two_variable& v)
01266 {
01267   df1b2variable * px=df3_two_variable::ind_var[0];
01268   df1b2variable * py=df3_two_variable::ind_var[1];
01269   df3_two_variable::num_ind_var=0;
01270   df3_two_variable::ind_var[0]=0;
01271   df3_two_variable::ind_var[1]=0;
01272   //df1b2variable * px=0;
01273   double  dfx= *v.get_u_x();
01274   double  dfy= *v.get_u_y();
01275   double * xd=px->get_u_dot();
01276   double * yd=py->get_u_dot();
01277   double * zd=get_u_dot();
01278   *get_u()=*v.get_u();
01279   for (int i=0;i<df1b2variable::nvar;i++)
01280   {
01281     *zd++ = dfx * *xd++ + dfy * *yd++;
01282   }
01283 
01284   f1b2gradlist->write_pass1(px,py,this,
01285     *v.get_u_x(),
01286     *v.get_u_y(),
01287     *v.get_u_xx(),
01288     *v.get_u_xy(),
01289     *v.get_u_yy(),
01290     *v.get_u_xxx(),
01291     *v.get_u_xxy(),
01292     *v.get_u_xyy(),
01293     *v.get_u_yyy());
01294   return *this;
01295 }
01296 
01301 df1b2variable div(const df1b2variable& x,const df1b2variable& y)
01302 {
01303   df1b2variable z;
01304   double xu=*x.get_u();
01305   double yu=*y.get_u();
01306   double yinv=1.0/yu;
01307   *z.get_u()=xu*yinv;
01308 
01309   double dfx= yinv;
01310   double dfy= -xu*yinv*yinv;
01311   double dfxx= 0.0;
01312   double dfxy=-yinv*yinv;
01313   double dfyy=2.0*xu*yinv*yinv*yinv;
01314   double dfxxx= 0.0;
01315   double dfxxy= 0.0;
01316   double dfxyy=2.0*yinv*yinv*yinv;
01317   double dfyyy=-6.0*xu*yinv*yinv*yinv*yinv;
01318 
01319   double * xd=x.get_u_dot();
01320   double * yd=y.get_u_dot();
01321   double * zd=z.get_u_dot();
01322 
01323   for (int i=0;i<df1b2variable::nvar;i++)
01324   {
01325     *zd++ = dfx * *xd++ + dfy * *yd++;
01326   }
01327 
01328   f1b2gradlist->write_pass1(&x,&y,&z,
01329     dfx,
01330     dfy,
01331     dfxx,dfxy,dfyy,
01332     dfxxx,dfxxy,dfxyy,dfyyy);
01333   return z;
01334 }
01335 
01340 df1b2variable mypow(const df1b2variable& x,double y)
01341 {
01342   df1b2variable z;
01343   double xu=*x.get_u();
01344   *z.get_u()=::pow(xu,y);
01345 
01346   double dfx= y*::pow(xu,y-1.0);
01347   double dfxx= y*(y-1.0)*::pow(xu,y-2.0);
01348   double dfxxx= y*(y-1.0)*(y-2.0)*::pow(xu,y-3.0);
01349   double * xd=x.get_u_dot();
01350   double * zd=z.get_u_dot();
01351 
01352   for (int i=0;i<df1b2variable::nvar;i++)
01353   {
01354     *zd++ = dfx * *xd++ ;
01355   }
01356 
01357   f1b2gradlist->write_pass1(&x,&z,
01358     dfx,
01359     dfxx,
01360     dfxxx);
01361   return z;
01362 }