Newton.cpp

#include "Newton.hpp"


void Newton::operator()(const IAssemble &A, Vector &u) {
  mout.StartBlock("Newton");
  double E, E_0;
  A.Initialize(u);
  vout(9) << "u(0)= " << endl << u << endl;
  Vector r(u);
  Matrix J(u);
  d = d_0 = A.Residual(u, r);
  E = E_0 = A.Energy(u);
  double eps = defaultEpsilon + defaultReduction * d;
  double d_previous = d;
  int LS_cnt = 0;

  int JU_cnt = 0; // counts iteration steps without Jacobian update
  Vector c(r);

  for (iter = 0; iter < max_iter; ++iter) {
    vout(3) << "r(" << iter << ")= " << endl << r << endl;
    if (d < eps) { break; }
    vout(1) << "d(" << iter << ")= " << d << endl;

    // Determines whether Jacobi Matrix is updated
    if (((iter - JU_cnt) >= JacobiUpdate) || (iter == 0)) {
      A.Jacobi(u, J);
      JU_cnt = iter;
      c = (*solver)(J) * r;
    } else c = (*solver)(J, 0) * r;

    vout(3) << "c(" << iter << ")= " << endl << c << endl;
    u -= c;
    vout(5) << "u-c " << endl << u << endl;
    d = A.Residual(u, r);
    E = A.Energy(u);
    if (d > d_previous) {
      for (int l = 1; l <= LS_iter; ++l) {
        if (iter == 0)
          if (suppressLS) {
            vout(2) << "line search suppressed" << endl;
            break;
          }
        vout(1) << "line search " << l << ": d(" << iter << ")= " << d << endl;
        c *= 0.5;
        u += c;
        d = A.Residual(u, r);
        E = A.Energy(u);
        if (d < d_previous) break;
      }
    }
    if (d > d_previous) {
      vout(5) << "line search unsuccessful." << endl;
      ++LS_cnt;
      if (LS_cnt == 3) {
        vout(1) << "too many line searches unsuccessful." << endl;
        iter = max_iter;
      }
    }
    d_previous = d;
  }
  endOut(1);
}


void ENewton::operator()(const IAssemble &A, Vector &u) {
  mout.StartBlock("ENewton");
  double E, E_0;
  A.Initialize(u);
  Vector r(u);
  Matrix J(u);
  d = d_0 = A.Residual(u, r);
  E = E_0 = A.Energy(u);
  double eps = defaultEpsilon + defaultReduction * d;
  int LS_cnt = 0;
  int JU_cnt = 0;
  Vector c(r);
  for (iter = 0; iter < max_iter; ++iter) {
    vout(3) << "r(" << iter << ")= " << r << endl;
    if (d < eps) { break; }
    if (d > VeryLarge) { break; }
    if (!(E < VeryLarge)) {
      if (!(E > VeryLarge)) {
        iter = max_iter;
        break;
      }
    }
    vout(1) << "d(" << iter << ")= " << d << "\t E(" << iter << ")= " << E << endl;
    A.Jacobi(u, J);
    JU_cnt = iter;
    c = (*solver)(J) * r;
//    double rc = std::real(r * c);
    double rc = r * c;
    vout(3) << "c(" << iter << ")= " << c << endl;
    u -= c;
    vout(5) << "u-c " << u << endl;
    double d_old = d;
    d = A.Residual(u, r);
    double E_old = E;
    E = A.Energy(u);
    if (d < d_old) continue;
    if (E > E_old + sigma * rc) {
      for (int l = 1; l <= LS_iter; ++l) {
        if (iter == 0)
          if (suppressLS) {
            vout(2) << "line search suppressed" << endl;
            break;
          }
        vout(1) << "line search " << l
                << ": d(" << iter << ")= " << d
                << "\t E(" << iter << ")= " << E << endl;
        c *= 0.5;
        rc *= 0.5;
        u += c;
        d = A.Residual(u, r);
        E = A.Energy(u);
        if (E < E_old + sigma * rc) break;
      }
    }
    if (E > E_old + sigma * rc) {
      vout(5) << "line search unsuccessful." << endl;
      ++LS_cnt;
      if (LS_cnt == 3) {
        vout(1) << "too many line searches unsuccessful." << endl;
        iter = max_iter;
      }
    }
  }
  vout(0) << "d(" << iter << ")= " << d << "\t E(" << iter << ")= " << E
          << " rate " << rate() << endl;
  mout.EndBlock();
}

bool NewtonMethod(IAssemble &assemble, Vector &u) {
  Newton newton;
  newton.PrintInfo();
  assemble.PrintInfo();
  newton.operator()(assemble, u);
  return newton.converged();
}

bool NewtonMethod(IAssemble *assemble, Vector &u) {
  return NewtonMethod(*assemble, u);
}