edit

meson.build

-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')

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;
             }
         }
     }

src/xtest.f95

+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