1 files changed, 82 insertions, 0 deletions
diff --git a/src/cormaddp.f90 b/src/cormaddp.f90
new file mode 100644
index 0000000..f9a0f69
--- /dev/null
+++ b/src/cormaddp.f90
@@ -0,0 +1,82 @@
+   subroutine cormadvecdp(matrix,nrow,ncol,res,ressize,correcteven)
+
+      use selectionalgo
+
+      implicit none
+      
+      !note kind dp is defined in selectionalgo
+      integer, intent(in) :: nrow, ncol,ressize, correcteven     !.Fortran R function will not deal with deferred size arrays
+      real(kind=dp), dimension(nrow,ncol), intent(inout) :: matrix
+      real(kind=dp), dimension(ressize), intent(out) :: res
+
+      integer :: i, j, n, p
+      real(kind=dp)  :: med,mad, fresh
+      real(kind=dp), dimension(:), allocatable :: U, V, A, B
+
+
+      p=ubound(matrix,2)
+      n=ubound(matrix,1)
+      allocate(U(n),V(n))
+
+      if (correcteven==1) then
+         do i=1,p-3,3
+            U=matrix(:,i)
+            med=iselect(U,evencorrection=.true.)
+            matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.))
+            U=matrix(:,i+1)
+            med=iselect(U,evencorrection=.true.)
+            matrix(:,i+1)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.))
+            U=matrix(:,i+2)
+            med=iselect(U,evencorrection=.true.)
+            matrix(:,i+2)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.))
+         end do
+         do i=i,p
+            U=matrix(:,i)
+            med=iselect(U,evencorrection=.true.)
+            matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.))
+         end do
+         !unrolled loops
+         n=1
+         do i=1,p-1
+            do j=i+1,p
+               U=matrix(:,i)+matrix(:,j)
+               V=-matrix(:,i)+matrix(:,j)
+               mad=(cost*iselect(abs(U-iselect(U,evencorrection=.true.)),evencorrection=.true.))**2
+               med=(cost*iselect(abs(V-iselect(V,evencorrection=.true.)),evencorrection=.true.))**2
+               res(n)=(mad-med)/(mad+med)
+               n=n+1
+            end do 
+         end do
+      else
+         do i=1,p-3,3
+            U=matrix(:,i)
+            med=iselect(U)
+            matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med)))
+            U=matrix(:,i+1)
+            med=iselect(U)
+            matrix(:,i+1)=(U-med)/(sqrt2*cost*iselect(abs(U-med)))
+            U=matrix(:,i+2)
+            med=iselect(U)
+            matrix(:,i+2)=(U-med)/(sqrt2*cost*iselect(abs(U-med)))
+         end do
+         do i=i,p
+            U=matrix(:,i)
+            med=iselect(U)
+            matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med)))
+         end do
+         !unrolled loops
+         n=1
+         do i=1,p-1
+            do j=i+1,p
+               U=matrix(:,i)+matrix(:,j)
+               V=-matrix(:,i)+matrix(:,j)
+               mad=(cost*iselect(abs(U-iselect(U))))**2
+               med=(cost*iselect(abs(V-iselect(V))))**2
+               res(n)=(mad-med)/(mad+med)
+               n=n+1
+            end do 
+         end do
+      end if
+
+
+   end subroutine cormadvecdp