From b5b78976e5cbf9dd18d67b2412f03fb28ef2fb87 Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Tue, 23 Aug 2022 18:16:49 +0300
Subject: [PATCH] edit
---
src/action.hpp | 13 ++++++---
src/core.cpp | 13 +++++----
src/xtest.f95 | 24 ++++++++++++----
usadelndsoc/solver.py | 65 ++++++++++++++++++++++++++++++++-----------
4 files changed, 84 insertions(+), 31 deletions(-)
diff --git a/src/action.hpp b/src/action.hpp
index c58cb3b..b93cade 100644
--- a/src/action.hpp
+++ b/src/action.hpp
@@ -60,7 +60,8 @@ inline Scalar S_2(const array::ArrayView<const Mask, Shape0>& mask,
const array::ArrayView<ConstScalar, Shape>& Qmat,
const array::ArrayView<const Complex, Shape2>& U,
const double& Lx,
- const double& Ly)
+ const double& Ly,
+ const Complex& D)
{
const int nx = Qmat.dim(0);
const int ny = Qmat.dim(1);
@@ -74,6 +75,9 @@ inline Scalar S_2(const array::ArrayView<const Mask, Shape0>& mask,
const size_t invdir[4] = {2, 3, 1, 0};
Scalar S{0.0};
+ if (D == Complex(0))
+ return S;
+
for (int i = 0; i < nx; ++i) {
for (int j = 0; j < ny; ++j) {
if (mask(i,j) == Mask::VACUUM)
@@ -102,7 +106,7 @@ inline Scalar S_2(const array::ArrayView<const Mask, Shape0>& mask,
}
}
- return (dd[0]*dd[1]*2.0/4.0) * S;
+ return (D*dd[0]*dd[1]*2.0/4.0) * S;
}
template <typename ConstScalar, typename Shape0, typename Shape, typename Shape2, typename Scalar = typename std::remove_const<ConstScalar>::type>
@@ -245,7 +249,8 @@ inline Scalar S(const array::ArrayView<const Mask, Shape0>& mask,
const array::ArrayView<const Complex, Shape3>& Omega,
const double& Lx,
const double& Ly,
- const Complex& alpha)
+ const Complex& alpha,
+ const Complex& D)
{
using array::Array;
using array::Shape;
@@ -283,7 +288,7 @@ inline Scalar S(const array::ArrayView<const Mask, Shape0>& mask,
}
}
- S = S + S_2(mask, Qmat, U, Lx, Ly) + S_Omega(mask, Qmat, omega, Omega, Lx, Ly) + S_H(mask, Qmat, U, Lx, Ly, alpha);
+ S = S + S_2(mask, Qmat, U, Lx, Ly, D) + S_Omega(mask, Qmat, omega, Omega, Lx, Ly) + S_H(mask, Qmat, U, Lx, Ly, alpha);
return S;
}
diff --git a/src/core.cpp b/src/core.cpp
index 022d90f..e98f4ba 100644
--- a/src/core.cpp
+++ b/src/core.cpp
@@ -52,6 +52,7 @@ private:
size_t nx_, ny_;
double Lx_, Ly_;
Complex alpha_;
+ Complex D_;
Complex omega_;
typedef Vector<size_t> SizeVector;
@@ -87,8 +88,8 @@ private:
}
public:
- Action(double Lx, double Ly, Complex alpha, MaskType&& mask, UType&& U, OmegaType&& Omega)
- : U_(std::move(U)), Omega_(std::move(Omega)), mask_(std::move(mask)), nx_(U_.dim(0)), ny_(U_.dim(1)), Lx_(Lx), Ly_(Ly), alpha_(alpha), omega_(0),
+ Action(double Lx, double Ly, Complex alpha, Complex D, MaskType&& mask, UType&& U, OmegaType&& Omega)
+ : U_(std::move(U)), Omega_(std::move(Omega)), mask_(std::move(mask)), nx_(U_.dim(0)), ny_(U_.dim(1)), Lx_(Lx), Ly_(Ly), alpha_(alpha), D_(D), omega_(0),
hes_pattern_map_()
{
if (Omega.dim(0) != nx_ || Omega.dim(1) != ny_)
@@ -107,7 +108,7 @@ public:
{
if (Q.dim(0) != nx_ || Q.dim(1) != ny_)
throw std::out_of_range("Q array wrong size");
- return S(mask_, Q, omega_, U_, Omega_, Lx_, Ly_, alpha_);
+ return S(mask_, Q, omega_, U_, Omega_, Lx_, Ly_, alpha_, D_);
}
CppAD::ADFun<Complex> eval_ad(const QType& Qval)
@@ -129,7 +130,7 @@ public:
CppAD::Independent(Q.storage(), dynamic);
Array<ADComplex, Shape<1> > res;
- res(0u) = S(mask_, Q, dynamic[0], U_, Omega_, Lx_, Ly_, alpha_);
+ res(0u) = S(mask_, Q, dynamic[0], U_, Omega_, Lx_, Ly_, alpha_, D_);
return {Q.storage(), res.storage()};
}
@@ -294,8 +295,8 @@ PYBIND11_MODULE(_core, m)
m.doc() = "usadelndsoc._core";
py::class_<Action>(m, "Action")
- .def(py::init<double, double, Complex, Action::MaskType, Action::UType, Action::OmegaType>(), "Action",
- py::arg("L_x"), py::arg("L_y"), py::arg("alpha"), py::arg("mask"), py::arg("U"), py::arg("Omega"))
+ .def(py::init<double, double, Complex, Complex, Action::MaskType, Action::UType, Action::OmegaType>(), "Action",
+ py::arg("L_x"), py::arg("L_y"), py::arg("alpha"), py::arg("D"), py::arg("mask"), py::arg("U"), py::arg("Omega"))
.def("eval", &Action::eval, "Evaluate value of action", py::arg("Q"))
.def("grad", &Action::grad, "Evaluate gradient of action", py::arg("Q"))
.def("hess", &Action::hess, "Evaluate hessian of action", py::arg("Q"))
diff --git a/src/xtest.f95 b/src/xtest.f95
index 930c8a1..d7df1da 100644
--- a/src/xtest.f95
+++ b/src/xtest.f95
@@ -1,7 +1,9 @@
program main
use adjac
implicit none
-
+
+ integer, parameter :: n = 35
+
type(adjac_complexan), dimension(:,:,:,:,:), allocatable :: Q
double complex, dimension(:,:,:,:,:), allocatable :: U
type(adjac_complexan), dimension(1) :: res
@@ -9,11 +11,16 @@ program main
double complex, dimension(:, :), allocatable ::jac
double complex :: v
+ double complex, dimension(:), allocatable :: hess_v
+ integer, dimension(:), allocatable :: hess_i
+ integer, dimension(:), allocatable :: hess_j
+ integer :: nnz
+
integer i1, i2, i3, i4, i5, p, p2
- allocate(Q(50,50,2,2,2))
- allocate(U(50,50,4,4,4))
- allocate(jac(1, 50*50*2*2*2))
+ allocate(Q(n,n,2,2,2))
+ allocate(U(n,n,4,4,4))
+ allocate(jac(1, n*n*2*2*2))
call adjac_reset()
@@ -44,12 +51,19 @@ program main
res(1) = S_2(Q, U, 0.1d0)
call adjac_get_dense_jacobian(res, jac)
- call adjac_reset(.false.)
do i1 = 1, 10
write(*,*) '->', jac(1,i1)
end do
+ call adjac_get_coo_hessian(res(1), nnz, hess_v, hess_i, hess_j)
+
+ do i1 = 1, nnz
+ write(*,*) hess_i(i1), hess_j(i1), '->>', hess_v(i1)
+ end do
+
+ call adjac_reset(.false.)
+
return
contains
function inverse(M) result(W)
diff --git a/usadelndsoc/solver.py b/usadelndsoc/solver.py
index fa32924..a7d50af 100644
--- a/usadelndsoc/solver.py
+++ b/usadelndsoc/solver.py
@@ -75,10 +75,11 @@ def _norm(z):
class Core:
- def __init__(self, Lx=1, Ly=1, alpha=0, nx=1, ny=1):
- self.Lx = 1
- self.Ly = 1
- self.alpha = 0
+ def __init__(self, Lx=1, Ly=1, alpha=0, D=1, nx=1, ny=1):
+ self.Lx = Lx
+ self.Ly = Ly
+ self.alpha = alpha
+ self.D = D
self.set_shape(nx, ny)
def set_shape(self, nx=1, ny=1):
@@ -95,7 +96,9 @@ class Core:
def new_action(self):
self.init_U()
- return Action(self.Lx, self.Ly, self.alpha, self._mask, self._U, self._Omega)
+ return Action(
+ self.Lx, self.Ly, self.alpha, self.D, self._mask, self._U, self._Omega
+ )
@property
def x(self):
@@ -185,7 +188,7 @@ class Solver:
if solver is None:
solver = "lgmres"
- if A_size > 20000:
+ if A_size > 40000:
_log_solve.debug(
f" Problem too large ({A_size}): cannot compute hessian"
)
@@ -377,7 +380,18 @@ class Solver:
return (abs(Phi + dx) < 1.0).all()
def _check_final(self, Phi):
- # XXX: check the branch here!
+ # Check branch
+ g = Phi[:, :, 0]
+ gt = Phi[:, :, 1]
+ I = np.eye(2)
+ s = np.sign(self._action.omega)
+
+ G = s * np.linalg.solve(I + g @ gt, I + gt @ g)
+ if G[:, :, 0, 0].real.min() < -0.05 or G[:, :, 1, 1].real.min() < -0.05:
+ Phi[:, :, 0] *= -0.5
+ Phi[:, :, 1] *= 0.5
+ _log_solve.info(" Force restart: wrong branch")
+ return False
return True
def _eval_rhs(self, Phi):
@@ -429,7 +443,7 @@ class Solver:
if sz > 1000 and pyamg is not None:
preconditioner = "pyamg"
elif umfpack is not None:
- preconditioner = "umfpack-ilu"
+ preconditioner = "umfpack" if sz < 100 else "umfpack-ilu"
else:
preconditioner = "splu"
@@ -461,7 +475,7 @@ class Solver:
def solve_many(self, omega, **kw):
omega = np.asarray(omega)
- Phi = np.zeros(self.Phi.shape + omega.shape, dtype=self.Phi.dtype)
+ Phi = np.zeros(omega.shape + self.Phi.shape, dtype=self.Phi.dtype)
self._M = None
self._Phi[...] = 0
@@ -474,13 +488,18 @@ class Solver:
for v in it:
self.solve(omega=v, action=action, **kw)
Phi_prev = self._Phi.copy()
- Phi[(Ellipsis,) + it.multi_index] = self._Phi
+ Phi[it.multi_index + (Ellipsis,)] = self._Phi
if np.isnan(self._Phi).any():
self._Phi[...] = Phi_prev
return Result(self, omega=omega, Phi=Phi)
+def _bcast_left(x, shape):
+ x = np.asarray(x)
+ return x[(Ellipsis,) + (None,) * (len(shape) - x.ndim)]
+
+
class Result:
def __init__(self, parent, omega=None, Phi=None, core=None):
self._core = core if core is not None else parent._core
@@ -494,26 +513,40 @@ class Result:
@property
def G(self):
- s = np.sign(self.omega)
+ g = self.Phi[..., 0, :, :]
+ gt = self.Phi[..., 0, :, :]
+ I = np.eye(2)
+ s = _bcast_left(np.sign(self.omega), g.shape)
with np.errstate(invalid="ignore"):
- return (1 - self.Phi[0] * self.Phi[1]) / (1 + self.Phi[0] * self.Phi[1]) * s
+ return s * np.linalg.solve(I + g @ gt, I - g @ gt)
@property
def F(self):
- s = np.sign(self.omega)
+ g = self.Phi[..., 0, :, :]
+ gt = self.Phi[..., 0, :, :]
+ I = np.eye(2)
+ s = _bcast_left(np.sign(self.omega), g.shape)
with np.errstate(invalid="ignore"):
- return 2 * self.Phi[0] / (1 + self.Phi[0] * self.Phi[1]) * s
+ return 2 * s * np.linalg.solve(I + g @ gt, g)
@property
def Fc(self):
+ g = self.Phi[..., 0, :, :]
+ gt = self.Phi[..., 0, :, :]
+ I = np.eye(2)
+ s = _bcast_left(np.sign(self.omega), g.shape)
with np.errstate(invalid="ignore"):
- return 2 * self.Phi[1] / (1 + self.Phi[0] * self.Phi[1])
+ return 2 * s * np.linalg.solve(I + gt @ g, gt)
@property
def g(self):
+ g = self.Phi[..., 0, :, :]
+ gt = self.Phi[..., 0, :, :]
+ I = np.eye(2)
+ s = _bcast_left(np.sign(self.omega), g.shape)
with np.errstate(invalid="ignore"):
G = self.G
- return np.array([[G, self.F], [self.Fc, -G]])
+ return s * np.linalg.solve(I + gt @ g, I - gt @ g)
@property
def J(self):
--
GitLab