🎉 Initial commit. Added source files

Author: elisno

perturb.jl

@@ -0,0 +1,41 @@
+using LinearAlgebra
+using DelimitedFiles
+
+N = 8;
+lambda = 0.1;
+H = Array(Diagonal(rand(ComplexF64,N)));
+
+# For comparison with Fortran, add anharmonic term in Julia
+begin
+    X = zeros(ComplexF64,N,N)
+    for j in 1:N
+        for i in 1:N
+            if abs(i-j) == 1
+                X[i,j]= complex(0.5 * sqrt(i+j+1))
+            end
+        end
+    end
+    native = H + lambda * X^4
+end
+
+function perturb!(N,lambda,H)
+    ccall((:perturb_,"./libperturb.so"), Cvoid, (Ref{Int64}, Ref{ComplexF64}, Ref{ComplexF64}), N, lambda, H)
+    return H
+end
+
+perturb!(N,lambda,H)
+
+
+# Write results to file
+begin
+    open("fortrancall.txt","w") do f
+        writedlm(f,H)
+    end 
+
+    open("julia.txt","w") do f
+        writedlm(f,native)
+    end
+end
+
+# Norm of residuals
+println("norm(H - native) = ", norm(H - native))

perturbhamiltonian.f90

@@ -0,0 +1,31 @@
+
+SUBROUTINE PERTURB(N, lambda, H)
+    integer*8, intent(in) :: N
+    complex*16, intent(in):: lambda
+    complex*16, intent(out) :: H(N,N)
+    complex*16, ALLOCATABLE, DIMENSION(:,:) :: X, X2, X4
+    complex*16                              :: czero=0.0_16
+
+    integer :: i, j
+
+
+    ALLOCATE(X(N,N), X2(N,N), X4(N,N))
+    X(:,:) = czero
+    X2(:,:) = czero
+    X4(:,:) = czero
+
+    do j=1,size(X,2) !n
+        do i=1,size(X,1) !n
+            if (abs(i-j) == 1) then
+                X(i,j)= complex(0.5* sqrt(real(i) + real(j) + 1.0), 0.0)
+            end if
+        end do
+    end do
+
+    X2 = matmul(X, X)
+    X4 = matmul(X2, X2)
+
+    H = H + lambda * X4
+
+    DEALLOCATE(X, X2, X4)
+END subroutine PERTURB