1 files changed, 28 insertions, 28 deletions

diff --git a/stokes.hpp b/stokes.hpp
index 9fd6f01..8c0c1fb 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,9 +1,10 @@
 #pragma once
 
 #include "p-spin.hpp"
-#include "complex_normal.hpp"
 #include "dynamics.hpp"
 
+#include <iostream>
+
 template <class Scalar, int ...ps>
 class Cord {
 public:
@@ -17,7 +18,7 @@ public:
     Scalar H2 = M.getHamiltonian(z2);
     Scalar H3 = M.getHamiltonian(z3);
 
-    if (real(H2) > real(H3)) {
+    if (std::real(H2) > std::real(H3)) {
       z0 = z2;
       z1 = z3;
     } else {
@@ -62,7 +63,7 @@ public:
     Vector<Scalar> ż = zDot(z, dH);
     Vector<Scalar> dz = df(t);
 
-    return 1 - real(ż.dot(dz)) / ż.norm() / dz.norm();
+    return 1 - std::real(ż.dot(dz)) / ż.norm() / dz.norm();
   }
 
   Real totalCost(unsigned nt) const {
@@ -117,24 +118,22 @@ public:
       Real dfdgn = dgCoeff(i, t);
       Vector<Scalar> dC = - 0.5 / ż.norm() / dz.norm() * (
           dfdgn * ż.conjugate() + fdgn * dżc * dz + fdgn * dż * dz.conjugate()
-          - real(dz.dot(ż)) * (
+          - std::real(dz.dot(ż)) * (
               dfdgn * dz.conjugate() / dz.squaredNorm() +
               fdgn * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm()
             )
         );
 
-      for (unsigned j = 0; j < dC.size(); j++) {
-        x(z.size() * i + j) = dC(j);
-        x(z.size() * gs.size() + z.size() * i + j) = conj(dC(j));
-      }
+      x.segment(z.size() * i, z.size()) = dC;
+      x.segment(z.size() * gs.size() + z.size() * i, z.size()) = dC.conjugate();
 
       for (unsigned j = 0; j < gs.size(); j++) {
         Real fdgm = gCoeff(j, t);
         Real dfdgm = dgCoeff(j, t);
-        Scalar CC = - real(dz.dot(ż)) / ż.norm() / dz.norm();
+        Scalar CC = - std::real(dz.dot(ż)) / ż.norm() / dz.norm();
         Vector<Scalar> dCm = - 0.5 / ż.norm() / dz.norm() * (
             dfdgm * ż.conjugate() + fdgm * dżc * dz + fdgm * dż * dz.conjugate()
-            - real(dz.dot(ż)) * (
+            - std::real(dz.dot(ż)) * (
                 dfdgm * dz.conjugate() / dz.squaredNorm() +
                 fdgm * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm()
               )
@@ -195,21 +194,17 @@ public:
              )
             ;
 
-        for (unsigned k = 0; k < z.size(); k++) {
-          for (unsigned l = 0; l < z.size(); l++) {
-            M(z.size() * i + k, z.size() * j + l) = dcdC(k, l);
-            M(z.size() * gs.size() + z.size() * i + k, z.size() * j + l) = conj(ddC(k, l));
-            M(z.size() * i + k, z.size() * gs.size() + z.size() * j + l) = ddC(k, l);
-            M(z.size() * gs.size() + z.size() * i + k, z.size() * gs.size() + z.size() * j + l) = conj(dcdC(k, l));
-          }
-        }
+        M.block(z.size() * i, z.size() * j, z.size(), z.size()) = dcdC;
+        M.block(z.size() * gs.size() + z.size() * i, z.size() * j, z.size(), z.size()) = ddC.conjugate();
+        M.block(z.size() * i, z.size() * gs.size() + z.size() * j, z.size(), z.size()) = ddC;
+        M.block(z.size() * gs.size() + z.size() * i, z.size() * gs.size() + z.size() * j, z.size(), z.size()) = dcdC.conjugate();
       }
     }
 
     return {x, M};
   }
 
-  std::pair<Real, Real> relaxStepNewton(unsigned nt, Real δ₀) {
+  std::tuple<Real, Real> relaxStepNewton(unsigned nt, Real δ₀) {
     Vector<Scalar> dgsTot = Vector<Scalar>::Zero(2 * z0.size() * gs.size());
     Matrix<Scalar> HessTot = Matrix<Scalar>::Zero(2 * z0.size() * gs.size(), 2 * z0.size() * gs.size());
 
@@ -231,32 +226,37 @@ public:
 
     Matrix<Scalar> I = abs(HessTot.diagonal().array()).matrix().asDiagonal();
 
+    Vector<Scalar> δg(gs.size() * z0.size());
+
     while (newCost > oldCost) {
-      Vector<Scalar> δg = (HessTot + δ * I).partialPivLu().solve(dgsTot);
+      δg = (HessTot + δ * I).partialPivLu().solve(dgsTot).tail(gs.size() * z0.size());
 
       for (unsigned i = 0; i < gs.size(); i++) {
-        for (unsigned j = 0; j < z0.size(); j++) {
-          cNew.gs[i](j) = gs[i](j) - δg[z0.size() * gs.size() + z0.size() * i + j];
-        }
+        cNew.gs[i] = gs[i] - δg.segment(z0.size() * i, z0.size());
       }
 
       newCost = cNew.totalCost(nt);
 
-      δ *= 2;
+      δ *= 1.5;
     }
 
+    std::cerr << newCost << " " << stepSize << " " << δ << std::endl;
+
     gs = cNew.gs;
 
-    return {δ / 2, stepSize};
+    return {δ / 1.5, stepSize};
   }
 
-  void relaxNewton(unsigned nt, Real δ₀, unsigned maxSteps) {
+  void relaxNewton(unsigned nt, Real δ₀, Real ε, unsigned maxSteps) {
     Real δ = δ₀;
     Real stepSize = std::numeric_limits<Real>::infinity();
     unsigned steps = 0;
-    while (stepSize > 1e-7 && steps < maxSteps) {
+    while (stepSize > ε && steps < maxSteps) {
       std::tie(δ, stepSize) = relaxStepNewton(nt, δ);
-      δ /= 3;
+      if (δ > 100) {
+        nt++;
+      }
+      δ /= 2;
       steps++;
     }
   }