Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62334579
dfsane_solver.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, May 12, 12:09
Size
5 KB
Mime Type
text/x-c
Expires
Tue, May 14, 12:09 (2 d)
Engine
blob
Format
Raw Data
Handle
17634342
Attached To
rTAMAAS tamaas
dfsane_solver.cpp
View Options
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "dfsane_solver.hh"
/* -------------------------------------------------------------------------- */
namespace
tamaas
{
/* -------------------------------------------------------------------------- */
DFSANESolver
::
DFSANESolver
(
Residual
&
residual
)
:
EPSolver
(
residual
),
search_direction
(
residual
.
getVector
().
dataSize
(),
residual
.
getVector
().
getNbComponents
()),
previous_residual
(
residual
.
getVector
().
dataSize
(),
residual
.
getVector
().
getNbComponents
()),
current_x
(
residual
.
getVector
().
dataSize
(),
residual
.
getVector
().
getNbComponents
()),
delta_x
(
residual
.
getVector
().
dataSize
(),
residual
.
getVector
().
getNbComponents
()),
delta_residual
(
residual
.
getVector
().
dataSize
(),
residual
.
getVector
().
getNbComponents
())
{}
void
DFSANESolver
::
solve
()
{
auto
&
x_prev
=
*
_x
;
UInt
k
=
0
,
nmax
=
500
;
Real
F_norm
=
0
;
previous_merits
.
clear
();
_residual
.
computeResidual
(
x_prev
);
previous_residual
=
_residual
.
getVector
();
const
auto
F_norm0
=
previous_residual
.
l2norm
();
eta
=
[
F_norm0
](
UInt
k
)
{
return
F_norm0
/
((
1
+
k
)
*
(
1
+
k
));
};
const
std
::
pair
<
Real
,
Real
>
sigma_extrema
{
1e-10
,
1e10
};
Real
sigma
=
1
;
do
{
computeSearchDirection
(
sigma
);
lineSearch
(
eta
(
k
));
// No need for a solution update
// it was done automatically in line search
sigma
=
computeSpectralCoeff
(
sigma_extrema
);
x_prev
=
current_x
;
previous_residual
=
_residual
.
getVector
();
F_norm
=
previous_residual
.
l2norm
();
Logger
().
get
(
LogLevel
::
debug
)
<<
k
<<
' '
<<
F_norm
<<
'\n'
;
}
while
(
F_norm
>
abs_tol
and
k
++
<
nmax
);
if
(
k
>=
nmax
)
TAMAAS_EXCEPTION
(
"DF-SANE did not converge"
);
Logger
().
get
(
LogLevel
::
info
)
<<
"DF-SANE: sucessful convergence ("
<<
k
<<
" iterations, "
<<
abs_tol
<<
")
\n
"
;
_residual
.
computeResidualDisplacement
(
*
_x
);
}
Real
DFSANESolver
::
computeSpectralCoeff
(
const
std
::
pair
<
Real
,
Real
>&
bounds
)
{
delta_x
=
current_x
;
delta_x
-=
*
_x
;
delta_residual
=
_residual
.
getVector
();
delta_residual
-=
previous_residual
;
auto
trial
=
delta_x
.
dot
(
delta_x
)
/
delta_x
.
dot
(
delta_residual
);
if
(
std
::
abs
(
trial
)
>=
bounds
.
first
and
std
::
abs
(
trial
)
<=
bounds
.
second
)
return
trial
;
else
{
auto
F
=
previous_residual
.
l2norm
();
// Rule from Cruz et al. (2006)
if
(
F
>
1
)
return
1
;
else
if
(
F
<
1e-5
)
return
1e5
;
else
return
1
/
F
;
}
}
void
DFSANESolver
::
computeSearchDirection
(
Real
sigma
)
{
_residual
.
computeResidual
(
*
_x
);
search_direction
=
_residual
.
getVector
();
search_direction
*=
-
sigma
;
}
void
DFSANESolver
::
lineSearch
(
Real
eta_k
)
{
const
UInt
M
=
10
;
const
Real
gamma
=
1e-4
;
Real
alphap
=
1
,
alpham
=
1
;
const
std
::
pair
<
Real
,
Real
>
tau_extrema
{
0.1
,
0.5
};
// Merit function (nexp = 2)
const
auto
f
=
[
&
](
Real
alpha
,
Real
sign
)
{
current_x
=
search_direction
;
current_x
*=
alpha
*
sign
;
current_x
+=
*
_x
;
_residual
.
computeResidual
(
current_x
);
return
_residual
.
getVector
().
dot
(
_residual
.
getVector
());
};
// Manage previous merits
const
auto
fk
=
_residual
.
getVector
().
dot
(
_residual
.
getVector
());
previous_merits
.
push_front
(
fk
);
if
(
previous_merits
.
size
()
>
M
)
previous_merits
.
pop_back
();
const
auto
fbar
=
*
std
::
max_element
(
previous_merits
.
begin
(),
previous_merits
.
end
());
do
{
Real
fp
=
f
(
alphap
,
+
1
);
if
(
fp
<
fbar
+
eta_k
-
gamma
*
alphap
*
alphap
*
fk
)
return
;
Real
fm
=
f
(
alpham
,
-
1
);
if
(
fm
<
fbar
+
eta_k
-
gamma
*
alpham
*
alpham
*
fk
)
return
;
alphap
=
computeAlpha
(
alphap
,
fp
,
fk
,
tau_extrema
);
alpham
=
computeAlpha
(
alpham
,
fm
,
fk
,
tau_extrema
);
if
(
std
::
isnan
(
alphap
)
or
std
::
isnan
(
alpham
))
TAMAAS_EXCEPTION
(
"NaN error in line search"
);
}
while
(
true
);
}
Real
DFSANESolver
::
computeAlpha
(
const
Real
alpha
,
Real
f
,
Real
fk
,
const
std
::
pair
<
Real
,
Real
>&
bounds
)
{
const
auto
alphat
=
alpha
*
alpha
*
fk
/
(
f
+
(
2
*
alpha
-
1
)
*
fk
);
if
(
alphat
<
bounds
.
first
*
alpha
)
return
bounds
.
first
*
alpha
;
else
if
(
alphat
>
bounds
.
second
*
alpha
)
return
bounds
.
second
*
alpha
;
else
return
alphat
;
}
}
// namespace tamaas
Event Timeline
Log In to Comment