diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makevars | 2 | ||||
| -rw-r--r-- | src/RSCdefines.h | 11 | ||||
| -rw-r--r-- | src/cormad_ptr.c | 180 | ||||
| -rw-r--r-- | src/cormaddp.f90 | 642 | ||||
| -rw-r--r-- | src/init.c | 74 | ||||
| -rw-r--r-- | src/selection_algos_ptr.c | 154 |
6 files changed, 406 insertions, 657 deletions
diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..0ee3615 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,2 @@ +PKG_CFLAGS = $(SHLIB_OPENMP_CFLAGS) +PKG_LIBS = $(SHLIB_OPENMP_CFLAGS) diff --git a/src/RSCdefines.h b/src/RSCdefines.h new file mode 100644 index 0000000..75203ac --- /dev/null +++ b/src/RSCdefines.h @@ -0,0 +1,11 @@ +/* Defines quickselect and cormad (double precision) */ + +double quickselect_recursive(double *vector_extract_k, int vec_size, int extract_this_element); + +void cormad(double *matrix, int n_row, int n_col, double *output, int evencorrect); + +#ifdef _OPENMP +void cormad_parallel(double *matrix, int n_row, int n_col, double *output, + int evencorrect, int num_threads); +#endif + diff --git a/src/cormad_ptr.c b/src/cormad_ptr.c new file mode 100644 index 0000000..670a71a --- /dev/null +++ b/src/cormad_ptr.c @@ -0,0 +1,180 @@ +#include "RSCdefines.h" +#include <math.h> +#include <stdio.h> +#ifdef _OPENMP +#include <omp.h> +#endif + +#define SQRT2 (sqrt(2)) +#define CONST 1.4826 + +void cormad(double *matrix, int n_row, int n_col, double *output, + int evencorrect) { + + int k = n_row / 2; // position of the median + int output_size = (n_col - 1) * n_col / 2; + double med, mad; // store variables for medians + double U[n_row], V[n_row]; // help vectors + + /* Transform matrix columns (CORMAD part 1) + * each column has n entries; + * matrix is assumed to be streamed in onevector + */ + int i = 0; + for (int l = 0; l < n_col * n_row; l++) { + if (i == (n_row - 1)) { + U[i] = *matrix; + med = quickselect_recursive(U, n_row, k); + if (evencorrect == 1) { // handle even correction + med = (med + quickselect_recursive(U, n_row, k - 1)) / 2; + } + for (int j = 0; j < n_row; j++) { + U[j] = *(matrix - (n_row - 1) + j) - med; + V[j] = fabs(U[j]); + } + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { // handle even correction + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + for (int j = 0; j < n_row; j++) // reassign + *(matrix - (n_row - 1) + j) = U[j] / (SQRT2 * CONST * med); + // prepare next iter + matrix++; + i = 0; + } else { + U[i] = *matrix; + // prepare next iter + matrix++; + i++; + } + } + + matrix = matrix - (n_col * n_row); /* reset pointer */ + + /* Operate on columns pairs (CORMAD part 2) */ + double *matrix_2 = matrix; // used to point at second column + int first_col = 0; /* Running first column */ + int second_col = 0; /* Running second column */ + for (int l = 0; l < output_size; l++) { + if (second_col == n_col - 1) { + first_col++; + second_col = first_col; + /* set pointers to columns */ + matrix += n_row; + matrix_2 = matrix; + } + second_col++; + matrix_2 += n_row; + + for (int i = 0; i < n_row; i++) { // auxiliary vectors from matrix + U[i] = *(matrix + i) + *(matrix_2 + i); + V[i] = -*(matrix + i) + *(matrix_2 + i); + } + mad = quickselect_recursive(U, n_row, k); + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { + mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2; + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + for (int i = 0; i < n_row; i++) { // reassign for new medians + U[i] = fabs(U[i] - mad); + V[i] = fabs(V[i] - med); + } + mad = quickselect_recursive(U, n_row, k); + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { + mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2; + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + mad = pow(CONST * mad, 2); + med = pow(CONST * med, 2); + + // Assign output + *output = (mad - med) / (mad + med); + output++; + } +} + +#ifdef _OPENMP +void cormad_parallel(double *matrix, int n_row, int n_col, double *output, + int evencorrect, int num_threads) { + int k = n_row / 2; // position of the median + int output_size = (n_col - 1) * n_col / 2; + double med, mad; // store variables for medians + double U[n_row], V[n_row]; // help vectors + + /* Transform matrix columns (CORMAD part 1) */ + double *help_matrix = matrix; /* help pointer for matrix */ +#pragma omp parallel for private(med, mad, U, V, help_matrix) \ + num_threads(num_threads) + for (int j = 0; j < n_col; j++) { // iterate on cols + help_matrix = matrix + n_row * j; /* set pointer at beg of col */ + for (int i = 0; i < n_row; i++) + U[i] = *(help_matrix + i); + med = quickselect_recursive(U, n_row, k); + if (evencorrect == 1) { // handle even correction + med = (med + quickselect_recursive(U, n_row, k - 1)) / 2; + } + for (int i = 0; i < n_row; i++) { + U[i] = *(help_matrix + i) - med; + V[i] = fabs(U[i]); + } + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { // handle even correction + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + for (int i = 0; i < n_row; i++) // reassign + *(help_matrix + i) = U[i] / (SQRT2 * CONST * med); + } + + // int l = 0; // used to iterate over output + /* Operate on columns pairs (CORMAD part 2) */ + help_matrix = matrix; /* help pointers for matrix */ + double *help_matrix_2 = matrix; +#pragma omp parallel for num_threads(num_threads) private( \ + med, mad, U, V, help_matrix, help_matrix_2) + for (int l = 0; l < output_size; l++) { + /* Detrmine columns pairs */ + int col1 = 0, col2 = 0; + int copy_l = l + 1; + for (int last_elem = n_col - 1; last_elem > 0; last_elem--) { + copy_l -= last_elem; + if (copy_l <= 0) { + col2 = (n_col - 1) + copy_l; + break; + } else { + col1++; + } + } + /* Set pointers to columns */ + help_matrix = matrix + n_row * col1; // col1 + help_matrix_2 = matrix + n_row * col2; // col2 + + for (int i = 0; i < n_row; i++) { // auxiliary vectors from matrix + U[i] = *(help_matrix + i) + *(help_matrix_2 + i); + V[i] = -*(help_matrix + i) + *(help_matrix_2 + i); + } + mad = quickselect_recursive(U, n_row, k); + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { + mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2; + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + for (int i = 0; i < n_row; i++) { // reassign for new medians + U[i] = fabs(U[i] - mad); + V[i] = fabs(V[i] - med); + } + mad = quickselect_recursive(U, n_row, k); + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { + mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2; + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + mad = pow(CONST * mad, 2); + med = pow(CONST * med, 2); + + // Assign output + *(output + l) = (mad - med) / (mad + med); + } +} +#endif diff --git a/src/cormaddp.f90 b/src/cormaddp.f90 deleted file mode 100644 index b0b6516..0000000 --- a/src/cormaddp.f90 +++ /dev/null @@ -1,642 +0,0 @@ -module selectionalgo - use, intrinsic :: iso_fortran_env - - implicit none - - real, parameter :: cost=1.4826, sqrt2=sqrt(2.) - integer, parameter :: sp = real32, dp=real64 - - interface qselect - procedure quickselectdp, quickselectscalardp - end interface qselect - - interface iselect - procedure introselectdp, introselectscalardp - end interface iselect - - contains - - !!!!!!! DOUBLE PRECISION - recursive subroutine quickselectrecursivedp(invector,tovector,k,output) - real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector - real(kind=dp), intent(inout), pointer, contiguous :: tovector(:) - integer, intent(inout) :: k - real(kind=dp), intent(out) :: output - real(kind=dp) :: swapper, pvt - integer :: i, r, subst - real(kind=dp), pointer :: a,b,c - - r=size(tovector) - - select case(r) - case(1) - output=tovector(1) - return - case(2) - select case(k) - case(1) - output=minval(tovector) - case(2) - output=maxval(tovector) - end select - return - case(3) - select case(k) - case(1) - output=minval(tovector) - case(3) - output=maxval(tovector) - case(2) - i=r/2+1 - output=sum(tovector)-minval(tovector)-maxval(tovector) - end select - return - case default - i=r/2+1 - - !pivoting section (pivot of 3) - a=>tovector(1) - b=>tovector(i) - c=>tovector(r) - pvt=a+b+c - swapper=max(a,b,c) - b=swapper - pvt=pvt-swapper - swapper=min(a,b,c) - a=swapper - pvt=pvt-swapper - c=pvt - - !one pass section - subst=1 - b=>tovector(subst) - do i=1,(r-1-3),3 - a => tovector(i) - if (a<pvt) then - swapper=b - b=a - a=swapper - subst=subst+1 - b=>tovector(subst) - end if - a => tovector(i+1) - if (a<pvt) then - swapper=b - b=a - a=swapper - subst=subst+1 - b=>tovector(subst) - end if - a => tovector(i+2) - if (a<pvt) then - swapper=b - b=a - a=swapper - subst=subst+1 - b=>tovector(subst) - end if - end do - do i=i,(r-1) - a => tovector(i) - if (a<pvt) then - swapper=b - b=a - a=swapper - subst=subst+1 - b=>tovector(subst) - end if - end do - tovector(r)=tovector(subst) - tovector(subst)=pvt - - !decide next step - select case(subst-k) - case(0) - output=tovector(k) - case(1:) - tovector => tovector(1:subst-1) - call quickselectrecursivedp(invector,tovector,k, output) - case(:-1) - k=k-subst - tovector => tovector(subst+1:) - call quickselectrecursivedp(invector,tovector,k, output) - end select - end select - end subroutine quickselectrecursivedp - - recursive subroutine introselectrecursivedp(invector,tovector,k,output) - real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector - real(kind=dp), intent(inout), pointer, contiguous :: tovector(:) - integer, intent(inout) :: k - real(kind=dp), intent(out) :: output - real(kind=dp) :: swapper, pvt - integer, save :: switch - integer :: i, r, subst, r_min, switch_after - real(kind=dp), pointer :: a,b,c - - r=ubound(tovector,1) - r_min = 3000 ! at least this bigger to switch to median of medians - switch_after = 5 ! at least fails this amount of time to reorder - ! one third of the vector to consider med of med - - select case(r) - case(1) - output=tovector(1) - switch=0 - return - case(2) - select case(k) - case(1) - output=minval(tovector) - case(2) - output=maxval(tovector) - end select - switch=0 - return - case(3) - select case(k) - case(1) - output=minval(tovector) - case(3) - output=maxval(tovector) - case(2) - i=r/2+1 - output=sum(tovector)-minval(tovector)-maxval(tovector) - end select - switch=0 - return - case default - i=r/2+1 - - !pivoting section (pivot of 3) - a=>tovector(1) - b=>tovector(i) - c=>tovector(r) - pvt=a+b+c - swapper=max(a,b,c) - b=swapper - pvt=pvt-swapper - swapper=min(a,b,c) - a=swapper - pvt=pvt-swapper - c=pvt - - !one pass section - subst=1 - b=>tovector(subst) - do i=1,(r-1-3),3 - a => tovector(i) - if (a<pvt) then - swapper=b - b=a - a=swapper - subst=subst+1 - b=>tovector(subst) - end if - a => tovector(i+1) - if (a<pvt) then - swapper=b - b=a - a=swapper - subst=subst+1 - b=>tovector(subst) - end if - a => tovector(i+2) - if (a<pvt) then - swapper=b - b=a - a=swapper - subst=subst+1 - b=>tovector(subst) - end if - end do - do i=i,(r-1) - a => tovector(i) - if (a<pvt) then - swapper=b - b=a - a=swapper - subst=subst+1 - b=>tovector(subst) - end if - end do - tovector(r)=tovector(subst) - tovector(subst)=pvt - - !decide next step - select case(subst-k) - case(0) - output=tovector(k) - switch=0 - case(1:) - tovector => tovector(1:subst-1) - if ((subst-1)>r*2/3) then - switch=switch+1 - if (switch==switch_after .and. r>r_min) then - !reset switch and call medofmeddp - switch=0 - call medofmeddp(invector,tovector,k, output) - else - call introselectrecursivedp(invector,tovector,k, output) - end if - else - switch=0 - call introselectrecursivedp(invector,tovector,k, output) - end if - case(:-1) - k=k-subst - tovector => tovector(subst+1:) - if ((r-subst)>r*2/3) then - switch=switch+1 - if (switch==switch_after .and. r>r_min) then - !reset switch and call medofmeddp - switch=0 - call medofmeddp(invector,tovector,k, output) - else - call introselectrecursivedp(invector,tovector,k, output) - end if - else - switch=0 - call introselectrecursivedp(invector,tovector,k, output) - end if - end select - end select - end subroutine introselectrecursivedp - - recursive subroutine medofmeddp(invector,tovector,k,output) - real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector - real(kind=dp), intent(inout), pointer, contiguous :: tovector(:) - real(kind=dp), intent(out) :: output - integer, intent(inout) :: k - - real(kind=dp), pointer, contiguous :: subsec(:) - real(kind=dp), dimension(:),allocatable :: meds - integer :: i,subst,r - real(kind=dp) :: pvt - !real(kind=dp) :: swapper - !real(kind=dp), pointer :: a,b - - r=ubound(tovector,1) - select case(mod(r,5)) - case(0) - allocate(meds(r/5)) - do i=1,r-4,5 - subsec=>tovector(i:i+4) - subst=minloc(subsec,1) - pvt=subsec(subst) - subsec(subst)=subsec(1) - subsec(1)=pvt - subst=minloc(subsec(2:),1)+1 - pvt=subsec(subst) - subsec(subst)=subsec(2) - subsec(2)=pvt - subst=minloc(subsec(3:),1)+2 - pvt=subsec(subst) - subsec(subst)=subsec(3) - subsec(3)=pvt - meds(int(i/5.0)+1)=pvt - end do - case(1) - allocate(meds(r/5+1)) - do i=1,r-4,5 - subsec=>tovector(i:i+4) - subst=minloc(subsec,1) - pvt=subsec(subst) - subsec(subst)=subsec(1) - subsec(1)=pvt - subst=minloc(subsec(2:),1)+1 - pvt=subsec(subst) - subsec(subst)=subsec(2) - subsec(2)=pvt - subst=minloc(subsec(3:),1)+2 - pvt=subsec(subst) - subsec(subst)=subsec(3) - subsec(3)=pvt - meds(int(i/5.0)+1)=pvt - end do - meds(r/5+1)=tovector(i) - case(2) - allocate(meds(r/5+1)) - do i=1,r-4,5 - subsec=>tovector(i:i+4) - subst=minloc(subsec,1) - pvt=subsec(subst) - subsec(subst)=subsec(1) - subsec(1)=pvt - subst=minloc(subsec(2:),1)+1 - pvt=subsec(subst) - subsec(subst)=subsec(2) - subsec(2)=pvt - subst=minloc(subsec(3:),1)+2 - pvt=subsec(subst) - subsec(subst)=subsec(3) - subsec(3)=pvt - meds(int(i/5.0)+1)=pvt - end do - subsec=>tovector(i:) - subst=minloc(subsec,1) - pvt=subsec(subst) - subsec(subst)=subsec(1) - subsec(1)=pvt - meds(r/5+1)=subsec(2) - case(3) - allocate(meds(r/5+1)) - do i=1,r-4,5 - subsec=>tovector(i:i+4) - subst=minloc(subsec,1) - pvt=subsec(subst) - subsec(subst)=subsec(1) - subsec(1)=pvt - subst=minloc(subsec(2:),1)+1 - pvt=subsec(subst) - subsec(subst)=subsec(2) - subsec(2)=pvt - subst=minloc(subsec(3:),1)+2 - pvt=subsec(subst) - subsec(subst)=subsec(3) - subsec(3)=pvt - meds(int(i/5.0)+1)=pvt - end do - subsec=>tovector(i:) - subst=minloc(subsec,1) - pvt=subsec(subst) - subsec(subst)=subsec(1) - subsec(1)=pvt - subst=minloc(subsec(2:),1)+1 - pvt=subsec(subst) - subsec(subst)=subsec(2) - subsec(2)=pvt - meds(r/5+1)=pvt - case(4) - allocate(meds(r/5+1)) - do i=1,r-4,5 - subsec=>tovector(i:i+4) - subst=minloc(subsec,1) - pvt=subsec(subst) - subsec(subst)=subsec(1) - subsec(1)=pvt - subst=minloc(subsec(2:),1)+1 - pvt=subsec(subst) - subsec(subst)=subsec(2) - subsec(2)=pvt - subst=minloc(subsec(3:),1)+2 - pvt=subsec(subst) - subsec(subst)=subsec(3) - subsec(3)=pvt - meds(int(i/5.0)+1)=pvt - end do - subsec=>tovector(i:) - subst=minloc(subsec,1) - pvt=subsec(subst) - subsec(subst)=subsec(1) - subsec(1)=pvt - subst=minloc(subsec(2:),1)+1 - pvt=subsec(subst) - subsec(subst)=subsec(2) - subsec(2)=pvt - subst=minloc(subsec(3:),1)+2 - pvt=subsec(subst) - subsec(subst)=subsec(3) - subsec(3)=pvt - meds(r/5+1)=pvt - end select - - !set pivot - pvt=qselect(meds) - do i=1,r-3,3 - if (tovector(i)==pvt) then - tovector(i)=tovector(r) - tovector(r)=pvt - exit - else if (tovector(i+1)==pvt) then - tovector(i+1)=tovector(r) - tovector(r)=pvt - exit - else if (tovector(i+2)==pvt) then - tovector(i+2)=tovector(r) - tovector(r)=pvt - exit - else - cycle - end if - end do - do i=i,r - if (tovector(i)==pvt) then - tovector(i)=tovector(r) - tovector(r)=pvt - exit - else - cycle - end if - end do - - call introselectrecursivedp(invector,tovector,k, output) - end subroutine medofmeddp - - function introselectdp(invector,ord,evencorrection) - real(kind=dp), intent(in), dimension(:) :: invector - logical, optional, intent(in) :: evencorrection - integer, intent(in), optional :: ord - real(kind=dp) :: introselectdp - - real(kind=dp), dimension(:), allocatable, target :: vector - real(kind=dp), pointer, contiguous :: tovector(:) => null() - integer :: k1, k - - k1=size(invector) - allocate(vector(k1)) - vector=invector - tovector => vector - - k1=ubound(invector,1) - if(present(ord)) then - k=ord - else - k=k1/2+1 - end if - - if(present(evencorrection)) then - if (k/=k1/2+1 .or. mod(k1,2)/=0) then - !write(*,'(a)') 'Warning: evencorrection argument ignored.' - k1=0 - else - k1=0 - if (evencorrection) k1=k-1 - end if - else - k1=0 - end if - - call introselectrecursivedp(vector,tovector,k, introselectdp) - - if (k1/=0) then - block - real(kind=dp) :: temp - temp=introselectdp - !tovector=>vector(1:k1+1) - tovector=>vector - call introselectrecursivedp(vector,tovector,k1, introselectdp) - introselectdp=(introselectdp+temp)/2.0 - end block - end if - end function introselectdp - - function quickselectdp(invector,ord,evencorrection) - real(kind=dp), intent(in), dimension(:) :: invector - logical, optional, intent(in) :: evencorrection - integer, intent(in), optional :: ord - real(kind=dp) :: quickselectdp - - real(kind=dp), dimension(:), allocatable, target :: vector - real(kind=dp), pointer, contiguous :: tovector(:) => null() - integer :: k1, k - - k1=size(invector) - allocate(vector(k1)) - vector=invector - tovector => vector - - k1=ubound(invector,1) - if(present(ord)) then - k=ord - else - k=k1/2+1 - end if - - if(present(evencorrection)) then - if (k/=k1/2+1 .or. mod(k1,2)/=0) then - !write(*,'(a)') 'Warning: evencorrection argument ignored.' - k1=0 - else - k1=0 - if (evencorrection) k1=k-1 - end if - else - k1=0 - end if - - call quickselectrecursivedp(vector,tovector,k, quickselectdp) - - if (k1/=0) then - block - real(kind=dp) :: temp - temp=quickselectdp - tovector=>vector(1:k1+1) - call quickselectrecursivedp(vector,tovector,k1, quickselectdp) - quickselectdp=(quickselectdp+temp)/2.0 - end block - end if - end function quickselectdp - - !deal with scalar input - function introselectscalardp(invector,ord, evencorrection) - real(kind=dp), intent(in) :: invector - integer, intent(in), optional :: ord - logical, optional, intent(in) :: evencorrection - real(kind=dp) :: introselectscalardp - - if (present(evencorrection) .or. present(ord)) then - continue - end if - - introselectscalardp=invector - end function introselectscalardp - - function quickselectscalardp(invector,ord, evencorrection) - real(kind=dp), intent(in) :: invector - integer, intent(in), optional :: ord - logical, optional, intent(in) :: evencorrection - real(kind=dp) :: quickselectscalardp - - if (present(evencorrection) .or. present(ord)) then - continue - end if - - quickselectscalardp=invector - end function quickselectscalardp - - -end module selectionalgo - -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 - !real(kind=dp) :: fresh - real(kind=dp), dimension(:), allocatable :: U, V - !real(kind=dp), dimension(:), allocatable :: 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 @@ -1,22 +1,66 @@ -#include <stdlib.h> -//#include <stdio.h> +#include "RSCdefines.h" #include <R.h> -#include <Rinternals.h> -#include <Rmath.h> -#include <Rdefines.h> #include <R_ext/RS.h> #include <R_ext/Rdynload.h> +#include <Rdefines.h> +#include <Rinternals.h> +#include <Rmath.h> +#include <stdio.h> +#include <stdlib.h> +#ifdef _OPENMP +#include <omp.h> +#endif + +// Wrapper for C cormad (with pointers) +SEXP cormad_C(SEXP R_matrix, SEXP len_rows, SEXP len_cols, SEXP R_evencorrect, + SEXP R_num_threads) { + // prepare input arguments + int n_row = asInteger(len_rows); + int n_col = asInteger(len_cols); + int output_size = (n_col - 1) * n_col / 2; + int evencorrect = + (n_row % 2 == 1) ? 0 : asInteger(R_evencorrect); // ingore for odd n + int num_threads = asInteger(R_num_threads); + SEXP output = PROTECT(allocVector(REALSXP, output_size)); // define output -extern void F77_NAME(cormadvecdp)(void *, void *, void *, void *, void *, void *); + // create pointers to (duplicate) matrix column + double *dupmat = + REAL(PROTECT(duplicate(R_matrix))); // copy input (to avoid modifications) + + // Call cormad +#ifdef _OPENMP + int max_threads = omp_get_max_threads(); + if (num_threads == 0) /* use max number of threads / 2 */ + num_threads = (max_threads / 2 == 0) ? 1 : max_threads / 2; + else if (num_threads > max_threads) + num_threads = max_threads; + else if (num_threads < 0) + num_threads = + (max_threads + num_threads > 0) ? max_threads + num_threads : 1; + + if (num_threads > 1) { + cormad_parallel(dupmat, n_row, n_col, REAL(output), evencorrect, + num_threads); + } else { + cormad(dupmat, n_row, n_col, REAL(output), evencorrect); + } +#else + if (num_threads != 1) + Rprintf("\nSpecified num_threads=%d, but OPENMP support not found; " + "switching to single core.\n\n", + num_threads); + cormad(dupmat, n_row, n_col, REAL(output), evencorrect); +#endif + // Remove protect to allow R grabage collect + UNPROTECT(2); + return output; +} -static const R_FortranMethodDef FortMethods[] = { - {"cormadvecdp", (DL_FUNC) &F77_NAME(cormadvecdp), 6}, - {NULL, NULL, 0} -}; +static const R_CallMethodDef callMethods[] = { + {"cormad_C", (DL_FUNC)&cormad_C, 5}, {NULL, NULL, 0}}; -void -R_init_RSC(DllInfo *dll) -{ - R_registerRoutines(dll, NULL, NULL, FortMethods, NULL); - R_useDynamicSymbols(dll, FALSE); +void R_init_RSC(DllInfo *dll) { + R_registerRoutines(dll, NULL, callMethods, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); + R_forceSymbols(dll, TRUE); } diff --git a/src/selection_algos_ptr.c b/src/selection_algos_ptr.c new file mode 100644 index 0000000..2abecac --- /dev/null +++ b/src/selection_algos_ptr.c @@ -0,0 +1,154 @@ +#include <stdio.h> +/* Implements functions quickselect and introselect */ + +static double select_corner_cases(double *vector_shorter_3, int size_of_vector, + int extract_this_element_Ccorrected); + +static double pivot_of_3(double *vector_to_be_pivoted, int size); + +static int quickselect_onepass(double *vector, int size); + +/* Functions */ + +double quickselect_recursive(double *vector, int size, int k) { + /* Arguments: + * *vector: pointer to vector + * idN: size of vector + * k: k-th element to be extracted from the vector output: the + * extracted element of order k */ + + /* Decide if we need to go with quickselect_recursive or return output + * (quickselect_recursive works only for vectors with at least 3 elems) */ + if (size < 3) { // select the returning value + return select_corner_cases(vector, size, k); + } + + int subst = quickselect_onepass(vector, size); + + /* Determine where to go next (left of subst or right?) */ + if (subst == k) { + return *(vector + subst); + } else if (subst > k) { // go left + size = subst; + return quickselect_recursive(vector, size, k); + } else { // go right + // readjust k + k -= subst + 1; // take it back to orginal value (for full vector) + vector += subst + 1; + size -= subst + 1; + return quickselect_recursive(vector, size, k); + } +} + +static double select_corner_cases(double *vector, int size, int k) { + /* Used to handel corner cases not handeled by recursive strategy. + * Recursive strategy needs vectors of at least 3 elements + + * Arguments: + * *vector: pointer to vector of doubles to sort + * size: size of the vector + * k: k-th element to be extracted from the vector output: the + * extracted element of order k + */ + + // total elements in vector are idN-id0+1 + double ret = -111; + switch (size) { + case 1: + return *vector; + case 2: + switch (k) { + case 0: // return the smallest of the two + if (*vector < *(vector + 1)) + return *vector; + else + return *(vector + 1); + case 1: // return the biggest of the two + if (*vector > *(vector + 1)) + return *vector; + else + return *(vector + 1); + } + } + return ret; +} + +static double pivot_of_3(double *vector, int size) { + /* perform pivot of 3 returning the median element and + * rearranging the vector so to have: + * id0, i(see below) , idN -> min, max, med + + * Arguments + * vector: vector where to make substitutions + * id0, i, idN: positions of first, median and last element to consider + */ + + int idN = size - 1; // index of last element of the vector + int i = idN / 2; // index of middle element of the vector + double a, b, c; // auxiliary variables (avoid dereferncing too much) + a = *vector; + b = *(vector + i); + c = *(vector + idN); + + /* Nota; vogliamo che l'elemento mediano si trovi alla fine */ + double swapper; + if ((a > b) ^ (a > c)) { // id0 is median element + swapper = c; + c = a; + if (swapper > b) { + a = b; + b = swapper; + } else { + a = swapper; + } + } else if ((b > a) ^ (b > c)) { // i is median element + swapper = c; + c = b; + if (swapper > a) { + b = swapper; + } else { + b = a; + a = swapper; + } + } else { // idN is median element + if (a > b) { + swapper = a; + a = b; + b = swapper; + } + } + *vector = a; + *(vector + i) = b; + *(vector + idN) = c; + return c; +} + +static int quickselect_onepass(double *vector, int size) { + /* perform one pass for quickselect and returns the element where the pivot + * was substituted: this is used to decide whether to go left or right */ + + /* Pivoting section (pivot of 3) */ + double pivot = pivot_of_3(vector, size); + + /* One pass on vector */ + double swapper; + double *subst_el = vector; // address of element to substitute (copy not to + // alter main pointer) + int subst = 0; + for (int i = 0; i < size - 1; i++) { + if (*subst_el < pivot) { + if (subst == i) { + subst++; + } else { + swapper = *subst_el; + *subst_el = *(vector + subst); + *(vector + subst) = swapper; + subst++; + } + } + subst_el++; + } + *subst_el = *(vector + subst); + *(vector + subst) = pivot; + return subst; +} |