1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5
6 #include <algorithm>
7 #include <math.h>
8 #include "GmshMessage.h"
9 #include "ConjugateGradients.h"
10 /*
11 min_a f(x+a*d);
12 f(x+a*d) = f(x) + f'(x) (
13 */
14
_norm(std::vector<double> & x)15 static double _norm(std::vector<double> &x)
16 {
17 double n = 0.0;
18 for(std::size_t i = 0; i < x.size(); i++) n += x[i] * x[i];
19 return sqrt(n);
20 }
scale(std::vector<double> & x,double s)21 static void scale(std::vector<double> &x, double s)
22 {
23 for(std::size_t i = 0; i < x.size(); i++) x[i] *= s;
24 }
25
gmshLineSearch(void (* func)(std::vector<double> & x,double & Obj,bool needGrad,std::vector<double> & gradObj,void *),void * data,std::vector<double> & x,std::vector<double> & p,std::vector<double> & g,double & f,double stpmax,int & check)26 static void gmshLineSearch(void (*func)(std::vector<double> &x, double &Obj,
27 bool needGrad,
28 std::vector<double> &gradObj,
29 void *), // the function
30 void *data, // eventual data
31 std::vector<double> &x, // variables
32 std::vector<double> &p, // search direction
33 std::vector<double> &g, // gradient
34 double &f, double stpmax, int &check)
35 {
36 int i;
37 double alam, alam2 = 1., alamin, f2 = 0., fold2 = 0., rhs1, rhs2, temp,
38 tmplam = 0.;
39
40 const double ALF = 1.e-4;
41 const double TOLX = 1.0e-9;
42
43 std::vector<double> xold(x), grad(x);
44 double fold;
45 (*func)(xold, fold, false, grad, data);
46
47 check = 0;
48 int n = x.size();
49 double norm = _norm(p);
50 if(norm > stpmax) scale(p, stpmax / norm);
51 double slope = 0.0;
52 for(i = 0; i < n; i++) slope += g[i] * p[i];
53 double test = 0.0;
54 for(i = 0; i < n; i++) {
55 temp = fabs(p[i]) / std::max(fabs(xold[i]), 1.0);
56 if(temp > test) test = temp;
57 }
58
59 alamin = TOLX / test;
60 alam = 1.;
61
62 while(1) {
63 for(i = 0; i < n; i++) x[i] = xold[i] + alam * p[i];
64 (*func)(x, f, false, grad, data);
65 if(f > 1.e280)
66 alam *= .5;
67 else
68 break;
69 }
70
71 while(1) {
72 for(i = 0; i < n; i++) x[i] = xold[i] + alam * p[i];
73 (*func)(x, f, false, grad, data);
74
75 // printf("alam = %12.5E alamin = %12.5E f = %12.5E fold = %12.5E slope
76 // %12.5E stuff %12.5E\n",alam,alamin,f,fold, slope, ALF * alam *
77 // slope);
78
79 // f = (*func)(x, data);
80 if(alam < alamin) {
81 for(i = 0; i < n; i++) x[i] = xold[i];
82 check = 1;
83 return;
84 }
85 else if(f <= fold + ALF * alam * slope)
86 return;
87 else {
88 if(alam == 1.0)
89 tmplam = -slope / (2.0 * (f - fold - slope));
90 else {
91 rhs1 = f - fold - alam * slope;
92 rhs2 = f2 - fold2 - alam2 * slope;
93 const double a =
94 (rhs1 / (alam * alam) - rhs2 / (alam2 * alam2)) / (alam - alam2);
95 const double b =
96 (-alam2 * rhs1 / (alam * alam) + alam * rhs2 / (alam2 * alam2)) /
97 (alam - alam2);
98 if(a == 0.0)
99 tmplam = -slope / (2.0 * b);
100 else {
101 const double disc = b * b - 3.0 * a * slope;
102 if(disc < 0.0)
103 Msg::Error("Roundoff problem in gmshLineSearch.");
104 else
105 tmplam = (-b + sqrt(disc)) / (3.0 * a);
106 }
107 if(tmplam > 0.5 * alam) tmplam = 0.5 * alam;
108 }
109 }
110 alam2 = alam;
111 f2 = f;
112 fold2 = fold;
113 alam = std::max(tmplam, 0.1 * alam);
114 }
115 }
116
117 // Simple Gradient Descent Minimization (use finite differences to compute the
118 // gradient)
119
GradientDescent(void (* func)(std::vector<double> & x,double & Obj,bool needGrad,std::vector<double> & gradObj,void *),std::vector<double> & x,void * data)120 double GradientDescent(void (*func)(std::vector<double> &x, double &Obj,
121 bool needGrad, std::vector<double> &gradObj,
122 void *), // its gradient
123 std::vector<double> &x, // The variables
124 void *data) // User data
125 {
126 const int MAXIT = 3;
127 const int N = x.size();
128
129 std::vector<double> grad(N);
130 std::vector<double> dir(N);
131 double f;
132
133 // printf("entering gradient descent (%d unknowns)...\n",N);
134
135 for(int iter = 0; iter < MAXIT; iter++) {
136 // compute gradient of func
137 double stpmax = 100000;
138 func(x, f, true, grad, data);
139 // printf("computing f ... %22.15E\n",f);
140 for(int i = 0; i < N; i++) dir[i] = -grad[i];
141 int check;
142 gmshLineSearch(func, data, x, dir, grad, f, stpmax, check);
143 // printf("line search is done, f reduces to %22.15E\n",f);
144 // printf("Line search done x = (%g %g) f = %g\n",x(0),x(1),f);
145 if(check == 1) break;
146 }
147 return f;
148 }
149