This file is indexed.

/usr/include/rheolef/newton-backtrack.h is in librheolef-dev 6.5-1build1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

 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
# ifndef _RHEO_NEWTON_BACKTRACK_H
# define _RHEO_NEWTON_BACKTRACK_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef 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 General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
namespace rheolef { 

template <class Problem, class Preconditioner, class Field, class Real>
int newton_backtrack (
  Problem P, Preconditioner T,
  Field u_old, Float Tu_old, Field delta_u, Real slope, Real norm_delta_u_max,
  Field& u, Field& Fu, Real& Tu, Real& lambda)
{
  const Float alpha    = 1e-4; // 1e-8 when strongly nonlinear
  const Float eps_mach = std::numeric_limits<Float>::epsilon();
  Float norm_delta_u = P.space_norm(delta_u);
  if (norm_delta_u > norm_delta_u_max) {
    Float c = norm_delta_u_max/norm_delta_u;
    warning_macro ("delta_u bounded by factor " << c);
    delta_u = c*delta_u;
  }
  // compute lambda_min
  Float norm_u = P.space_norm(u);
  if (norm_u < norm_delta_u) norm_u = norm_delta_u;
  Float lambda_min = eps_mach*norm_u/norm_delta_u;
  if (lambda_min > 1) { // machine precision problem detected
    u = u_old;
    return 1;
  }
  lambda = 1;
  Float lambda_prev = 0;
  Float Tu_prev = 0;
  Float Tu_prev_old = 0;
  for (size_t k = 0; true; k++) {
    u = u_old + lambda*delta_u;
    Fu = P.residue(u);
    Tu = T(P,Fu);
    Float lambda_next;
    if (lambda < lambda_min) { // machine precision problem detected
      u = u_old;
      return 1;
    } else if (Tu <= Tu_old + alpha*lambda*slope) {
      return 0; // have a valid lambda
    } else if (lambda == 1) {
      // first iteration: first order recursion
      lambda_next = - 0.5*slope/(Tu - Tu_old - slope);
    } else {
      // second and more iterations: second order recursion
      Float z      = Tu      - Tu_old      - lambda*slope;
      Float z_prev = Tu_prev - Tu_prev_old - lambda_prev*slope;
      Float a = (  z/sqr(lambda) - z_prev/sqr(lambda_prev))/(lambda - lambda_prev);
      Float b = (- z*lambda_prev/sqr(lambda) + z_prev*lambda/sqr(lambda_prev))
               /(lambda - lambda_prev);
      if (a == 0) {
        lambda_next = - slope/(2*b);
      } else {
        Float Delta = sqr(b) - 3*a*slope;
        if (Delta < 0) {
          warning_macro("newton/backtrack: machine precision problem detected");
          return 1;
        }
        lambda_next = (-b + ::sqrt(Delta))/(3*a);
      }
      lambda_next = std::min (lambda/2, lambda_next);
    }
    lambda_prev = lambda;
    Tu_prev = Tu;
    Tu_prev_old = Tu_old;
    lambda = std::max (lambda/10, lambda_next);
  }
  return 0;
}
}// namespace rheolef
# endif // _RHEO_NEWTON_BACKTRACK_H