Revision 1109 trunk/src/df1b2separable/df1b2lp11.cpp
df1b2lp11.cpp (revision 1109)  

2  2 
* $Id$ 
3  3 
* 
4  4 
* Author: David Fournier 
5 
* Copyright (c) 20082012 Regents of the University of California


5 
* Copyright (c) 20082012 Regents of the University of California 

6  6 
*/ 
7  7 
/** 
8  8 
* \file 
...  ...  
86  86 
int minsave=M.indexmin(); 
87  87 
M.shift(1); 
88  88 
int n=M.indexmax(); 
89 


89  
90  90 
int bw=M.bandwidth(); 
91  91 
banded_lower_triangular_dmatrix L(1,n,bw); 
92  92 
#ifndef SAFE_INITIALIZE 
...  ...  
168  168 
const banded_lower_triangular_dmatrix& C=quiet_choleski_decomp(m,ierr); 
169  169 
if (ierr==0) 
170  170 
{ 
171 
id/=2.0;


171 
id/=2.0; 

172  172 
w=solve_trans(C,::solve(C,v)); 
173  173 
dvector delta=m*w; 
174  174 
dvector err=solve_trans(C,::solve(C,vdelta)); 
...  ...  
183  183 
} 
184  184 
else 
185  185 
{ 
186 
id*=2.0;


186 
id*=2.0; 

187  187 
positivize(m,id); 
188  188 
ierr=0; 
189  189 
dirty=1; 
...  ...  
205  205 
do_newton_raphson_state_space(function_minimizer * pfmin,double f_from_1, 
206  206 
int& no_converge_flag) 
207  207 
{ 
208 
laplace_approximation_calculator::where_are_we_flag=2;


208 
laplace_approximation_calculator::where_are_we_flag=2; 

209  209 
double fbest=1.e+100; 
210  210 
double fval=1.e+100; 
211  211 
double maxg=fabs(evaluate_function(fbest,uhat,pfmin)); 
212  212  
213  213  
214 
laplace_approximation_calculator::where_are_we_flag=0;


214 
laplace_approximation_calculator::where_are_we_flag=0; 

215  215 
dvector uhat_old(1,usize); 
216  216 
safe_choleski_solver scs(0.1); 
217  217 
//for(int ii=1;ii<=num_nr_iters;ii++) 
218  218 
int ii=0; 
219  219 
do 
220 
{


220 
{ 

221  221 
bHess>initialize(); 
222  222  
223  223 
grad.initialize(); 
...  ...  
251  251 
{ 
252  252 
cout << " grad compare " << norm(g1grad) << endl; 
253  253 
} 
254 
step=scs.solve(*bHess,g1);


255 
//step=scs.solve(*bHess,grad);


254 
step=scs.solve(*bHess,g1); 

255 
//step=scs.solve(*bHess,grad); 

256  256 
if (nr_debug==1) 
257  257 
{ 
258  258 
cout << " angle = " << step*grad/(norm(step)*norm(grad)) << endl; 
...  ...  
265  265 
int smallshrink=0; 
266  266 
do 
267  267 
{ 
268 
if (++iic>10)


268 
if (++iic>10) 

269  269 
{ 
270  270 
ibad=1; 
271  271 
break; 
...  ...  
282  282 
maxg=fabs(evaluate_function(fval,utry,g,pfmin)); 
283  283 
if (nr_debug==1) 
284  284 
{ 
285 
cout << " fbestfval = " << setprecision(15)


285 
cout << " fbestfval = " << setprecision(15) 

286  286 
<< fbestfval << endl; 
287  287 
} 
288  288 
if (fval>fbest && maxg>1.e10) 
...  ...  
293  293 
smallshrink=2; 
294  294 
else if (maxg<1.e8) 
295  295 
smallshrink=1; 
296 


296  
297  297 
if (nr_debug==1) 
298  298 
{ 
299  299 
testangle=g*step/(norm(g)*norm(step)); 
Also available in: Unified diff