ADMB Documentation  11.1.1932
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp4.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp4.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 #if defined(USE_LAPLACE)
00012 #  include <admodel.h>
00013 #  include <df1b2fun.h>
00014 #  include <adrndeff.h>
00015 void get_second_ders(int xs,int us,const init_df1b2vector y,dmatrix& Hess,
00016   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
00017   laplace_approximation_calculator* lap);
00018 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
00019   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00020   const dmatrix& _Hessadjoint,function_minimizer * pmin);
00021 
00026 void get_second_ders_master(int xs,int us,const init_df1b2vector _y,
00027   dmatrix& Hess,
00028   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
00029   laplace_approximation_calculator * lpc)
00030 {
00031   // Note that xs is the number of active non random effects
00032   // parameters
00033   // us is the number of random effects parameters
00034   int j;
00035   int i;
00036   ADUNCONST(init_df1b2vector,y)
00037   int num_der_blocks=lpc->num_der_blocks;
00038   int xsize=lpc->xsize;
00039   int usize=lpc->usize;
00040 
00041   for (int ip=1;ip<=num_der_blocks;ip++)
00042   {
00043     df1b2variable::minder=lpc->minder(ip);
00044     df1b2variable::maxder=lpc->maxder(ip);
00045     lpc->set_u_dot(ip);
00046     df1b2_gradlist::set_yes_derivatives();
00047     (*re_objective_function_value::pobjfun)=0;
00048     df1b2variable pen=0.0;
00049     df1b2variable zz=0.0;
00050     initial_df1b2params::reset(y,pen);
00051     if (initial_df1b2params::separable_flag)
00052     {
00053       initial_df1b2params::separable_calculation_type=2;
00054       Hess.initialize();
00055       Dux.initialize();
00056     }
00057 
00058     //cout << "2D" << endl;
00059     pfmin->user_function();
00060 
00061     //pfmin->user_function(y,zz);
00062     (*re_objective_function_value::pobjfun)+=pen;
00063     (*re_objective_function_value::pobjfun)+=zz;
00064 
00065     if (!initial_df1b2params::separable_flag)
00066     {
00067       set_dependent_variable(*re_objective_function_value::pobjfun);
00068       df1b2_gradlist::set_no_derivatives();
00069       df1b2variable::passnumber=1;
00070       df1b2_gradcalc1();
00071 
00072       int mind=y(1).minder;
00073       int jmin=max(mind,xsize+1);
00074       int jmax=min(y(1).maxder,xsize+usize);
00075       for (i=1;i<=usize;i++)
00076         for (j=jmin;j<=jmax;j++)
00077           Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00078 
00079       jmin=max(mind,1);
00080       jmax=min(y(1).maxder,xsize);
00081       for (i=1;i<=usize;i++)
00082         for (j=jmin;j<=jmax;j++)
00083           Dux(i,j)=y(i+xsize).u_bar[j-1];
00084     }
00085     if (ip<num_der_blocks)
00086       f1b2gradlist->reset();
00087   }
00088   if (ad_comm::pvm_manager)
00089   {
00090     if (ad_comm::pvm_manager->timing_flag)
00091       ad_comm::pvm_manager->tm.get_elapsed_time_and_reset();
00092   }
00093 
00094   d3_array M1=get_dmatrix_from_slaves();
00095 
00096   if (ad_comm::pvm_manager)
00097   {
00098     if (ad_comm::pvm_manager->timing_flag)
00099     {
00100       double time=ad_comm::pvm_manager->tm.get_elapsed_time_and_reset();
00101       if (ad_comm::global_logfile)
00102       {
00103         (*ad_comm::global_logfile) << "Time in get dmatrix from slaves "
00104           << time << endl;
00105       }
00106     }
00107   }
00108   int k;
00109   int mmin=M1.indexmin();
00110   int mmax=M1.indexmax();
00111   for (k=mmin;k<=mmax;k++)
00112   {
00113     int imin=M1(k).indexmin();
00114     int imax=M1(k).indexmax();
00115     for (i=imin;i<=imax;i++)
00116     {
00117       int jmin=M1(k,i).indexmin();
00118       int jmax=M1(k,i).indexmax();
00119       dvector & M1_ki=M1(k,i);
00120       for (int j=jmin;j<=jmax;j++)
00121       {
00122         Hess(i,j-xsize)=M1_ki(j);
00123       }
00124     }
00125   }
00126   if (ad_comm::pvm_manager)
00127   {
00128     if (ad_comm::pvm_manager->timing_flag)
00129       ad_comm::pvm_manager->tm.get_elapsed_time_and_reset();
00130   }
00131 
00132   imatrix flags=get_int_from_slaves();
00133   d3_array Dux1=get_dmatrix_from_slaves(flags);
00134   if (ad_comm::pvm_manager)
00135   {
00136     if (ad_comm::pvm_manager->timing_flag)
00137     {
00138       double time=ad_comm::pvm_manager->tm.get_elapsed_time_and_reset();
00139       if (ad_comm::global_logfile)
00140       {
00141         (*ad_comm::global_logfile) << "Time in get dmatrix Dux from slaves "
00142           << time << endl;
00143       }
00144     }
00145   }
00146 
00147   mmin=Dux1.indexmin();
00148   mmax=Dux1.indexmax();
00149   for (k=mmin;k<=mmax;k++)
00150   {
00151     if (allocated(Dux1(k)))
00152     {
00153       int imin=Dux1(k).indexmin();
00154       int imax=Dux1(k).indexmax();
00155       for (i=imin;i<=imax;i++)
00156       {
00157         int jmin=Dux1(k,i).indexmin();
00158         int jmax=Dux1(k,i).indexmax();
00159         dvector & Dux1_ki=Dux1(k,i);
00160         for (int j=jmin;j<=jmax;j++)
00161         {
00162           Dux(i,j)=Dux1_ki(j);
00163         }
00164       }
00165     }
00166   }
00167   //cout << "Hess" << endl;
00168   //cout << norm2(Hess-trans(Hess)) << endl;
00169   //cout << "Dux" << endl;
00170   //cout << setprecision(16) << Dux << endl;
00171 }
00172 
00177 dvector laplace_approximation_calculator::default_calculations_parallel_master
00178   (const dvector& _x,const double& _f,function_minimizer * pfmin)
00179 {
00180   // for use when there is no separability
00181   ADUNCONST(dvector,x)
00182   ADUNCONST(double,f)
00183   int i,j;
00184 
00185   send_dvector_to_slaves(x);
00186   initial_params::set_inactive_only_random_effects();
00187   gradient_structure::set_NO_DERIVATIVES();
00188   initial_params::reset(x);    // get current x values into the model
00189   gradient_structure::set_YES_DERIVATIVES();
00190 
00191   initial_params::set_active_only_random_effects();
00192   //int lmn_flag=0;
00193   if (ad_comm::time_flag)
00194   {
00195     if (ad_comm::ptm1)
00196     {
00197       ad_comm::ptm1->get_elapsed_time_and_reset();
00198     }
00199     if (ad_comm::ptm)
00200     {
00201       ad_comm::ptm->get_elapsed_time_and_reset();
00202     }
00203   }
00204   if (ad_comm::time_flag)
00205   {
00206     if (ad_comm::ptm)
00207     {
00208       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00209       if (ad_comm::global_logfile)
00210       {
00211         (*ad_comm::global_logfile) << endl << " Time pos 0 "
00212           << time << endl;
00213       }
00214     }
00215   }
00216 
00217   if (!inner_lmnflag)
00218   {
00219     if (!ADqd_flag)
00220       uhat=get_uhat_quasi_newton(x,pfmin);
00221     else
00222       uhat=get_uhat_quasi_newton_qd(x,pfmin);
00223   }
00224   else
00225   {
00226     uhat=get_uhat_lm_newton(x,pfmin);
00227   }
00228 
00229   if (ad_comm::time_flag)
00230   {
00231     if (ad_comm::ptm)
00232     {
00233       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00234       if (ad_comm::global_logfile)
00235       {
00236         (*ad_comm::global_logfile) << " Time in inner optimization "
00237           << time << endl;
00238       }
00239     }
00240   }
00241 
00242   for (i=1;i<=xsize;i++)
00243   {
00244     y(i)=x(i);
00245   }
00246   for (i=1;i<=usize;i++)
00247   {
00248     y(i+xsize)=uhat(i);
00249   }
00250 
00251   for(int ii=1;ii<=num_nr_iters;ii++)
00252   {
00253     {
00254       // test newton raphson
00255       Hess.initialize();
00256      cout << "Newton raphson " << ii << endl;
00257       get_newton_raphson_info_master(pfmin);
00258 
00259       if (ad_comm::time_flag)
00260       {
00261         if (ad_comm::ptm)
00262         {
00263           double time=ad_comm::ptm->get_elapsed_time_and_reset();
00264           if (ad_comm::global_logfile)
00265           {
00266             (*ad_comm::global_logfile) << " Time in newton-raphson " <<  ii
00267               << "  " << time << endl;
00268           }
00269         }
00270       }
00271         dvector step;
00272 #if defined(USE_ATLAS)
00273       if (!ad_comm::no_atlas_flag)
00274         step=-atlas_solve_spd(Hess,grad);
00275       else
00276         step=-solve(Hess,grad);
00277 #else
00278       step=-solve(Hess,grad);
00279 #endif
00280       if (ad_comm::time_flag)
00281       {
00282         if (ad_comm::ptm)
00283         {
00284           double time=ad_comm::ptm->get_elapsed_time_and_reset();
00285           if (ad_comm::global_logfile)
00286           {
00287             (*ad_comm::global_logfile) << " time_in solve " <<  ii << "  "
00288               << time << endl;
00289           }
00290         }
00291       }
00292       send_dvector_to_slaves(step);
00293       //cout << norm(step1-step) << endl;
00294 
00295       if (ad_comm::time_flag)
00296       {
00297         if (ad_comm::ptm)
00298         {
00299           double time=ad_comm::ptm->get_elapsed_time_and_reset();
00300           if (ad_comm::global_logfile)
00301           {
00302             (*ad_comm::global_logfile) << " Time in send dvector to slaves "
00303               <<  ii <<  "   " << time << endl;
00304           }
00305         }
00306       }
00307 
00308       f1b2gradlist->reset();
00309       f1b2gradlist->list.initialize();
00310       f1b2gradlist->list2.initialize();
00311       f1b2gradlist->list3.initialize();
00312       f1b2gradlist->nlist.initialize();
00313       f1b2gradlist->nlist2.initialize();
00314       f1b2gradlist->nlist3.initialize();
00315 
00316       uhat+=step;
00317 
00318       evaluate_function(uhat,pfmin);
00319     }
00320 
00321     for (i=1;i<=usize;i++)
00322     {
00323       y(i+xsize)=uhat(i);
00324     }
00325   }
00326 
00327   if (ad_comm::time_flag)
00328   {
00329     if (ad_comm::ptm)
00330     {
00331       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00332       if (ad_comm::global_logfile)
00333       {
00334         (*ad_comm::global_logfile) << " Time in reset and evaluate function"
00335           << time << endl;
00336       }
00337     }
00338   }
00339   get_second_ders_master(xsize,usize,y,Hess,Dux,f1b2gradlist,pfmin,this);
00340   //int sgn=0;
00341 
00342   if (ad_comm::time_flag)
00343   {
00344     if (ad_comm::ptm)
00345     {
00346       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00347       if (ad_comm::global_logfile)
00348       {
00349         (*ad_comm::global_logfile) << " Time in dget second ders "
00350           << time << endl;
00351       }
00352     }
00353   }
00354 
00355   f=calculate_laplace_approximation(x,uhat,Hess,xadjoint,uadjoint,
00356     Hessadjoint,pfmin);
00357 
00358   if (ad_comm::time_flag)
00359   {
00360     if (ad_comm::ptm)
00361     {
00362       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00363       if (ad_comm::global_logfile)
00364       {
00365         (*ad_comm::global_logfile) <<
00366           " Time in calculate laplace approximation "
00367           << time << endl;
00368       }
00369     }
00370   }
00371 
00372 
00373   for (int ip=num_der_blocks;ip>=1;ip--)
00374   {
00375     df1b2variable::minder=minder(ip);
00376     df1b2variable::maxder=maxder(ip);
00377     int mind=y(1).minder;
00378     int jmin=max(mind,xsize+1);
00379     int jmax=min(y(1).maxder,xsize+usize);
00380     for (i=1;i<=usize;i++)
00381     {
00382       for (j=jmin;j<=jmax;j++)
00383       {
00384         //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00385         y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00386       }
00387     }
00388 
00389     int mmin=minder.indexmin();
00390     int mmax=minder.indexmax();
00391     ivector jjmin(mmin,mmax-1);
00392     ivector jjmax(mmin,mmax-1);
00393     for (i=mmin+1;i<=mmax;i++)
00394     {
00395       jjmin(i-1)=max(minder(i),xsize+1)-xsize;
00396       jjmax(i-1)=min(maxder(i),xsize+usize)-xsize;
00397     }
00398 
00399     if (ad_comm::pvm_manager)
00400     {
00401       if (ad_comm::pvm_manager->timing_flag)
00402         ad_comm::pvm_manager->tm.get_elapsed_time_and_reset();
00403     }
00404     send_dmatrix_to_slaves(Hessadjoint,jjmin,jjmax);
00405     if (ad_comm::pvm_manager)
00406     {
00407       if (ad_comm::pvm_manager->timing_flag)
00408       {
00409         double time=ad_comm::pvm_manager->tm.get_elapsed_time_and_reset();
00410         if (ad_comm::global_logfile)
00411         {
00412           (*ad_comm::global_logfile) <<
00413             "Time in send dmatrix to slaves Hessadjoint "
00414             << time << endl;
00415         }
00416       }
00417     }
00418 
00419     if (initial_df1b2params::separable_flag)
00420     {
00421       for (j=1;j<=xsize+usize;j++)
00422       {
00423         *y(j).get_u_tilde()=0;
00424       }
00425       Hess.initialize();
00426       initial_df1b2params::separable_calculation_type=3;
00427       pfmin->user_function();
00428     }
00429     else
00430     {
00431       if (ip<num_der_blocks)
00432       {
00433         f1b2gradlist->reset();
00434         set_u_dot(ip);
00435         df1b2_gradlist::set_yes_derivatives();
00436         (*re_objective_function_value::pobjfun)=0;
00437         df1b2variable pen=0.0;
00438         df1b2variable zz=0.0;
00439 
00440         initial_df1b2params::reset(y,pen);
00441         pfmin->user_function();
00442         (*re_objective_function_value::pobjfun)+=pen;
00443         (*re_objective_function_value::pobjfun)+=zz;
00444 
00445         set_dependent_variable(*re_objective_function_value::pobjfun);
00446         df1b2_gradlist::set_no_derivatives();
00447         df1b2variable::passnumber=1;
00448         df1b2_gradcalc1();
00449       }
00450 
00451       for (i=1;i<=usize;i++)
00452       {
00453         for (j=jmin;j<=jmax;j++)
00454         {
00455           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00456           y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00457         }
00458       }
00459 
00460       //int mind=y(1).minder;
00461       df1b2variable::passnumber=2;
00462       df1b2_gradcalc1();
00463 
00464       df1b2variable::passnumber=3;
00465       df1b2_gradcalc1();
00466 
00467       f1b2gradlist->reset();
00468       f1b2gradlist->list.initialize();
00469       f1b2gradlist->list2.initialize();
00470       f1b2gradlist->list3.initialize();
00471       f1b2gradlist->nlist.initialize();
00472       f1b2gradlist->nlist2.initialize();
00473       f1b2gradlist->nlist3.initialize();
00474     }
00475 
00476     if (ad_comm::time_flag)
00477     {
00478       if (ad_comm::ptm)
00479       {
00480         double time=ad_comm::ptm->get_elapsed_time_and_reset();
00481         if (ad_comm::global_logfile)
00482         {
00483           (*ad_comm::global_logfile) << " time for 3rd derivatives "
00484             << time << endl;
00485         }
00486       }
00487     }
00488     dvector dtmp(1,xsize);
00489     for (i=1;i<=xsize;i++)
00490     {
00491       dtmp(i)=*y(i).get_u_tilde();
00492     }
00493     if (initial_df1b2params::separable_flag)
00494     {
00495       dvector scale(1,nvar);   // need to get scale from somewhere
00496       /*int check=*/initial_params::stddev_scale(scale,x);
00497       dvector sscale=scale(1,Dux(1).indexmax());
00498       for (i=1;i<=usize;i++)
00499       {
00500         Dux(i)=elem_prod(Dux(i),sscale);
00501       }
00502       dtmp=elem_prod(dtmp,sscale);
00503     }
00504     if (ad_comm::time_flag)
00505     {
00506       if (ad_comm::ptm)
00507       {
00508         ad_comm::ptm->get_elapsed_time_and_reset();
00509       }
00510     }
00511 
00512     dmatrix slave_dtmps=get_dvector_from_slaves();
00513     if (ad_comm::time_flag)
00514     {
00515       if (ad_comm::ptm)
00516       {
00517         double time=ad_comm::ptm->get_elapsed_time_and_reset();
00518         if (ad_comm::global_logfile)
00519         {
00520           (*ad_comm::global_logfile) <<
00521             " Time in get dvector from slaves slave_dtmp "
00522             << time << endl;
00523         }
00524       }
00525     }
00526     mmin=slave_dtmps.indexmin();
00527     mmax=slave_dtmps.indexmax();
00528 
00529     for (j=mmin;j<=mmax;j++)
00530     {
00531       dvector& s_j=slave_dtmps(j);
00532       for (i=1;i<=xsize;i++)
00533       {
00534         //xadjoint(i)+=s_j(i);
00535         dtmp(i)+=s_j(i);
00536       }
00537     }
00538 
00539     for (i=1;i<=xsize;i++)
00540     {
00541       xadjoint(i)+=dtmp(i);
00542     }
00543     dmatrix slave_utmps=get_dvector_from_slaves();
00544     if (ad_comm::ptm)
00545     {
00546       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00547       if (ad_comm::global_logfile)
00548       {
00549         (*ad_comm::global_logfile) <<
00550           " Time in get dvector from slaves slave_utmps "
00551           << time << endl;
00552       }
00553     }
00554 
00555     mmin=slave_utmps.indexmin();
00556     mmax=slave_utmps.indexmax();
00557 
00558     for (j=mmin;j<=mmax;j++)
00559     {
00560       dvector& s_j=slave_utmps(j);
00561       for (i=1;i<=usize;i++)
00562       {
00563         //xadjoint(i)+=s_j(i);
00564         uadjoint(i)+=s_j(i);
00565       }
00566     }
00567 
00568     for (i=1;i<=usize;i++)
00569       uadjoint(i)+=*y(xsize+i).get_u_tilde();
00570   }
00571 #if defined(USE_ATLAS)
00572   if (!ad_comm::no_atlas_flag)
00573   {
00574     //xadjoint -= uadjoint*atlas_solve_spd_trans(Hess,Dux);
00575     xadjoint -= atlas_solve_spd_trans(Hess,uadjoint)*Dux;
00576   }
00577   else
00578   {
00579     //xadjoint -= uadjoint*inv(Hess)*Dux;
00580     xadjoint -= solve(Hess,uadjoint)*Dux;
00581   }
00582 #else
00583   //xadjoint -= uadjoint*inv(Hess)*Dux;
00584   xadjoint -= solve(Hess,uadjoint)*Dux;
00585 #endif
00586 
00587   if (ad_comm::ptm)
00588   {
00589     double time=ad_comm::ptm->get_elapsed_time_and_reset();
00590     if (ad_comm::global_logfile)
00591     {
00592       (*ad_comm::global_logfile) << " Time in second solve "
00593         << time << endl;
00594     }
00595   }
00596   if (ad_comm::ptm1)
00597   {
00598     double time=ad_comm::ptm1->get_elapsed_time_and_reset();
00599     if (ad_comm::global_logfile)
00600     {
00601       (*ad_comm::global_logfile) << " Total time in function evaluation "
00602         << time << endl << endl;
00603     }
00604   }
00605 
00606   return xadjoint;
00607 }
00608 
00613 void laplace_approximation_calculator::get_newton_raphson_info_master
00614   (function_minimizer * pfmin)
00615 {
00616   int i,j,ip;
00617 
00618   if (ad_comm::time_flag)
00619   {
00620     if (ad_comm::ptm)
00621     {
00622         (*ad_comm::global_logfile) << " Starting Newton-Raphson "
00623           <<  endl;
00624     }
00625   }
00626   for (ip=1;ip<=num_der_blocks;ip++)
00627   {
00628     df1b2variable::minder=minder(ip);
00629     df1b2variable::maxder=maxder(ip);
00630 
00631     set_u_dot(ip);
00632 
00633     // do we need to reallocate memory for df1b2variables?
00634     check_for_need_to_reallocate(ip);
00635     df1b2_gradlist::set_yes_derivatives();
00636     //cout << re_objective_function_value::pobjfun << endl;
00637     //cout << re_objective_function_value::pobjfun->ptr << endl;
00638     (*re_objective_function_value::pobjfun)=0;
00639     df1b2variable pen=0.0;
00640     df1b2variable zz=0.0;
00641     initial_df1b2params::reset(y,pen);
00642     if (initial_df1b2params::separable_flag)
00643     {
00644       Hess.initialize();
00645       grad.initialize();
00646     }
00647 
00648     double time1;
00649     if (ad_comm::time_flag)
00650     {
00651       if (ad_comm::ptm)
00652       {
00653         /*double time1=*/ad_comm::ptm->get_elapsed_time();
00654       }
00655     }
00656 
00657     pfmin->user_function();
00658 
00659     if (ad_comm::time_flag)
00660     {
00661       if (ad_comm::ptm)
00662       {
00663         double time=ad_comm::ptm->get_elapsed_time();
00664         if (ad_comm::global_logfile)
00665         {
00666           (*ad_comm::global_logfile) << "       Time in user_function() "
00667             <<  ip << "  " << time-time1 << endl;
00668         }
00669       }
00670     }
00671 
00672 
00673     (*re_objective_function_value::pobjfun)+=pen;
00674     (*re_objective_function_value::pobjfun)+=zz;
00675 
00676     if (!initial_df1b2params::separable_flag)
00677     {
00678       set_dependent_variable(*re_objective_function_value::pobjfun);
00679       df1b2_gradlist::set_no_derivatives();
00680       df1b2variable::passnumber=1;
00681       df1b2_gradcalc1();
00682       int mind=y(1).minder;
00683       int jmin=max(mind,xsize+1);
00684       int jmax=min(y(1).maxder,xsize+usize);
00685       for (i=1;i<=usize;i++)
00686         for (j=jmin;j<=jmax;j++)
00687           Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00688 
00689       for (j=jmin;j<=jmax;j++)
00690         grad(j-xsize)= re_objective_function_value::pobjfun->u_dot[j-mind];
00691     }
00692     if (ip<num_der_blocks)
00693       f1b2gradlist->reset();
00694   }
00695   double time1;
00696   if (ad_comm::time_flag)
00697   {
00698     if (ad_comm::ptm)
00699     {
00700       time1=ad_comm::ptm->get_elapsed_time();
00701     }
00702   }
00703 
00704   d3_array M1=get_dmatrix_from_slaves();
00705 
00706   if (ad_comm::time_flag)
00707   {
00708     if (ad_comm::ptm)
00709     {
00710       double time=ad_comm::ptm->get_elapsed_time();
00711       if (ad_comm::global_logfile)
00712       {
00713         (*ad_comm::global_logfile) <<
00714           "       Time to getdmatrix from slaves in "
00715           " get newton raphson " <<  ip << "  " << time-time1 << endl;
00716       }
00717     }
00718   }
00719   int k;
00720   int mmin=M1.indexmin();
00721   int mmax=M1.indexmax();
00722   for (k=mmin;k<=mmax;k++)
00723   {
00724     int imin=M1(k).indexmin();
00725     int imax=M1(k).indexmax();
00726     for (i=imin;i<=imax;i++)
00727     {
00728       int jmin=M1(k,i).indexmin();
00729       int jmax=M1(k,i).indexmax();
00730       dvector & M1_ki=M1(k,i);
00731       for (int j=jmin;j<=jmax;j++)
00732       {
00733         Hess(i,j-xsize)=M1_ki(j);
00734       }
00735     }
00736   }
00737 
00738   if (ad_comm::time_flag)
00739   {
00740     if (ad_comm::ptm)
00741     {
00742       /*double time1=*/ad_comm::ptm->get_elapsed_time();
00743     }
00744   }
00745   dmatrix g1=get_dvector_from_slaves();
00746 
00747   if (ad_comm::time_flag)
00748   {
00749     if (ad_comm::ptm)
00750     {
00751       double time=ad_comm::ptm->get_elapsed_time();
00752       if (ad_comm::global_logfile)
00753       {
00754         (*ad_comm::global_logfile) << "       Time in get dvector "
00755             " from slaves in get newton raphson " <<  ip << "  "
00756           << time-time1 << endl;
00757       }
00758     }
00759   }
00760   mmin=g1.indexmin();
00761   mmax=g1.indexmax();
00762   for (k=mmin;k<=mmax;k++)
00763   {
00764     int jmin=g1(k).indexmin();
00765     int jmax=g1(k).indexmax();
00766     dvector & g1_k=g1(k);
00767     for (int j=jmin;j<=jmax;j++)
00768     {
00769       grad(j-xsize)=g1_k(j);
00770     }
00771   }
00772 }
00773 #endif