edit

src/core.cpp

@@ -132,9 +132,7 @@ public:
                 compute(Q);
 
             auto jac = f_.Jacobian(x);
-
-            std::array<size_t, 1> shape = {jac.size()/2};
-            return Array<Complex, Shape<Dynamic>, Vector<double> >(std::move(jac), shape);
+            return Array<double, Shape<Dynamic>, Vector<double> >(std::move(jac), {jac.size()});
         }
 
     auto hess(const QType& Q)
@@ -184,7 +182,7 @@ public:
             Array<double, Shape<Dynamic> > jac({n});
             stack.jacobian(jac.data());
 
-            return std::move(jac);
+            return jac;
         }
 };
 

src/xtest.f95

@@ -2,11 +2,11 @@ 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(50,50,2,2,2) :: Q
+  double complex, dimension(50,50,4,4,4) :: U
   type(adjac_complexan), dimension(1) :: res
 
-  double complex, dimension(1, 5*5*2*2*2) ::jac
+  double complex, dimension(1, 50*50*2*2*2) ::jac
   double complex :: v
 
   integer i1, i2, i3, i4, i5, p, p2
@@ -118,6 +118,9 @@ contains
     integer, dimension(2), parameter :: Nj = (/ 1, 0 /)
     integer, dimension(4), parameter :: invk = (/ 4, 3, 2, 1 /)
 
+    type(adjac_complexan), dimension(:,:,:,:), allocatable :: Qmat
+    logical, dimension(:,:), allocatable :: initQ
+
     type(adjac_complexan), dimension(4,4) :: M, Q1, Q2
 
     integer :: i, j, k, i2, j2, nx, ny
@@ -133,11 +136,18 @@ contains
        stop
     end if
 
+    allocate(Qmat(nx,ny,4,4))
+    allocate(initQ(ny,nx))
+    initQ = .false.
+    
     res = 0
 
     do i = 1, nx
        do j = 1, ny
-          Q1 = get_Q_matrix(Q(i,j,:,:,:))
+          if (.not. initQ(i,j)) then
+             Qmat(i,j,:,:) = get_Q_matrix(Q(i,j,:,:,:))
+             initQ(i,j) = .true.
+          end if
 
           do k = 1, 2
              i2 = i + Ni(k)
@@ -146,7 +156,13 @@ contains
              if (i2 < 1 .or. i2 > nx) cycle
              if (j2 < 1 .or. j2 > ny) cycle
 
-             Q2 = get_Q_matrix(Q(i2,j2,:,:,:))
+             if (.not. initQ(i2,j2)) then
+                Qmat(i2,j2,:,:) = get_Q_matrix(Q(i2,j2,:,:,:))
+                initQ(i2,j2) = .true.
+             end if
+
+             Q1 = Qmat(i,j,:,:)
+             Q2 = Qmat(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))