diff options
| author | Luca Coraggio <luca.coraggio@unina.it> | 2020-07-04 09:50:03 +0000 |
|---|---|---|
| committer | cran-robot <csardi.gabor+cran@gmail.com> | 2020-07-04 09:50:03 +0000 |
| commit | 511e3ca9e5235e018f772693907d9ec10002b02a (patch) | |
| tree | c7cb699babfa439e6bfbe47007e3916867517f76 /src/cormaddp.f90 | |
version 1.0
Diffstat (limited to 'src/cormaddp.f90')
| -rw-r--r-- | src/cormaddp.f90 | 82 |
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 |