remove the contigious keyword. triggers multiple copies.

Libraries/libmscbdecomp/Exponentials.F90

@@ -214,7 +214,7 @@ end subroutine FullExp_vecmult_T
 
 subroutine FullExp_lmult(this, mat)
     class(FullExp) :: this
-    complex(kind=kind(0.D0)), intent(inout), contiguous :: mat(:,:)
+    complex(kind=kind(0.D0)), intent(inout) :: mat(:,:)
     integer :: i
     do i = this%evals-1, 1, -2
        call this%stages(i+1)%lmult_T(mat)
@@ -234,7 +234,7 @@ end subroutine FullExp_adjoint_over_two
 
 subroutine FullExp_lmultinv(this, mat)
     class(FullExp) :: this
-    complex(kind=kind(0.D0)), intent(inout), contiguous :: mat(:,:)
+    complex(kind=kind(0.D0)), intent(inout) :: mat(:,:)
     integer :: i
     do i = 1, this%evals, 2
        call this%stages(i)%lmultinv(mat)
@@ -290,7 +290,7 @@ end subroutine EulerExp_dealloc
 !> @brief 
 !> This function multiplies this Euler exponential with a vector.
 !
-!> @param[in] this The exponential opbject
+!> @param[in] this The exponential object
 !> @param[in] vec The vector that we multiply
 !--------------------------------------------------------------------
 subroutine EulerExp_vecmult(this, vec)
@@ -324,7 +324,7 @@ end subroutine EulerExp_vecmult_T
 
 subroutine EulerExp_lmultinv(this, mat)
     class(EulerExp) :: this
-    complex(kind=kind(0.D0)), dimension(:, :), contiguous :: mat
+    complex(kind=kind(0.D0)), dimension(:, :) :: mat
     integer :: i
     do i = 1, this%nrofcols
         call this%singleexps(i)%dat%lmultinv(mat)
@@ -333,7 +333,7 @@ end subroutine EulerExp_lmultinv
 
 subroutine EulerExp_lmult(this, mat)
     class(EulerExp) :: this
-    complex(kind=kind(0.D0)), dimension(:, :), contiguous :: mat
+    complex(kind=kind(0.D0)), dimension(:, :) :: mat
     integer :: i
     do i = this%nrofcols, 1, -1
         call this%singleexps(i)%dat%lmult(mat)

Libraries/libmscbdecomp/GeneralSingleColExp.F90

@@ -63,7 +63,7 @@ contains
 !--------------------------------------------------------------------
 subroutine GeneralSingleColExp_lmultinv(this, mat)
     class(GeneralSingleColExp), intent(in) :: this
-    complex(kind=kind(0.D0)), dimension(:, :), intent(inout), contiguous :: mat
+    complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
     
     call lmultthreeelementbase(this%cinv, this%sinv, this%xy, this%nrofentries, mat)
 end subroutine GeneralSingleColExp_lmultinv

Libraries/libmscbdecomp/HomogeneousSingleColExp.F90

@@ -65,7 +65,7 @@ contains
 !--------------------------------------------------------------------
 subroutine HomogeneousSingleColExp_lmultinv(this, mat)
     class(HomogeneousSingleColExp), intent(in) :: this
-    complex(kind=kind(0.D0)), dimension(:, :), intent(inout), contiguous :: mat
+    complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
     
     call lmultbase(this%cinv, this%sinv, this%xy, this%nrofentries, mat)
 end subroutine HomogeneousSingleColExp_lmultinv

Libraries/libmscbdecomp/SingleColExpBase.F90

@@ -84,7 +84,7 @@ module SingleColExpBase_mod
       subroutine lmultinterface(this, mat)
          import SingleColExpBase
          class(SingleColExpBase), intent(in) :: this
-         Complex(kind=kind(0.d0)), intent(inout),  dimension(:,:), contiguous :: mat
+         Complex(kind=kind(0.d0)), intent(inout),  dimension(:,:) :: mat
       end subroutine
 
     !--------------------------------------------------------------------
@@ -121,7 +121,7 @@ module SingleColExpBase_mod
       subroutine lmultinvinterface(this, mat)
          import SingleColExpBase
          class(SingleColExpBase), intent(in) :: this
-         Complex(kind=kind(0.d0)), intent(inout),  dimension(:,:), contiguous :: mat
+         Complex(kind=kind(0.d0)), intent(inout),  dimension(:,:) :: mat
       end subroutine
 
     !--------------------------------------------------------------------
@@ -228,7 +228,6 @@ end subroutine
 !
 !> Notes: unifying x and y into one array gave some speedup.
 !> Unifying c and s did not...
-!> FIXME: ndim divisible by two...
 !> This is an internal helper function that finds reuse in multiple places.
 !
 !> @param[in] c the diagonal data
@@ -243,20 +242,11 @@ pure subroutine lmultbase(c, s, xy, nrofentries, mat)
     complex (kind=kind(0.d0)), allocatable, intent(in) :: s(:)
     integer, allocatable, intent(in) :: xy(:)
     integer, intent(in) ::nrofentries
-    complex(kind=kind(0.D0)), dimension(:, :), intent(inout), contiguous :: mat
+    complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
     
     integer :: i, j, k, ndim, loopend
     integer, parameter :: step = 2 ! determined to be fastest on 6x6 hubbard
     complex(kind=kind(0.D0)) :: t1(step), t2(step)
-    integer, allocatable, dimension(:) :: xyarray
-    complex(kind=kind(0.D0)), allocatable, dimension(:) :: snh
-    real(kind=kind(0.D0)), allocatable, dimension(:) :: csh
-
-! The intel compiler is really helped by using these temporary arrays
-    allocate(xyarray(size(xy)), csh(nrofentries), snh(nrofentries) )
-    xyarray = xy
-    csh = c
-    snh = s
 
     ndim = size(mat,1)
     loopend = (ndim/step)*step
@@ -266,26 +256,25 @@ pure subroutine lmultbase(c, s, xy, nrofentries, mat)
     do j = 1, loopend, step
         do i = 1, nrofentries! for every matrix
             do k = 1,step
-                t1(k) = mat(xyarray(2*i-1), j+k-1)
-                t2(k) = mat(xyarray(2*i), j+k-1)
+                t1(k) = mat(xy(2*i-1), j+k-1)
+                t2(k) = mat(xy(2*i), j+k-1)
             enddo
             do k = 1, step
-                mat(xyarray(2*i-1), j+k-1) = csh(i) * t1(k) + snh(i) * t2(k)
-                mat(xyarray(2*i), j+k-1) = csh(i) * t2(k) + conjg(snh(i)) * t1(k)
+                mat(xy(2*i-1), j+k-1) = c(i) * t1(k) + s(i) * t2(k)
+                mat(xy(2*i), j+k-1) = c(i) * t2(k) + conjg(s(i)) * t1(k)
             enddo
         enddo
     enddo
-    
+
     ! remainder loop
     if ((ndim - loopend) .ne. 0) then
         do i = 1, nrofentries! for every matrix
-            t1(1) = mat(xyarray(2*i-1), ndim)
-            t2(1) = mat(xyarray(2*i), ndim)
-            mat(xyarray(2*i-1), ndim) = csh(i) * t1(1) + snh(i) * t2(1)
-            mat(xyarray(2*i), ndim) = csh(i) * t2(1) + conjg(snh(i)) * t1(1)
+            t1(1) = mat(xy(2*i-1), ndim)
+            t2(1) = mat(xy(2*i), ndim)
+            mat(xy(2*i-1), ndim) = c(i) * t1(1) + s(i) * t2(1)
+            mat(xy(2*i), ndim) = c(i) * t2(1) + conjg(s(i)) * t1(1)
         enddo
     endif
-    deallocate(xyarray, csh, snh)
 end subroutine
 
 !--------------------------------------------------------------------

Libraries/libmscbdecomp/TraceLessSingleColExp.F90

@@ -87,7 +87,7 @@ pure subroutine lmultthreeelementbase(c, s, x, nrofentries, mat)
     complex (kind=kind(0.d0)), allocatable, intent(in) :: s(:)
     integer, allocatable, intent(in) :: x(:)
     integer, intent(in) ::nrofentries
-    complex(kind=kind(0.D0)), dimension(:, :), intent(inout), contiguous :: mat
+    complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
     
     integer :: i, j, k, ndim, loopend
     integer, parameter :: step = 2 ! determined to be fastest on 6x6 hubbard
@@ -144,7 +144,7 @@ end subroutine
 !--------------------------------------------------------------------
 subroutine TraceLessSingleColExp_lmult(this, mat)
     class(TraceLessSingleColExp), intent(in) :: this
-    complex(kind=kind(0.D0)), dimension(:, :), intent(inout), contiguous :: mat
+    complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
 
     call lmultthreeelementbase(this%c, this%s, this%xy, this%nrofentries, mat)
 end subroutine TraceLessSingleColExp_lmult
@@ -157,12 +157,12 @@ end subroutine TraceLessSingleColExp_lmult
 !> Perform the multiplication of this inverted exponential with a matrix:
 !>  out = this*mat
 !
-!> @param[in] this The exponential that we consider
+!> @param[in] this The exponential that we consider.
 !> @param[inout] mat the matrix that we modify.
 !--------------------------------------------------------------------
 subroutine TraceLessSingleColExp_lmultinv(this, mat)
     class(TraceLessSingleColExp), intent(in) :: this
-    complex(kind=kind(0.D0)), dimension(:, :), intent(inout), contiguous :: mat
+    complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
     integer :: i, j, k, ndim, loopend
     integer, parameter :: step = 2
     complex(kind=kind(0.D0)) :: t1(step), t2(step)

Libraries/libmscbdecomp/ZeroDiagSingleColExp.F90

@@ -75,14 +75,14 @@ end subroutine ZeroDiagSingleColExp_vecmult
 !--------------------------------------------------------------------
 subroutine ZeroDiagSingleColExp_lmult(this, mat)
     class(ZeroDiagSingleColExp), intent(in) :: this
-    complex(kind=kind(0.D0)), dimension(:, :), intent(inout), contiguous :: mat
+    complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
 
     call lmultbase(this%c, this%s, this%xy, this%nrofentries, mat)
 end subroutine ZeroDiagSingleColExp_lmult
 
 subroutine ZeroDiagSingleColExp_lmultinv(this, mat)
     class(ZeroDiagSingleColExp), intent(in) :: this
-    complex(kind=kind(0.D0)), dimension(:, :), intent(inout), contiguous :: mat
+    complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
     integer :: i, j, k, ndim, loopend
     integer, parameter :: step = 2
     complex(kind=kind(0.D0)) :: t1(step), t2(step)