ADMB Documentation  11.1.2364
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp7.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp7.cpp 1935 2014-04-26 02:02:58Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 
00017 void laplace_approximation_calculator::
00018    do_separable_stuff_laplace_approximation_banded_adjoint
00019    (const df1b2variable& ff)
00020 {
00021   //***********************************************************
00022   //***********************************************************
00023   set_dependent_variable(ff);
00024   df1b2_gradlist::set_no_derivatives();
00025   df1b2variable::passnumber=1;
00026   df1b2_gradcalc1();
00027 
00028   init_df1b2vector & locy= *funnel_init_var::py;
00029   imatrix& list=*funnel_init_var::plist;
00030   int num_local_re=0;
00031   int num_local_fe=0;
00032 
00033   int i;
00034   ivector lre_index(1,funnel_init_var::num_active_parameters);
00035   ivector lfe_index(1,funnel_init_var::num_active_parameters);
00036 
00037   for (i=1;i<=funnel_init_var::num_active_parameters;i++)
00038   {
00039     if (list(i,1)>xsize)
00040     {
00041       lre_index(++num_local_re)=i;
00042     }
00043     else if (list(i,1)>0)
00044     {
00045       lfe_index(++num_local_fe)=i;
00046     }
00047   }
00048 
00049   if (num_local_re > 0)
00050   {
00051     int j;
00052     if (hesstype==3)
00053     {
00054       for (i=1;i<=num_local_re;i++)
00055       {
00056         int lrei=lre_index(i);
00057         for (j=1;j<=num_local_re;j++)
00058         {
00059           int lrej=lre_index(j);
00060           int i1=list(lrei,1)-xsize;
00061           int i2=list(lrei,2);
00062           int j1=list(lrej,1)-xsize;
00063           int j2=list(lrej,2);
00064           if (i1>=j1) (*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00065         }
00066       }
00067     }
00068     else if (hesstype==4)
00069     {
00070       if (sparse_hessian_flag==0)
00071       {
00072         for (i=1;i<=num_local_re;i++)
00073         {
00074           int lrei=lre_index(i);
00075           for (j=1;j<=num_local_re;j++)
00076           {
00077             int lrej=lre_index(j);
00078             int i1=list(lrei,1)-xsize;
00079             int i2=list(lrei,2);
00080             int j1=list(lrej,1)-xsize;
00081             int j2=list(lrej,2);
00082             Hess(i1,j1)+=locy(i2).u_bar[j2-1];
00083           }
00084         }
00085       }
00086       else
00087       {
00088         for (i=1;i<=num_local_re;i++)
00089         {
00090           int lrei=lre_index(i);
00091           for (j=1;j<=num_local_re;j++)
00092           {
00093             int lrej=lre_index(j);
00094             int i1=list(lrei,1)-xsize;
00095             int i2=list(lrei,2);
00096             int j1=list(lrej,1)-xsize;
00097             int j2=list(lrej,2);
00098 
00099             if (i1<=j1)
00100             {
00101               sparse_count++;
00102               (*sparse_triplet2)((*sparse_iterator)(sparse_count))
00103                 +=locy(i2).u_bar[j2-1];
00104             }
00105           }
00106         }
00107       }
00108     }
00109   }
00110 
00111   // Now the adjoint code
00112 
00113   if (num_local_re > 0)
00114   {
00115     int j;
00116     if (hesstype==3)
00117     {
00118       for (i=1;i<=num_local_re;i++)
00119       {
00120         int lrei=lre_index(i);
00121         for (j=1;j<=num_local_re;j++)
00122         {
00123           int lrej=lre_index(j);
00124           int i1=list(lrei,1)-xsize;
00125           int i2=list(lrei,2);
00126           int j1=list(lrej,1)-xsize;
00127           int j2=list(lrej,2);
00128           if (i1>=j1)
00129           {
00130             //(*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00131             locy(i2).get_u_bar_tilde()[j2-1]=(*bHessadjoint)(i1,j1);
00132           }
00133         }
00134       }
00135     }
00136     else if (hesstype==4)
00137     {
00138       if (pmin->lapprox->sparse_hessian_flag==0)
00139       {
00140         for (i=1;i<=num_local_re;i++)
00141         {
00142           int lrei=lre_index(i);
00143           for (j=1;j<=num_local_re;j++)
00144           {
00145             int lrej=lre_index(j);
00146             int i1=list(lrei,1)-xsize;
00147             int i2=list(lrei,2);
00148             int j1=list(lrej,1)-xsize;
00149             int j2=list(lrej,2);
00150             {
00151               //(*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00152               locy(i2).get_u_bar_tilde()[j2-1]=Hessadjoint(i1,j1);
00153             }
00154           }
00155         }
00156       }
00157       else
00158       {
00159         dcompressed_triplet * vsparse_triplet_adjoint =
00160           pmin->lapprox->vsparse_triplet_adjoint;
00161         for (i=1;i<=num_local_re;i++)
00162         {
00163           int lrei=lre_index(i);
00164           for (j=1;j<=num_local_re;j++)
00165           {
00166             int lrej=lre_index(j);
00167             int i1=list(lrei,1)-xsize;
00168             int i2=list(lrei,2);
00169             int j1=list(lrej,1)-xsize;
00170             int j2=list(lrej,2);
00171             {
00172               //(*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00173               //locy(i2).get_u_bar_tilde()[j2-1]=Hessadjoint(i1,j1);
00174               if (i1<=j1)
00175               {
00176                 //(*sparse_triplet2)((*sparse_iterator)(sparse_count_adjoint))
00177                 //  +=locy(i2).u_bar[j2-1];
00178                 sparse_count_adjoint++;
00179                 locy(i2).get_u_bar_tilde()[j2-1]=
00180                   (*vsparse_triplet_adjoint)
00181                     ((*sparse_iterator)(sparse_count_adjoint));
00182               }
00183             }
00184           }
00185         }
00186       }
00187     }
00188   }
00189 
00190   df1b2variable::passnumber=2;
00191   df1b2_gradcalc1();
00192 
00193   df1b2variable::passnumber=3;
00194   df1b2_gradcalc1();
00195 
00196   f1b2gradlist->reset();
00197   f1b2gradlist->list.initialize();
00198   f1b2gradlist->list2.initialize();
00199   f1b2gradlist->list3.initialize();
00200   f1b2gradlist->nlist.initialize();
00201   f1b2gradlist->nlist2.initialize();
00202   f1b2gradlist->nlist3.initialize();
00203 
00204 
00205   for (i=1;i<=num_local_fe;i++)
00206   {
00207     int lfei=lfe_index(i);
00208     int i1=list(lfei,1);
00209     //if (lfei - list(lfei,2))
00210      // cout << lfei << " " <<  list(lfei,2) << endl;
00211     local_dtemp(i1)+=*locy(lfei).get_u_tilde();
00212 #if !defined(OPT_LIB)
00213     int mmin=xadjoint.indexmin();
00214     int mmax=xadjoint.indexmax();
00215     if (i1<mmin || i1> mmax)
00216     {
00217        cerr << "this can't happen" << endl;
00218       ad_exit(1);
00219     }
00220 #endif
00221   }
00222 
00223   for (i=1;i<=num_local_re;i++)
00224   {
00225     int i1=list(lre_index(i),1)-xsize;
00226     int i2=list(lre_index(i),2);
00227     uadjoint(i1)+=*locy(i2).get_u_tilde();
00228   }
00229 
00230   f1b2gradlist->reset();
00231   f1b2gradlist->list.initialize();
00232   f1b2gradlist->list2.initialize();
00233   f1b2gradlist->list3.initialize();
00234   f1b2gradlist->nlist.initialize();
00235   f1b2gradlist->nlist2.initialize();
00236   f1b2gradlist->nlist3.initialize();
00237   funnel_init_var::num_vars=0;
00238   funnel_init_var::num_active_parameters=0;
00239   funnel_init_var::num_inactive_vars=0;
00240 }