From eca3bd047f9a53e4f92f51ab753edf195a937e85 Mon Sep 17 00:00:00 2001
From: b92paul <b92paul@gmail.com>
Date: Thu, 8 Mar 2018 15:02:51 +0800
Subject: [PATCH] Fixed bug and optimize code in TRON::tron

Variable "iter" and "w" won't change in the
while loop of function TRON::tron when the
direction does not lead to sufficient
function-value decrease, so using expression
"iter == 1" as a flag to initialize variables
will cause problems in the first step.
We introduce a new variable "delta_adjusted" to
fix the bug.

The function call "fun_obj->get_diagH(M)" and the
calculation of "(1-alpha_pcg) I + alpha_pcg M"
can be saved if "iter" and "w" remains the same.
We move some lines of code to save the
calculation.
---
 tron.cpp | 19 ++++++++++++-------
 1 file changed, 12 insertions(+), 7 deletions(-)

diff --git a/tron.cpp b/tron.cpp
index d662d99..2d93dc3 100644
--- a/tron.cpp
+++ b/tron.cpp
@@ -103,17 +103,16 @@ void TRON::tron(double *w)
 	if (gnorm <= eps*gnorm0)
 		search = 0;
 
-	iter = 1;
+	fun_obj->get_diagH(M);
+	for(i=0; i<n; i++)
+		M[i] = (1-alpha_pcg) + alpha_pcg*M[i];
+	delta = sqrt(uTMv(n, g, M, g));
 
 	double *w_new = new double[n];
 	bool reach_boundary;
+	bool delta_adjusted = false;
 	while (iter <= max_iter && search)
 	{
-		fun_obj->get_diagH(M);
-		for(i=0; i<n; i++)
-			M[i] = (1-alpha_pcg) + alpha_pcg*M[i];
-		if (iter == 1)
-			delta = sqrt(uTMv(n, g, M, g));
 		cg_iter = trpcg(delta, g, M, s, r, &reach_boundary);
 
 		memcpy(w_new, w, sizeof(double)*n);
@@ -128,8 +127,11 @@ void TRON::tron(double *w)
 
 		// On the first iteration, adjust the initial step bound.
 		sMnorm = sqrt(uTMv(n, s, M, s));
-		if (iter == 1)
+		if (iter == 1 && !delta_adjusted)
+		{
 			delta = min(delta, sMnorm);
+			delta_adjusted = true;
+		}
 
 		// Compute prediction alpha*sMnorm of the step.
 		if (fnew - f - gs <= 0)
@@ -160,6 +162,9 @@ void TRON::tron(double *w)
 			memcpy(w, w_new, sizeof(double)*n);
 			f = fnew;
 			fun_obj->grad(w, g);
+			fun_obj->get_diagH(M);
+			for(i=0; i<n; i++)
+				M[i] = (1-alpha_pcg) + alpha_pcg*M[i];
 
 			gnorm = dnrm2_(&n, g, &inc);
 			if (gnorm <= eps*gnorm0)
-- 
2.40.0