From a415926866ae1a0f563a95e192d88115b9d19926 Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Tue, 12 Jun 2018 12:47:06 +0200
Subject: [PATCH] removal of redundant memory copies, stopping when grad norm
 too low

---
 src/lib/Dirac/lbfgs.c               | 42 ++++++++++++++++++++---------
 src/lib/Dirac/lbfgs_nocuda.c        |  4 +--
 src/lib/Dirac/robust_lbfgs_nocuda.c |  4 +--
 3 files changed, 33 insertions(+), 17 deletions(-)

diff --git a/src/lib/Dirac/lbfgs.c b/src/lib/Dirac/lbfgs.c
index 2486dda..d853fba 100644
--- a/src/lib/Dirac/lbfgs.c
+++ b/src/lib/Dirac/lbfgs.c
@@ -351,8 +351,9 @@ cubic_interp(
 //  f0d=(p01*p01-p02*p02)/(2.0*step);
   f0d=(p01-p02)/(2.0*step);
 
-  my_dcopy(m,xk,1,xp,1); /* xp<=xk */
-  my_daxpy(m,pk,b,xp); /* xp<=xp+(b)*pk */
+  //my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp=xk+(a-step)*pk */
+  //my_daxpy(m,pk,b,xp); /* xp<=xp+(b)*pk */
+  my_daxpy(m,pk,-a+step+b,xp); /* xp<=xp+(b)*pk */
 //  func(xp,x,m,n,adata);
 //  my_daxpy(n,xo,-1.0,x);
 //  f1=my_dnrm2(n,x);
@@ -416,8 +417,9 @@ cubic_interp(
     fz0=f0+f1;
    } else {
     /* evaluate function for this root */
-    my_dcopy(m,xk,1,xp,1); /* xp<=xk */
-    my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(z0)*pk */
+    //my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp=xk+(b-step)*pk */
+    //my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
+    my_daxpy(m,pk,-b+step+a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
 //    func(xp,x,m,n,adata);
 //    my_daxpy(n,xo,-1.0,x);
 //    fz0=my_dnrm2(n,x);
@@ -510,8 +512,9 @@ linesearch_zoom(
 
 
     /* evaluate phi(aj) */
-    my_dcopy(m,xk,1,xp,1); /* xp<=xk */
-    my_daxpy(m,pk,aj,xp); /* xp<=xp+(alphaj)*pk */
+    //my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp = xk+alphaj*pk */
+    //my_daxpy(m,pk,aj,xp); /* xp<=xp+(aj)*pk */
+    my_daxpy(m,pk,-alphaj+aj,xp); /* xp<=xp+(aj)*pk */
 //    func(xp,x,m,n,adata);
     /* calculate x<=x-xo */
 //    my_daxpy(n,xo,-1.0,x);
@@ -530,8 +533,9 @@ linesearch_zoom(
       bj=alphaj; /* aj unchanged */
     } else {
      /* evaluate grad(alphaj) */
-     my_dcopy(m,xk,1,xp,1); /* xp<=xk */
-     my_daxpy(m,pk,alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
+     //my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp = xk+aj*pk */
+     //my_daxpy(m,pk,alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
+     my_daxpy(m,pk,-aj+alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
 //     func(xp,x,m,n,adata);
      /* calculate x<=x-xo */
 //     my_daxpy(n,xo,-1.0,x);
@@ -696,6 +700,16 @@ linesearch(
 #ifdef DEBUG
   printf("cost=%lf grad=%lf mu=%lf, alpha1=%lf\n",phi_0,gphi_0,mu,alpha1);
 #endif
+  /* catch if not finite (deltaphi=0 or nan) */
+  if (!isnormal(mu)) {
+    free(x);
+    free(xp);
+#ifdef DEBUG
+    printf("line interval too small\n");
+#endif
+    return mu;
+  }
+
 
   ci=1;
   alphai=alpha1; /* initial value for alpha(i) : check if 0<alphai<=mu */
@@ -737,8 +751,9 @@ linesearch(
    } 
 
    /* evaluate grad(phi(alpha(i))) */
-   my_dcopy(m,xk,1,xp,1); /* NOT NEEDED here?? xp<=xk */
-   my_daxpy(m,pk,alphai+step,xp); /* xp<=xp+(alphai+step)*pk */
+   //my_dcopy(m,xk,1,xp,1); /* NOT NEEDED here because already xp=xk+alphai*pk */
+   //my_daxpy(m,pk,alphai+step,xp); /* xp<=xp+(alphai+step)*pk */
+   my_daxpy(m,pk,step,xp); /* xp<=xp+(alphai+step)*pk */
 //   func(xp,x,m,n,adata);
    /* calculate x<=x-xo */
 //   my_daxpy(n,xo,-1.0,x);
@@ -998,7 +1013,7 @@ lbfgs_fit_common(
   cm=0;
   ci=0;
  
-  while (ck<itmax) {
+  while (ck<itmax && isnormal(gradnrm) && gradnrm>CLM_STOP_THRESH) {
    /* mult with hessian  pk=-H_k*gk */
    if (ck<M) {
     mult_hessian(m,pk,gk,s,y,rho,ck,ci);
@@ -1011,10 +1026,11 @@ lbfgs_fit_common(
    /* parameters alpha1=10.0,sigma=0.1, rho=0.01, t1=9, t2=0.1, t3=0.5 */
    /* FIXME: update paramters for GPU gradient */
    alphak=linesearch(func,xk,pk,10.0,0.1,0.01,9,0.1,0.5,x,m,n,step,adata,&tp, &tpg);
-   /* check if step size is too small, then stop */
-   if (fabs(alphak)<CLM_STOP_THRESH) {
+   /* check if step size is too small, or nan, then stop */
+   if (!isnormal(alphak) || fabs(alphak)<CLM_EPSILON) {
     break;
    }
+
    /* update parameters xk1=xk+alpha_k *pk */
    my_dcopy(m,xk,1,xk1,1);
    my_daxpy(m,pk,alphak,xk1);
diff --git a/src/lib/Dirac/lbfgs_nocuda.c b/src/lib/Dirac/lbfgs_nocuda.c
index 5eec415..d0e986e 100644
--- a/src/lib/Dirac/lbfgs_nocuda.c
+++ b/src/lib/Dirac/lbfgs_nocuda.c
@@ -460,8 +460,8 @@ cubic_interp(
    } else {
     /* evaluate function for this root */
     //my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp=xk+(b-step)*pk */
-    //my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(z0)*pk */
-    my_daxpy(m,pk,-b+step+a+z0*(b-a),xp); /* xp<=xp+(z0)*pk */
+    //my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
+    my_daxpy(m,pk,-b+step+a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
     func(xp,x,m,n,adata);
     my_daxpy(n,xo,-1.0,x);
     fz0=my_dnrm2(n,x);
diff --git a/src/lib/Dirac/robust_lbfgs_nocuda.c b/src/lib/Dirac/robust_lbfgs_nocuda.c
index e8d2a8e..53e61c7 100644
--- a/src/lib/Dirac/robust_lbfgs_nocuda.c
+++ b/src/lib/Dirac/robust_lbfgs_nocuda.c
@@ -314,8 +314,8 @@ cubic_interp(
    } else {
     /* evaluate function for this root */
     //my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp=xk+(b-step)*pk */
-    //my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(z0)*pk */
-    my_daxpy(m,pk,-b+step+a+z0*(b-a),xp); /* xp<=xp+(z0)*pk */
+    //my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
+    my_daxpy(m,pk,-b+step+a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
     func(xp,x,m,n,adata);
     //my_daxpy(n,xo,-1.0,x);
     //fz0=my_dnrm2(n,x);