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