ext/numo/liblinear/src/newton.cpp in numo-liblinear-2.0.0 vs ext/numo/liblinear/src/newton.cpp in numo-liblinear-2.1.0

@@ -115,59 +115,60 @@
 	fun_obj->grad(w0, g);
 	double gnorm0 = dnrm2_(&n, g, &inc);
 	delete [] w0;
 
 	f = fun_obj->fun(w);
-	info("init f %5.3e\n", f);
 	fun_obj->grad(w, g);
 	double gnorm = dnrm2_(&n, g, &inc);
+	info("init f %5.3e |g| %5.3e\n", f, gnorm);
 
 	if (gnorm <= eps*gnorm0)
 		search = 0;
 
-	double *w_new = new double[n];
 	while (iter <= max_iter && search)
 	{
 		fun_obj->get_diag_preconditioner(M);
 		for(i=0; i<n; i++)
 			M[i] = (1-alpha_pcg) + alpha_pcg*M[i];
 		cg_iter = pcg(g, M, s, r);
 
 		fold = f;
-		step_size = fun_obj->linesearch_and_update(w, s, & f, g, init_step_size);
+		step_size = fun_obj->linesearch_and_update(w, s, &f, g, init_step_size);
 
 		if (step_size == 0)
 		{
 			info("WARNING: line search fails\n");
 			break;
 		}
 
-		info("iter %2d f %5.3e |g| %5.3e CG %3d step_size %4.2e \n", iter, f, gnorm, cg_iter, step_size);
-
-		actred = fold - f;
-		iter++;
-
 		fun_obj->grad(w, g);
-
 		gnorm = dnrm2_(&n, g, &inc);
+
+		info("iter %2d f %5.3e |g| %5.3e CG %3d step_size %4.2e \n", iter, f, gnorm, cg_iter, step_size);
+		
 		if (gnorm <= eps*gnorm0)
 			break;
 		if (f < -1.0e+32)
 		{
 			info("WARNING: f < -1.0e+32\n");
 			break;
 		}
+		actred = fold - f;
 		if (fabs(actred) <= 1.0e-12*fabs(f))
 		{
 			info("WARNING: actred too small\n");
 			break;
 		}
+
+		iter++;
 	}
 
+	if(iter >= max_iter)
+		info("\nWARNING: reaching max number of Newton iterations\n");
+
 	delete[] g;
 	delete[] r;
-	delete[] w_new;
 	delete[] s;
 	delete[] M;
 }
 
 int NEWTON::pcg(double *g, double *M, double *s, double *r)
@@ -175,11 +176,11 @@
 	int i, inc = 1;
 	int n = fun_obj->get_nr_variable();
 	double one = 1;
 	double *d = new double[n];
 	double *Hd = new double[n];
-	double zTr, znewTrnew, alpha, beta, cgtol;
+	double zTr, znewTrnew, alpha, beta, cgtol, dHd;
 	double *z = new double[n];
 	double Q = 0, newQ, Qdiff;
 
 	for (i=0; i<n; i++)
 	{
@@ -196,13 +197,18 @@
 	int max_cg_iter = max(n, 5);
 
 	while (cg_iter < max_cg_iter)
 	{
 		cg_iter++;
-		fun_obj->Hv(d, Hd);
 
-		alpha = zTr/ddot_(&n, d, &inc, Hd, &inc);
+		fun_obj->Hv(d, Hd);
+		dHd = ddot_(&n, d, &inc, Hd, &inc);
+		// avoid 0/0 in getting alpha
+		if (dHd <= 1.0e-16)
+			break;
+		
+		alpha = zTr/dHd;
 		daxpy_(&n, &alpha, d, &inc, s, &inc);
 		alpha = -alpha;
 		daxpy_(&n, &alpha, Hd, &inc, r, &inc);
 
 		// Using quadratic approximation as CG stopping criterion
@@ -234,10 +240,10 @@
 
 	delete[] d;
 	delete[] Hd;
 	delete[] z;
 
-	return(cg_iter);
+	return cg_iter;
 }
 
 void NEWTON::set_print_string(void (*print_string) (const char *buf))
 {
 	newton_print_string = print_string;