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