diff --git a/meson.build b/meson.build
index f553bbca6006b5f90636a528cf93fb9a245b2015..677286518d4391d2f32773ba471218be1bd280c3 100644
--- a/meson.build
+++ b/meson.build
@@ -1,4 +1,4 @@
-project('usadelndsoc', ['cpp'],
+project('usadelndsoc', ['cpp', 'fortran'],
version : '0.1',
meson_version: '>= 0.55.0',
default_options : ['cpp_std=c++17', 'build.cpp_std=c++17'])
@@ -43,6 +43,8 @@ pybind11_dep = declare_dependency(
dependencies : [py3_dep]
)
deps = [cppad_dep, eigen3_dep, adept_dep, adolc_dep]
-executable('main', ['src/main.cpp'], dependencies : deps)
+#executable('main', ['src/main.cpp'], dependencies : deps)
+
+executable('main2', ['src/xtest.f95', 'src/adjac/adjac.f95'])
subdir('usadelndsoc')
diff --git a/src/action.hpp b/src/action.hpp
index 63d377bd5dbd35ea283f3f424cf9d66e6e41fa37..62567a7cee3b7afbc0c4b3ce1f58f2358d5d53ed 100644
--- a/src/action.hpp
+++ b/src/action.hpp
@@ -33,8 +33,8 @@ get_Q_matrix(const array::ArrayView<Scalar, array::Shape<2,2,2> >& Q)
Qm.template block<2,2>(0,0) = N * (I - g * gt);
Qm.template block<2,2>(0,2) = N * 2 * g;
- Qm.template block<2,2>(2,0) = -Nt * 2 * gt;
- Qm.template block<2,2>(2,2) = -Nt * (I + gt * g);
+ Qm.template block<2,2>(2,0) = Nt * 2 * gt;
+ Qm.template block<2,2>(2,2) = -Nt * (I - gt * g);
return Qm;
}
@@ -42,8 +42,8 @@ get_Q_matrix(const array::ArrayView<Scalar, array::Shape<2,2,2> >& Q)
template <typename Scalar, typename Shape, typename Shape2, typename ScalarW = typename std::remove_const<Scalar>::type>
inline ScalarW S_2(const array::ArrayView<Scalar, Shape>& Q, const array::ArrayView<const Complex, Shape2>& U, double h)
{
- const int N_i[4] = {0, 1, -1, 0};
- const int N_j[4] = {1, 0, 0, -1};
+ const int N_i[2] = {0, 1};
+ const int N_j[2] = {1, 0};
const size_t invk[4] = {3, 2, 1, 0};
ScalarW S = Scalar(0);
@@ -55,7 +55,7 @@ inline ScalarW S_2(const array::ArrayView<Scalar, Shape>& Q, const array::ArrayV
auto Q1 = get_Q_matrix(Q.part(i,j));
/* Sum over neighbors */
- for (size_t k = 0; k < 4; ++k) {
+ for (size_t k = 0; k < 2; ++k) {
int di = N_i[k];
int dj = N_j[k];
@@ -75,9 +75,12 @@ inline ScalarW S_2(const array::ArrayView<Scalar, Shape>& Q, const array::ArrayV
auto U12 = U.part(i,j,k).matrix();
auto U21 = U.part(i2,j2,invk[k]).matrix();
- auto M = (Q1 - U12*Q2*U21);
+ //auto M = (Q1 - U12*Q2*U21);
+ //S += 2/(4*h) * (M * M).trace();
- S += 2/(4*h) * (M * M).trace();
+ auto M = Q1 * Q2;
+ S += M(2,2);
+ return S;
}
}
}
diff --git a/src/xtest.f95 b/src/xtest.f95
new file mode 100644
index 0000000000000000000000000000000000000000..3e90c5b0f1f0719e9d7bc1e430fe7f20608dc813
--- /dev/null
+++ b/src/xtest.f95
@@ -0,0 +1,160 @@
+program main
+ use adjac
+ implicit none
+
+ type(adjac_complexan), dimension(5,5,2,2,2) :: Q
+ double complex, dimension(5,5,4,4,4) :: U
+ type(adjac_complexan), dimension(1) :: res
+
+ double complex, dimension(1, 5*5*2*2*2) ::jac
+ double complex :: v
+
+ integer i1, i2, i3, i4, i5, p, p2
+
+ call adjac_reset()
+
+ p = 1
+ p2 = 1
+ do i1 = 1, size(Q,1)
+ do i2 = 1, size(Q,2)
+ do i3 = 1, size(Q,3)
+ do i4 = 1, size(Q,4)
+ do i5 = 1, size(Q,5)
+ v = p - 1
+ call adjac_set_independent(Q(i1,i2,i3,i4,i5), v, p)
+ p = p + 1
+ end do
+ end do
+ end do
+ do i3 = 1, size(U,3)
+ do i4 = 1, size(U,4)
+ do i5 = 1, size(U,5)
+ U(i1,i2,i3,i4,i5) = p2 - 1
+ p2 = p2 + 1
+ end do
+ end do
+ end do
+ end do
+ end do
+
+ 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
+
+ return
+contains
+ function inverse(M) result(W)
+ use adjac
+ implicit none
+ type(adjac_complexan), dimension(2,2), intent(in) :: M
+ type(adjac_complexan), dimension(2,2) :: W
+
+ type(adjac_complexan) :: det
+
+ det = M(1,1)*M(2,2) - M(1,2)*M(2,1)
+
+ W(1,1) = M(2,2) / det
+ W(2,2) = M(1,1) / det
+ W(1,2) = -M(1,2) / det
+ W(2,1) = -M(2,1) / det
+ end function inverse
+
+ function get_Q_matrix(Q) result(Qm)
+ use adjac
+ implicit none
+ type(adjac_complexan), dimension(2,2,2), target, intent(in) :: Q
+ type(adjac_complexan), dimension(:,:), pointer :: g
+ type(adjac_complexan), dimension(:,:), pointer :: gt
+ type(adjac_complexan), dimension(4,4) :: Qm
+
+ double complex, dimension(2,2) :: II
+ type(adjac_complexan), dimension(2,2) :: N, Nt
+
+ g => Q(1,:,:)
+ gt => Q(2,:,:)
+
+ II = 0
+ II(1,1) = 1
+ II(2,2) = 1
+
+ N = inverse(II + matmul(g,gt))
+ Nt = inverse(II + matmul(gt,g))
+
+ Qm(1:2,1:2) = matmul(N, II - matmul(g, gt))
+ Qm(1:2,3:4) = 2 * matmul(N, g)
+ Qm(3:4,1:2) = 2 * matmul(Nt, gt)
+ Qm(3:4,3:4) = -matmul(Nt, II - matmul(g, gt))
+ end function get_Q_matrix
+
+ function trace(M) result(res)
+ use adjac
+ implicit none
+ type(adjac_complexan), dimension(4,4), intent(in) :: M
+ type(adjac_complexan) :: res
+
+ integer :: i
+
+ res = 0
+ do i = 1, size(M,1)
+ res = res + M(i,i)
+ end do
+ end function trace
+
+ function S_2(Q, U, h) result(res)
+ use adjac
+ implicit none
+
+ type(adjac_complexan), dimension(:,:,:,:,:), intent(in) :: Q
+ double complex, dimension(:,:,:,:,:), intent(in) :: U
+ double precision, intent(in) :: h
+ type(adjac_complexan) :: res
+
+ integer, dimension(2), parameter :: Ni = (/ 0, 1 /)
+ integer, dimension(2), parameter :: Nj = (/ 1, 0 /)
+ integer, dimension(4), parameter :: invk = (/ 4, 3, 2, 1 /)
+
+ type(adjac_complexan), dimension(4,4) :: M, Q1, Q2
+
+ integer :: i, j, k, i2, j2, nx, ny
+
+ nx = size(Q, 1)
+ ny = size(Q, 2)
+ if (size(Q, 3) .ne. 2 .or. size(Q, 4) .ne. 2 .or. size(Q, 5) .ne. 2) then
+ write(*,*) 'wrong size for Q'
+ stop
+ end if
+ if (size(U, 3) .ne. 4 .or. size(U, 4) .ne. 4 .or. size(U, 5) .ne. 4) then
+ write(*,*) 'wrong size for U'
+ stop
+ end if
+
+ res = 0
+
+ do i = 1, nx
+ do j = 1, ny
+ Q1 = get_Q_matrix(Q(i,j,:,:,:))
+
+ do k = 1, 2
+ i2 = i + Ni(k)
+ j2 = j + Nj(k)
+
+ if (i2 < 1 .or. i2 > nx) cycle
+ if (j2 < 1 .or. j2 > ny) cycle
+
+ Q2 = get_Q_matrix(Q(i2,j2,:,:,:))
+
+ M = Q1 - matmul(U(i,j,k,:,:), matmul(Q2, U(i2,j2,invk(k),:,:)))
+ !res = res + 2/(4*h) * trace(matmul(M, M))
+ M = matmul(Q1, Q2)
+ res = res + M(3,3)
+ return
+ end do
+ end do
+ end do
+ end function S_2
+end program main