Column sorting

[john_h_16]

Given matrix A of size MxN, I would like to sort column-wise. In Matlab this can be solved easily and quickly by sort(A) and it can even sort row-wise by sort(A,2).

In fortran, I don't find this and thus far I just iteratively sort each column using dlsart2 but this is definitely much slower than MATLAB. I am sure someone must have run into this issue and I hope you can help me speed this up, perhaps along these lines:

1. Is there any column sorting subroutine that works faster then iterating dlsart2, which costs MN log M assuming that each iterate is M log M.

2. Is dlsart2 faster than any other sorting functions? Have anyone compare this it with dpsort.f90 of slatec or some sorting functions from orderpack? I found that dlsart2 is two times faster than qsortd.f90.

3. Is there anything else that is faster than dlsart2 that can beat MATLAB sort function?

John 

[mecej4]

What is 'dlsart2'? Neither Google nor Bing search turns up anything about it. What do you mean by " I just iteratively sort each column"? You can only sort by one selected column -- the column with the sort keys-- if you want to avoid mixing up rows. You can, of course, have secondary, etc. keys in other columns if the primary keys can occur more than once. 

Matlab contains thousands of functions; Fortran has hundreds. There is no intrinsic sorting function in Fortran. Sorting is the topic of an entire volume of Knuth's book set "The Art of Computer Programming", and different sorting algorithms are suited to different sorting problems. You need to define your needs accurately and completely before it makes sense to compare sorting utilities regarding their suitability for your application.

[john_h_16]

Thanks for your response.

dlsart2 is a subroutine of scalapack part of mkl library. Here is the detail documentation: https://software.intel.com/en-us/node/521657

Sorry for the confusion about iterative. What I mean was sort each column one-by-one with a loop. Just to give an example (written in Matlab for simplicity) with matrix A of size M x N:

for j=1:N,  out(:,j) = sort(A(:,j)); end

So, I basically replicate this loop with the dlsart2 whereas in Matlab one can avoid this loop by just doing sort(A,1).

Now what I was asking was really simple: Say if one can't sort column-wise like in MATLAB which means loop is unavoidable. Is there any array sorting subroutine that can be faster than dlsart2? 

[mecej4]

Ah, I see that you spelled the Lapack routine 'dlasrt2' as 'dlsart2', which is why my search failed.

One reason that routines such as dlasrt2 can be slow is that, in addition to sorting the input array, they also return the sort index array. If you have no use for the index array, the sort could be made faster by not building this index.

How big are the columns that you sort, and how many columns do you wish to sort? Are the columns completely random, or do they have some partial order before you sort?

[john_h_16]

The number of columns and the size of the columns will be the same for my application, they can easily be of order 100,000 to 10,000,000.

I am using it to find the k-nearest neighbor, so there should not be in any particular order. For this application, I will need the index too but I don't need all of them. Instead, only the first k sorted array and indices are enough. 

[Andrey_N_Intel]

Hi John,

Please, have a look at the Summary Statistics component of Intel(R) Math Kernel Library, https://software.intel.com/en-us/node/521928. You will need to compute order statistics with routines of this component of the library. For input data set of size n x p where n is number of observations/feature vectors, and p is number of features, the algorithm will sort the columns of the library. Also, if your data is in a different  layout ( p x n ), the algorithm will sort rows of your matrix. 

Additional details how to compute order stats can be found in the Application Notes for Summary Statistics, https://software.intel.com/en-us/mkl_11.1_sslnotes, see section Algorithms and Interfaces in Summary Statistics, subsections Calculating Order Statistics and Estimating Quantiles.

Examples\vslc subdirectory of the library installation directory contains the vsldquantile.c example, which shows how to compute quantiles and order stats. Computation or quantiles can be easily disabled by respective configuration of the example. If you, however, need first k elements in each sorted column of your dataset, then computation of the quantiles might work as well.

Please, let us know, if you have additional questions on the feature of the library.

Thanks,

Andrey

[john_h_16]

Dear Andrey, 

Thanks for your suggestion, it looks like the order statistics is a possibility as you suggested. However, I look at the link you gave me and I can't tell what is computed and how to use it; I am not so familiar with quantile too, I think all I need is just the order statistics. Unfortunately, I don't read C so it wasn't clear to me how to even use this library. Do you have any simple example written in f90?

John

[Andrey_N_Intel]

Hi John,

as the first step to get a general idea how to use API of the library, can you have a quick look at the Application Notes for Summary Statistics, Section "Common Usage Model of Summary Statistics Algorithms"? While the sample code is in C it explains the computational flow with Summary Statistics. Intel(R) MKL provides Fortran examples as well: vsldquantile.f in the examples\vslf subdirectory of the library installation directory is Fortran counterpart of that C example and allows you to compute quantiles and order statistics. Please, learn, build and run it as the next step. If necessary, we will help simplify the example so it would be used to compute order statistics only.

You also might want to have a look at other Fortran examples for Summary Statistics such as vsldstatsum.f or vsldcp.f.

Thanks,

Andrey

[john_h_16]

Andrey,

I finally found this example code vsldquantile.f that you are talking about since we the computer that I am working with has mkl installed. Now, I am having trouble to compile this file, how do you link the library? 

John 

[mecej4]

Please be specific in your problem reports. There is no universal solution for "I am having trouble to compile".

The -mkl -I. option to ifort will suffice on Linux (/Qmkl on Windows), if you copy the example source and the three include files used in it from the .../mkl/examples/vslf/source directory to a working directory.

[Andrey_N_Intel]

Another option to build the example is to use makefile available under vslf folder. You can configure the building process for the examples by providing input parameters into that make file. Please, have a look at the help section in the makefile.

Thanks,

Andrey

[john_h_16]

Dear Andrey and Mecej4

After one day, I finally manage to compile it with 

export MKLROOT=/usr/global/intel/mkl/11.0.0.079/
ifort -mkl -I /usr/global/intel/mkl/11.0.0.079/examples/vslf/source/ vsldquantile.f

Having this means I still have to figure out how to get rid of the quantile that I don't need. It looks like this example includes many .inc files that I am not sure whether I will need it in the end. Looking at the code and the poor documentation, I would not even try to use it now since I am not so sure I can even get the ordering index and I am betting this will be slower than MATLAB sorting function. 

I decide to just use kdtree2 for k-nearest neighbor application which is a well documented code so it is easier to implement (at least for me). In case someone else is reading this thread in the future, I can tell you one drawback with kdtree2 is that it is still slower than MATLAB if the number of neighbors is large but one can just just stick MPI into it. This at least will resolve my problem rather than using the ordering statistics software which I am sure would be useful if one understands how to use it.

John

[mecej4]

The focus of the discussion has changed substantially. For example, the first mention of k-nearest neighbors is in #5. There was a long thread on this topic three years ago in the Usenet group comp.lang.fortran .

See https://groups.google.com/forum/#!topic/comp.lang.fortran/Gad1-auDYuU .

If approximate nearest neighbors will suffice, there are some interesting possibilities; see http://www.cs.umd.edu/~mount/ANN/ .