Find the median value of an array.

Hello,

Given an array with 10,000 real values, I want to find the median. And
I want to do this about 100,000 times. My current plan is to implement
Quicksort or Merge Sort in Fortran and pick the middle point of the
array to get the median. I was just wondering Fortran already has a
way to do this. I am sure I am not the first person in history to want
the median of a group of values in Fortran. But I cannot find an
intrinsic function called "median".

Any thoughts?

Thanks.
Daniel.
0
Daniel
11/11/2010 4:49:36 PM
comp.lang.fortran 11387 articles. 0 followers. Post Follow

22 Replies
3501 Views

Similar Articles

[PageSpeed] 54
"Daniel Carrera" <dcarrera@gmail.com> wrote in message 
news:d510fdac-68cd-4637-811d-985d58ec0541@40g2000vbn.googlegroups.com...

> Given an array with 10,000 real values, I want to find the median. And
> I want to do this about 100,000 times. My current plan is to implement
> Quicksort or Merge Sort in Fortran and pick the middle point of the
> array to get the median. I was just wondering Fortran already has a
> way to do this. I am sure I am not the first person in history to want
> the median of a group of values in Fortran. But I cannot find an
> intrinsic function called "median".

I wrote up some stuff a while back:

http://home.comcast.net/~kmbtib/Fortran_stuff/order_stat.i90
http://home.comcast.net/~kmbtib/Fortran_stuff/order_stat_test.f90

Can't recall whether it's exactly what you need, but it's faster
than sorting, anyhow.  Write back if you encounter problems.

-- 
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end


0
James
11/11/2010 4:56:53 PM
In article <ibh78q$l28$1@news.eternal-september.org>,
James Van Buskirk <not_valid@comcast.net> wrote:
>"Daniel Carrera" <dcarrera@gmail.com> wrote in message 
>news:d510fdac-68cd-4637-811d-985d58ec0541@40g2000vbn.googlegroups.com...
>
>> Given an array with 10,000 real values, I want to find the median. And
>> I want to do this about 100,000 times. My current plan is to implement
>> Quicksort or Merge Sort in Fortran and pick the middle point of the
>> array to get the median. I was just wondering Fortran already has a
>> way to do this. I am sure I am not the first person in history to want
>> the median of a group of values in Fortran. But I cannot find an
>> intrinsic function called "median".
>
>I wrote up some stuff a while back:
>
>http://home.comcast.net/~kmbtib/Fortran_stuff/order_stat.i90
>http://home.comcast.net/~kmbtib/Fortran_stuff/order_stat_test.f90
>
>Can't recall whether it's exactly what you need, but it's faster
>than sorting, anyhow.  Write back if you encounter problems.

There is a much faster algorithm, incidentally.

1) Estimate a range around its value using a small, fixed sample.

2) Pass through the array, copying the values within the range,
and counting those to the left and right.

3) Repeat on the subset if the counts were good, and using the
fact that you now know which side it was if not.

4) Repeat until the sample is small enough to sort.

You can use a limited size array for the subset, incidentally,
and do a modification of the above if you run out of space.

The theory of this is standard statistics.

I have the code somewhere, but wouldn't be able to find it,
because I can't even rememember which decade I did it in and
scanning through all my archives (in half a dozen formats) is
too time-consuming.  Sorry.


Regards,
Nick Maclaren.
0
nmm1
11/11/2010 5:40:48 PM
nmm1@cam.ac.uk wrote:
(snip)
>>"Daniel Carrera" <dcarrera@gmail.com> wrote in message 
>>news:d510fdac-68cd-4637-811d-985d58ec0541@40g2000vbn.googlegroups.com...

>>> Given an array with 10,000 real values, I want to find the median. 
(snip)

> There is a much faster algorithm, incidentally.

Is that the one in Knuth, "The Art of Computer Programming", my
guess is volume 3?   It doesn't sound familiar, but I could have
forgotten how it went.
 
> 1) Estimate a range around its value using a small, fixed sample.
 
> 2) Pass through the array, copying the values within the range,
> and counting those to the left and right.
 
> 3) Repeat on the subset if the counts were good, and using the
> fact that you now know which side it was if not.
 
> 4) Repeat until the sample is small enough to sort.
 
> You can use a limited size array for the subset, incidentally,
> and do a modification of the above if you run out of space.
 
> The theory of this is standard statistics.
 
> I have the code somewhere, but wouldn't be able to find it,
> because I can't even rememember which decade I did it in and
> scanning through all my archives (in half a dozen formats) is
> too time-consuming.  Sorry.

-- glen
0
glen
11/11/2010 9:17:12 PM
In article <ibhmgo$b1q$1@news.eternal-september.org>,
glen herrmannsfeldt  <gah@ugcs.caltech.edu> wrote:
>
>>>> Given an array with 10,000 real values, I want to find the median. 
>
>> There is a much faster algorithm, incidentally.
>
>Is that the one in Knuth, "The Art of Computer Programming", my
>guess is volume 3?   It doesn't sound familiar, but I could have
>forgotten how it went.

No.  Knuth is not strong on statistical algorithms.  It includes
the one that James Van Buskirk posted.


Regards,
Nick Maclaren.
0
nmm1
11/11/2010 9:28:11 PM
nmm1@cam.ac.uk wrote:
> In article <ibhmgo$b1q$1@news.eternal-september.org>,
> glen herrmannsfeldt  <gah@ugcs.caltech.edu> wrote:

>>>>> Given an array with 10,000 real values, I want to find the median. 

>>> There is a much faster algorithm, incidentally.

>>Is that the one in Knuth, "The Art of Computer Programming", my
>>guess is volume 3?   It doesn't sound familiar, but I could have
>>forgotten how it went.
 
> No.  Knuth is not strong on statistical algorithms.  

The one in Knuth is supposed to be O(N), which doesn't mean
that it is fastest for any, or even most, sets of input data.

Some of this seems to be in the wikipedia page Selection_algorithm.

The wikipedia page for Median_graph has a reference to Knuth's
volume 4, I believe available on the web, but not in print.

> It includes the one that James Van Buskirk posted.

I don't always read posted links.

-- glen

0
glen
11/11/2010 9:43:13 PM
In article <ibho1g$b1q$3@news.eternal-september.org>,
glen herrmannsfeldt  <gah@ugcs.caltech.edu> wrote:
>
>>>>>> Given an array with 10,000 real values, I want to find the median. 
>
>>>> There is a much faster algorithm, incidentally.
>
>>>Is that the one in Knuth, "The Art of Computer Programming", my
>>>guess is volume 3?   It doesn't sound familiar, but I could have
>>>forgotten how it went.
> 
>> No.  Knuth is not strong on statistical algorithms.  
>
>The one in Knuth is supposed to be O(N), which doesn't mean
>that it is fastest for any, or even most, sets of input data.

The one I posted is also O(N) with probability one, but with a
much lower factor.  It is very common that replacing a pure
mathematical worst case with a probability zero one produces
a much faster algorithm, that is every bit as reliable in real
use.


Regards,
Nick Maclaren.
0
nmm1
11/11/2010 9:47:00 PM
Thanks for the help.

Now I know that I need to write my own function. Thanks for the
suggestions for algorithms. I'll try to pick something that I
understand, but is efficient. I have an idea for an algorithm based on
the quicksort / mergesort idea. Probably not as efficient as James'
code, but it should be good enough, and it'll be something I
understand.

Thanks for the help.

Cheers,
Daniel.


0
Daniel
11/12/2010 9:52:08 AM
On 12.11.2010. 10:52, Daniel Carrera wrote:
> Thanks for the help.
>
> Now I know that I need to write my own function. Thanks for the
> suggestions for algorithms. I'll try to pick something that I
> understand, but is efficient. I have an idea for an algorithm based on
> the quicksort / mergesort idea. Probably not as efficient as James'
> code, but it should be good enough, and it'll be something I
> understand.
>
> Thanks for the help.

Take a look at Michel Olagnon's (who used to frequent clf):

http://www.fortran-2000.com/rank/

In particular, INDMED and VALMED are based on Knuth. I used some
of other Michel's routines, but not these ones.

--
Jugoslav
www.xeffort.com
0
Jugoslav
11/12/2010 3:02:51 PM
On Nov 12, 4:02=A0pm, Jugoslav Dujic <jdu...@yahoo.com> wrote:
> Take a look at Michel Olagnon's (who used to frequent clf):
>
> http://www.fortran-2000.com/rank/
>
> In particular, INDMED and VALMED are based on Knuth. I used some
> of other Michel's routines, but not these ones.

Thanks. I'll take a look at VALMED. In fact, I have a question
already. VALMED starts with:

Recursive Function D_valmed (XDONT) Result (res_med)

I thought that function definitions always needed a return type like
"real function foo()". When I compile my own code with gfortran it
always complains if I don't include a return type, but Michel's module
compiles just fine the way it is. What's Michel doing that allows him
to get away with a function without a return type?

Anyway, I'll try to understand how Michel's code works. I did come up
with a nice, simple algorithm to get the median in O(n log n), and is
fast enough for my purposes (will add ~2.7 mins to a program that runs
8 hours). But I'll probably learn something interesting by looking at
Michel's code.

For reference, my code is below. It sorts an array of real numbers and
is inspired quicksort / mergesort:


real function qmedian(list)
    real, dimension(:), intent(in) :: list
    integer :: n
    real :: dummy

    n =3D size(list)
    dummy =3D rand(1)
    if (mod(n,2) =3D=3D 0) then
        qmedian =3D qmedian_helper(list, n/2)
    else
        qmedian =3D qmedian_helper(list, n/2 + 1)
    end if
end function

real recursive function qmedian_helper(list, q) result(res)
    ! Parameters
    real, dimension(:), intent(in) :: list
    integer, intent(in) :: q

    ! Internal variables
    real, allocatable, dimension(:) :: left, right
    real :: pivot, lmax, lmin
    integer :: n, j

    ! Algorithm
    n =3D size(list)
    call bounds(list, lmin, lmax)
    if (q =3D=3D 1) then
        res =3D lmin
    else if (q =3D=3D n) then
        res =3D lmax
    else
        !!! Split 'list' into two parts and recurse.

        ! Pick a random element from the list.
        pivot =3D list(int(rand(0)*n)+1)

        j =3D count(list <=3D pivot)

        ! The previous conditions should remove corner cases.
        allocate(left(j))
        allocate(right(n-j))

        left =3D pack(list, list <=3D pivot)
        right=3D pack(list, list  > pivot)

        if (j >=3D q) then
            ! Left side contains kth element.
            res =3D qmedian_helper(left, q)
        else
            ! Right side contains kth element.
            res =3D qmedian_helper(right, q-j)
        end if
    end if
end function

subroutine bounds(list, lmin, lmax)
    real, dimension(:), intent(in) :: list
    real, intent(out) :: lmin, lmax
    integer :: n, i

    n =3D size(list)
    lmin =3D list(1)
    lmax =3D list(1)
    do i =3D 2,n
        if (list(i) < lmin) lmin =3D list(i)
        if (list(i) > lmax) lmax =3D list(i)
    end do
end subroutine


Cheers,
Daniel.
0
Daniel
11/12/2010 7:47:41 PM
Daniel Carrera wrote:
> On Nov 12, 4:02 pm, Jugoslav Dujic <jdu...@yahoo.com> wrote:
>> Take a look at Michel Olagnon's (who used to frequent clf):
>>
>> http://www.fortran-2000.com/rank/
>>
>> In particular, INDMED and VALMED are based on Knuth. I used some
>> of other Michel's routines, but not these ones.
> 
> Thanks. I'll take a look at VALMED. In fact, I have a question
> already. VALMED starts with:
> 
> Recursive Function D_valmed (XDONT) Result (res_med)
> 
> I thought that function definitions always needed a return type like
> "real function foo()". 

I suspect the answer is made clear if you look at the first point under
"programming style" in the above link,

Ian
0
Ian
11/12/2010 7:57:41 PM
On Nov 12, 11:47=A0am, Daniel Carrera <dcarr...@gmail.com> wrote:
> On Nov 12, 4:02=A0pm, Jugoslav Dujic <jdu...@yahoo.com> wrote:
>
> > Take a look at Michel Olagnon's (who used to frequent clf):
>
> >http://www.fortran-2000.com/rank/
>
> > In particular, INDMED and VALMED are based on Knuth. I used some
> > of other Michel's routines, but not these ones.
>
> Thanks. I'll take a look at VALMED. In fact, I have a question
> already. VALMED starts with:
>
> Recursive Function D_valmed (XDONT) Result (res_med)
>
> I thought that function definitions always needed a return type like
> "real function foo()". When I compile my own code with gfortran it
> always complains if I don't include a return type, but Michel's module
> compiles just fine the way it is. What's Michel doing that allows him
> to get away with a function without a return type?

Removing the large comment from Michel's code, I finds
the first 3 lines:

Recursive Function D_valmed (XDONT) Result (res_med)
      Real (kind=3Dkdp), Dimension (:), Intent (In) :: XDONT
      Real (kind=3Dkdp) :: res_med

The return type is declared in line 3.

--
steve
0
steve
11/12/2010 7:59:39 PM
Daniel Carrera <dcarrera@gmail.com> wrote:

> Recursive Function D_valmed (XDONT) Result (res_med)
> 
> I thought that function definitions always needed a return type like
> "real function foo()". When I compile my own code with gfortran it
> always complains if I don't include a return type, but Michel's module
> compiles just fine the way it is. What's Michel doing that allows him
> to get away with a function without a return type?

He presumably declares the type later in the code in a separate
statement. Not only is that allowed, it is what I prefer and recommend.
There can be problems with declaring the type directly on the function
statement; it works in most cases, but not all of them. Declaring the
type in a separate type declaration statement works in *ALL* cases. I
prefer to do things that work all the time so that I don't need to worry
about checking when I have crossed the boundary between the cases that
work and the ones that don't.

Also note that when you declare the result variable to have a separate
name, it is the result variable whose type gets declared - not the
function.

As an aside, no, you don't have to explicitly declare the type. If you
don't declare the type, implicit typing applies, just as with any other
variable. I disrecommend the use of implicit typing, and if you use
impkicit none, the compiler will bitch about it, but it is allowed. If
your compiler complains when you omit the type, then probably you are
using implicit none (which I do recommend). Either that or you are doing
something that is incompatible with the implicit type. Be aware that
implicit typing does *NOT* mean that an entity has no type; it just
means that the type is determined by the implicit typing rules instead
of by an explicit type declaration. I wouldn't mention this except that
I wonder why you would ever be trying to write code that doesn't declare
the return type and also isn't deliberately using implicit typing; the
only reason I can think of (other than just an accident) is that one
might think that failing to explicitly declare the type somehow left it
generic.

-- 
Richard Maine                    | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle           |  -- Mark Twain
0
nospam
11/12/2010 8:05:11 PM
On Nov 12, 8:59=A0pm, steve <kar...@comcast.net> wrote:
> On Nov 12, 11:47=A0am, Daniel Carrera <dcarr...@gmail.com> wrote:
>
>
>
>
>
> > On Nov 12, 4:02=A0pm, Jugoslav Dujic <jdu...@yahoo.com> wrote:
>
> > > Take a look at Michel Olagnon's (who used to frequent clf):
>
> > >http://www.fortran-2000.com/rank/
>
> > > In particular, INDMED and VALMED are based on Knuth. I used some
> > > of other Michel's routines, but not these ones.
>
> > Thanks. I'll take a look at VALMED. In fact, I have a question
> > already. VALMED starts with:
>
> > Recursive Function D_valmed (XDONT) Result (res_med)
>
> > I thought that function definitions always needed a return type like
> > "real function foo()". When I compile my own code with gfortran it
> > always complains if I don't include a return type, but Michel's module
> > compiles just fine the way it is. What's Michel doing that allows him
> > to get away with a function without a return type?
>
> Removing the large comment from Michel's code, I finds
> the first 3 lines:
>
> Recursive Function D_valmed (XDONT) Result (res_med)
> =A0 =A0 =A0 Real (kind=3Dkdp), Dimension (:), Intent (In) :: XDONT
> =A0 =A0 =A0 Real (kind=3Dkdp) :: res_med
>
> The return type is declared in line 3.

Interesting. I didn't know you could do that. I was beginning to think
that he was relying on implicit types in some way. I notice that he
doesn't use "IMPLICIT NONE".

Thanks. I learned something new today.

Daniel.
0
Daniel
11/12/2010 8:16:08 PM
On Nov 12, 9:05=A0pm, nos...@see.signature (Richard Maine) wrote:
> He presumably declares the type later in the code in a separate
> statement. Not only is that allowed, it is what I prefer and recommend.

Yes, he does. I didn't realize that was possible until Steve pointed
it out.


> There can be problems with declaring the type directly on the function
> statement; it works in most cases, but not all of them. Declaring the
> type in a separate type declaration statement works in *ALL* cases. I
> prefer to do things that work all the time so that I don't need to worry
> about checking when I have crossed the boundary between the cases that
> work and the ones that don't.

Sounds reasonable. I'm surprised that declaring the type on the
function statement can fail. Can you give me an example? Anyway, I
certainly want to use methods that work all the time.


> If your compiler complains when you omit the type, then probably you are
> using implicit none (which I do recommend). Either that or you are doing
> something that is incompatible with the implicit type.

Yes, I always use implicit none. I only have a few weeks knowledge of
Fortran 95, but I learned to use implicit none from the very
beginning. Implicit typing scares me. I'm surprised that Michael
doesn't use implicit none.


> I wouldn't mention this except that
> I wonder why you would ever be trying to write code that doesn't declare
> the return type and also isn't deliberately using implicit typing; the
> only reason I can think of (other than just an accident) is that one
> might think that failing to explicitly declare the type somehow left it
> generic.


It was just an accident. I just forgot to enter the type once and saw
the compiler error. I just "understood" that every function needs a
return type, and I was surprised when I didn't see one in Michael's
code. I never thought of implicit typing (a concept that still seems
very odd to me).

Cheers,
Daniel.
0
Daniel
11/12/2010 8:30:35 PM
Daniel Carrera <dcarrera@gmail.com> wrote:
> On Nov 12, 9:05�pm, nos...@see.signature (Richard Maine) wrote:
>> He presumably declares the type later in the code in a separate
>> statement. Not only is that allowed, it is what I prefer and recommend.
 
> Yes, he does. I didn't realize that was possible until Steve pointed
> it out.
 
>> There can be problems with declaring the type directly on the function
>> statement; it works in most cases, but not all of them. Declaring the
>> type in a separate type declaration statement works in *ALL* cases. I
>> prefer to do things that work all the time so that I don't need to worry
>> about checking when I have crossed the boundary between the cases that
>> work and the ones that don't.
 
> Sounds reasonable. I'm surprised that declaring the type on the
> function statement can fail. Can you give me an example? Anyway, I
> certainly want to use methods that work all the time.

The separate return variable is required for direct recursion.
(When a function calls itself directly.)  Inside a function without
a separate return variable, either declared on the FUNCTION statement
or in a separate statement, the function name is a variable, and
not a function reference.

-- glen
0
glen
11/12/2010 8:43:13 PM
Daniel Carrera <dcarrera@gmail.com> wrote:

> On Nov 12, 9:05 pm, nos...@see.signature (Richard Maine) wrote:

> > There can be problems with declaring the type directly on the function
> > statement; it works in most cases, but not all of them. Declaring the
> > type in a separate type declaration statement works in *ALL* cases. I
> > prefer to do things that work all the time so that I don't need to worry
> > about checking when I have crossed the boundary between the cases that
> > work and the ones that don't.
> 
> Sounds reasonable. I'm surprised that declaring the type on the
> function statement can fail. Can you give me an example? Anyway, I
> certainly want to use methods that work all the time.

Looking at the code in question, I see that it isn't actualy too far
from the main case that I think fails. Basically, the problems come from
ordering requirements. The function statement is always the first
statement in the function; there are no exceptions to that.

There are some things that you can't use until after the statement that
defines them. If you want to use any of those things on the function
statement, there can be a problem in that the function statement is not
likely to be after the statement where they are defined.

A derived type is one of those things. Normally you can't use a derived
type until after its definition. I think (not worth checking, but I
think it is right) that a special-case exception was made for using a
derived type in a function statement just in order to allow it. I don't
tend to like such special cases and I'd have probably voted against
making one for this (don't recall whether I ever got a chance to make
such a vote - might have). Special cases tend to grow as one discovers
that you didn't catch quite everything. Indeed, this special case didn't
catch everything.

I think the missing case is type parameters. One often (and preferably)
specifies type parameter values using named constants. I do no think
there is a simillar exception to allow using a named constant in a
function statement before its definition. So if you try to do something
like

  real(my_kind) function f(x)
    use module_with_a_definition_of_my_kind
    ...

I think you are out of luck. It is possible that I've gotten confused on
the details and that this particular case works. But if it isn't exactly
this one, it is something simillar. One of the advatanges of doing the
type declararation in a separate statement is that I don't actually have
to keep track of the fine points about exactly which cases work.

Posted excerpts of the posted code do show that it uses a named constant
for a kind type parameter, so it is at least close to this case. I don't
see a USE in the posted extract though, so perhaps the kind type
parameter might be accessed via host association. That would be a case
that was ok.... well as long as it wasn't an interface body.

-- 
Richard Maine                    | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle           |  -- Mark Twain
0
nospam
11/12/2010 9:21:31 PM
On Nov 12, 10:21=A0pm, nos...@see.signature (Richard Maine) wrote:
> A derived type is one of those things. Normally you can't use a derived
> type until after its definition...
>
> I think the missing case is type parameters. One often (and preferably)
> specifies type parameter values using named constants. I do no think
> there is a simillar exception to allow using a named constant in a
> function statement before its definition. So if you try to do something
> like
>
> =A0 real(my_kind) function f(x)
> =A0 =A0 use module_with_a_definition_of_my_kind
> =A0 =A0 ...
>
> I think you are out of luck.


I just found another example. Suppose you want a function that returns
an array:

    function myfunc() result(res)
        real, dimension(3) :: res
        res =3D (/ 2.3, 2.4, 2.5 /)
    end function

This compiles an runs correctly. But when I tried move the type
declaration to the function definition, the compilation failed.

Cheers,
Daniel.
0
Daniel
11/12/2010 10:38:25 PM
Daniel Carrera <dcarrera@gmail.com> wrote:

> I just found another example. Suppose you want a function that returns
> an array:
> 
>     function myfunc() result(res)
>         real, dimension(3) :: res
>         res = (/ 2.3, 2.4, 2.5 /)
>     end function
> 
> This compiles an runs correctly. But when I tried move the type
> declaration to the function definition, the compilation failed.

Good point, thought I'll make a minor terminology correction. You
actually can still put the type into the function statement. You just
can't put the dimension attribute there. The dimension is not part of
the type even though you normally put it in the type declaration
statement; it is an attribute.

You can do (and I just checked to make sure that I hadn't missed
something - worked fine).

    real function myfunc() result(res)
        dimension :: res(3)
        res = (/ 2.3, 2.4, 2.5 /)
    end function

I do not recommend that at all. I personally find it abominable. But it
is standard conforming.

Now that you mention attributes, that reminds me that things like the
pointer attribute have simillar issues in not being specifiable on the
function statement. But I don't even want to mention that it is a
possibility for a function to return a pointer because that might prompt
someone to write such functions, which I advise against. Pretend I
didn't say anything. :-)

Functions that return allocatables are a much better example in that I
think they can be good ideas.

-- 
Richard Maine                    | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle           |  -- Mark Twain
0
nospam
11/13/2010 12:54:35 AM
Hi, I didn't read everything above.
But a search did not find the term "partial".
IMHO, quantiles, median belongs to it, are fastest found with "partial
sorting".
"Partial sorting" means that the data are not fully sorted, it is only
sorted so far that all smaller values are before the median and all
larger
values are behind, but they are not exactly sorted.
There are good algorithms available if you just look for
"partial sorting".
Wolfgang

On Nov 12, 7:54=A0pm, nos...@see.signature (Richard Maine) wrote:
> Daniel Carrera <dcarr...@gmail.com> wrote:
> > I just found another example. Suppose you want a function that returns
> > an array:
>
> > =A0 =A0 function myfunc() result(res)
> > =A0 =A0 =A0 =A0 real, dimension(3) :: res
> > =A0 =A0 =A0 =A0 res =3D (/ 2.3, 2.4, 2.5 /)
> > =A0 =A0 end function
>
> > This compiles an runs correctly. But when I tried move the type
> > declaration to the function definition, the compilation failed.
>
> Good point, thought I'll make a minor terminology correction. You
> actually can still put the type into the function statement. You just
> can't put the dimension attribute there. The dimension is not part of
> the type even though you normally put it in the type declaration
> statement; it is an attribute.
>
> You can do (and I just checked to make sure that I hadn't missed
> something - worked fine).
>
> =A0 =A0 real function myfunc() result(res)
> =A0 =A0 =A0 =A0 dimension :: res(3)
> =A0 =A0 =A0 =A0 res =3D (/ 2.3, 2.4, 2.5 /)
> =A0 =A0 end function
>
> I do not recommend that at all. I personally find it abominable. But it
> is standard conforming.
>
> Now that you mention attributes, that reminds me that things like the
> pointer attribute have simillar issues in not being specifiable on the
> function statement. But I don't even want to mention that it is a
> possibility for a function to return a pointer because that might prompt
> someone to write such functions, which I advise against. Pretend I
> didn't say anything. :-)
>
> Functions that return allocatables are a much better example in that I
> think they can be good ideas.
>
> --
> Richard Maine =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0| Good judgment come=
s from experience;
> email: last name at domain . net | experience comes from bad judgment.
> domain: summertriangle =A0 =A0 =A0 =A0 =A0 | =A0-- Mark Twain

0
wolfgang
11/13/2010 9:15:56 AM
On 11/12/2010 09:30 PM, Daniel Carrera wrote:
> On Nov 12, 9:05 pm, nos...@see.signature (Richard Maine) wrote:
>> He presumably declares the type later in the code in a separate
>> statement. Not only is that allowed, it is what I prefer and recommend.
>
> Yes, he does. I didn't realize that was possible until Steve pointed
> it out.
>
>
>> There can be problems with declaring the type directly on the function
>> statement; it works in most cases, but not all of them. Declaring the
>> type in a separate type declaration statement works in *ALL* cases. I
>> prefer to do things that work all the time so that I don't need to worry
>> about checking when I have crossed the boundary between the cases that
>> work and the ones that don't.
>
> Sounds reasonable. I'm surprised that declaring the type on the
> function statement can fail. Can you give me an example? Anyway, I
> certainly want to use methods that work all the time.
>
>
>> If your compiler complains when you omit the type, then probably you are
>> using implicit none (which I do recommend). Either that or you are doing
>> something that is incompatible with the implicit type.
>
> Yes, I always use implicit none. I only have a few weeks knowledge of
> Fortran 95, but I learned to use implicit none from the very
> beginning. Implicit typing scares me. I'm surprised that Michael
> doesn't use implicit none.
>

Well, I write a lot of quick and dirty programs, but when it is code 
that I want to keep or publish, I do check for undeclared variables.
In the case of the sorting and ranking routines, I compiled them with an 
F compiler, that enforces implicit none, so everything is declared.

Regarding the way declarations are arranged, some features come from the 
fact that I used a pre-processor to build the various type versions of 
the routines, by lack of a suitable mechanism in the language (at least 
at that time, I don't know if something has appeared now).

The code looked at first like:
Module m_mrgrnk
Integer, Parameter :: kdp = selected_real_kind(15)
public :: mrgrnk
private :: kdp
private :: R_mrgrnk, I_mrgrnk, D_mrgrnk
interface mrgrnk
   module procedure D_mrgrnk, R_mrgrnk, I_mrgrnk
end interface mrgrnk
contains

Subroutine D_mrgrnk (XDONT, IRNGT)
$define MYTYPE "Real (kind=kdp)"
$include "mrgrnk.inc"
End Subroutine D_mrgrnk

Subroutine R_mrgrnk (XDONT, IRNGT)
$define MYTYPE "Real"
$include "mrgrnk.inc"
End Subroutine R_mrgrnk

Subroutine I_mrgrnk (XDONT, IRNGT)
$define MYTYPE "Integer"
$include "mrgrnk.inc"
End Subroutine I_mrgrnk
End Module m_mrgrnk



0
Michel
11/15/2010 8:20:52 AM
On 11/12/2010 04:02 PM, Jugoslav Dujic wrote:
>
> On 12.11.2010. 10:52, Daniel Carrera wrote:
>> Thanks for the help.
>>
>> Now I know that I need to write my own function. Thanks for the
>> suggestions for algorithms. I'll try to pick something that I
>> understand, but is efficient. I have an idea for an algorithm based on
>> the quicksort / mergesort idea. Probably not as efficient as James'
>> code, but it should be good enough, and it'll be something I
>> understand.
>>
>> Thanks for the help.
>
> Take a look at Michel Olagnon's (who used to frequent clf):
>
> http://www.fortran-2000.com/rank/
>
> In particular, INDMED and VALMED are based on Knuth. I used some
> of other Michel's routines, but not these ones.
>

I programmed Knuth's algorithm for the sake of sport, but in practice it 
is so complex that you would have to go to huge values of N to see some 
advantage over the simple pivoting strategy.
0
Michel
11/15/2010 8:27:08 AM
In article <8kc96tFfv4U1@mid.individual.net>,
Michel Olagnon  <molagnon@ifremer-a-oter.fr> wrote:
>
>I programmed Knuth's algorithm for the sake of sport, but in practice it 
>is so complex that you would have to go to huge values of N to see some 
>advantage over the simple pivoting strategy.

That is true, but the statistical one is enough simpler and faster
that you start to see gains at quite low values - if I recall, it's
about 20-25.


Regards,
Nick Maclaren.
0
nmm1
11/15/2010 8:49:44 AM
Reply:
Similar Artilces:

Rearranging arrays.
Hi, thanks in advance for your input. I'm a super newbie Matlab user, so also thanks for your patience! :) I have a collection of data files: "z.txt" that are x,y pairs. The x column in each data file is the same. I have written a program that will choose a particular x value. I need to find the y value for that x in each z.txt file and make a new array that is y vs. z. Then I want to be able to do this with many x values and save all the y vs. z arrays for each x to the same file. Ex: 100.txt: 0.10,0.44 0.20,0.32 0.30,0.68 0.40,0.84 0.50, 0.24 0.60, 0.4...

Find and Grep together
Hi All, I am sorry if this is not the right group for this question. I am using find and grep together to find a string in a file from the top level directories. find . -name "*.c" -exec grep "string" {} \; -print | more. | \/ top level directorie I am sure that the string I am searching is in the lower directories where I am searching but this command is unable to find the string which I am looking. Can anyone please help me. My version control tool is clear case and my config spec also supports all the files....

Multiply 3D array by 2D array
Sorry, brain tired... can someone remind me what the fastest way to do this is: a = rand(m,m,n); b = rand(m,m); for j=1:n a(:,:,j) = a(:,:,j) .* b; end spasmous: <SNIP <nightmare on elmstreet>, featuring robert englund and spasmous... one of the many solutions m=4; n=3; a=rand(m,m,n); b=rand(m,m); c=a.*repmat(b,[1 1 n]) us us wrote: > spasmous: > <SNIP <nightmare on elmstreet>, featuring robert englund and > spasmous... > > one of the many solutions > > m=4; > n=3; > a=rand(m,m,n); > b=rand...

How to find out the current version of the MySQL installation ?
How to find out (from the command line under Linux) the current version of the MySQL installation ? Or is there a config file which contains the version number ? Werner Werner Sammer <wersam@yahoo.de> wrote: > How to find out (from the command line under Linux) > the current version of the MySQL installation ? > Or is there a config file which contains the version number ? > Werner Command line: mysql --version Regards, Johan -- _____________________________________ Ing. Johan van Oostrum chaos geordend - www.chaosgeordend.nl _____________________________________ >...

finding max
Hello all I have a matrix a = 1 2 3 4 5 6 7 8 9 10 11 12 now I call max(a) ans = 9 10 11 12 to find the max, I have to call max(max(a)) Is there any way I can get the max of a matrix in one command? Thanks sore wrote: > > Hello all > > I have a matrix > > a = > > 1 2 3 4 > 5 6 7 8 > 9 10 11 12 > > now I call > > max(a) > > ans = > > 9 10 11 12 > > to find the max, I have to call max(max(a)) > > Is there any way I can get the max of a matrix > in one command? >...

Find it all
Find it all http://asb-comm-logicx.net/ ...

Narrowed down bottlneck to disk, how can I find out which FS is being hammered? 100% utilization
In my process of drilling down to get to the bottom of our performance issues, I am seeing huge I/O wait on a particular disk. What tools are available to dig deeper, find out which file system is hammering the disk and/or which processes are hogging the CPUS? We are talking 100% utilization here: mpstat output: CPU minf mjf xcal intr ithr csw icsw migr smtx srw syscl usr sys wt idl 0 199 0 259 12 2 176 9 22 30 0 1882 11 6 83 0 1 205 0 259 204 194 176 8 23 34 0 1801 11 5 84 0 2 203 0 289 34 23 176 7 23 31 ...

DigitalWrite sequence, Filtering digital data arrays
Hi all, I�ve been trying to solve the following issue the last couple of days with no success, so I hope one of you can help me out.The first problem is that I�ve changing waveform data and sometimes there are two distinct digital data at one time point. As shown in the brief example attached, for timepoint 7500 the input is either '1000' or '0000' (corresponding to row 3 and 4). In this particular case I�d like to delete the first entry (row 3), since the 'one' pulse of the first entry sends out a short high-pulse although it should just contain zeros. Unfortunately, t...

Complicated Find Expression
Hi, I am having problems with my find expression. Basically I want to find to find all files that are _not_ in a test directory and then execute a grep statement on them. The result of that grep statement is further refined by another grep. Finally I would like to echo the result. My grep expression on works fine by itself. The problem is with the find integration. Here is my expression, followed by the output. Thanks for any suggestions/pointers, Mike find ./ -path '*test*' -prune -o -print -exec "grep -PRi 'out.print.*?[^/][^/].*?the' {} | grep -Pi 'out...

[ace-users] Can't find ACEd.dll
This is a multi-part message in MIME format. ------_=_NextPart_001_01C7A8F9.1452B078 Content-Type: text/plain; charset="us-ascii" Content-Transfer-Encoding: quoted-printable Hi=20 =20 ACE VERSION: 5.5.8 =20 HOST MACHINE and OPERATING SYSTEM: Win XP=20 If on Windows based OS's, which version of WINSOCK do you use?: don't know =20 TARGET MACHINE and OPERATING SYSTEM, if different from HOST: COMPILER NAME AND VERSION (AND PATCHLEVEL): =20 THE $ACE_ROOT/ace/config.h FILE [if you use a link to a platform- specific f...

ow-find without sun 'find' key
I've just moved over to a windows platform from solaris and have installed the windows version of xemacs 21.4. I'm having severe withdrawal symptoms from the loss of the sun 'find' key. So to correct this I would like to bind the M-s key combo to the ow-find function and M-S-s to the ow-find-backward function. I have 2 problems trying to achieve this. 1. I have sucessfully bound the M-s keys to the 'ow-find function using the following code in the init.el file. (define-key global-map [(meta ?s)] 'ow-find) However when I execute the function it returns the follo...

Find File Attributes
Hello, As a part of my application,I need to compare files generated everyday.I need not need to compare the file contents.I just need to compare the file size,the line count.The files are named in the format,Hyyyymmdd,H denoting History. At present,I use the Dir function to search for the physical existence of the file and I use the debug.print for functions such as filelen(path),GetAttr() etc. My question is how would I put this as a module so that I can use a date function to compare files generated today with that of yesterday (mondays files to be compared against fridays). Also is it pos...

SVM feature and value pair
Hi, I think I am a new comer here. I am trying to use svm_light software. From the example, I can see that the training data is in pairs. feature and value related to it. I have a question about the value, how is the value computed? Is there a special system to calculate the value ? In other words, how can I find the value associated with the feature? I will be very appreciated if someone has experiences on this. Thanks zhenyu [ comp.ai is moderated. To submit, just post and be patient, or if ] [ that fails mail your article to <comp-ai@moderators.isc.org>, and ] [ ask your news adm...

Table of contents, {rl} in array speak
Because playing with formatting is more fun than... say writing content for the thesis, I've been trying to come up with a way to make the TOC come out as suggested in Elements of Typographical style. instead of 1. Foo ............................... 3 1.2 Blarg ............................ 10 2. Grok .............................. 11 (which seems to be all you can play with tocloft) I'd like: 3 -- 1. Foo 10 -- 1.2 Blarg 11 -- 2 Grok or some variation thereof. I've tried redefining \l@X, f.ex: \renewcommand*\l@chapter[2]{% #2 \hbox{#1} \par } Th...

I need to find the rows that exist in one table but not in the other with condition
I need to find the rows that exist in one table but not in the other with this condition: (prod_name exist in table1 and not in table2.prod_name ) AND (prod_name exist in table1 and not in table2.'S'+prod_name ) explanation: i want to know if the product not exit and if the combination of the charachter "S" with the product Name also not exist at the other table B.R yuvi SELECT prod_name FROM table1 as A WHERE NOT EXISTS (select * from table2 as B where A.prod_name = B.prod_name) AND NOT EXISTS (select * from table2 as C where A....

find breakpoint programmatically
Hi All, &nbsp; I accidentally posted this question in the special interest forum by mistake so sorry for the duplicate post but i thought it might never get read over there. <a href="http://forums.ni.com/ni/board/message?board.id=BreakPoint&amp;message.id=4162" target="_blank">http://forums.ni.com/ni/board/message?board.id=BreakPoint&amp;message.id=4162</a> &nbsp; Hi all, &nbsp; I was wondering if there is a way to search the vi hierachy for breakpoints programmatically?&nbsp; Why you may ask? &nbsp; I have a piece of code that take...

Newbie advice and installing on a RAID-0 Array
Built a new system over the holidays, and I figure why not dedicate a partition to Linux so I can start learning it and hopefully find something to use it for. I have never messed with it before and have basically zero knowledge of the OS itself. I have no doubt I'll be able to figure out whatever installation directions there are, but the problem lies with selecting which distribution of Linux to use. I'm pretty quick on picking up things once I get into it so I don't want something completely dumbed down, but just something where I'll be able to operate simple progr...

counting repetitions of values in a vector
Hi, I need to write a function that takes a row vector as input, and counts the number of repetitions of each value in the vector, returning a structure with each value and the corresponding number of repetitions. For example: A = [1 1 1 2 2 2 2 3 3 3 3 3]; B = my_function(A) val = [1 2 3] rep = [3 4 5] I am stuck on how to count the repetitions of each value. Advice please? NOTE: This is for a piece of coursework I have to complete. Therefore, I do NOT want a code solution, I just want some advice as to how to tackle this problem. Thanks in advance. rupture wrote: > Hi, > >...

Re: Find unique records and recode
data one; input old $; cards; 5269-7c 5269-8c 5269-8c 5269-3z 5222-7z 5222-7z 5222-7z ; data two; set one; by old notsorted; new + first.old; run; On Mon, 5 Sep 2005 07:24:43 -0700, tortoise <cychen9@GMAIL.COM> wrote: >Hi, >I have one character variable (old). >There are many combinations of "old" variable. >How could I find unique and recode those data to "new"????? > > old new >5269-7c 1 >5269-8c 2 >5269-8c 2 >5269-3z 3 >5222-7z 4 >5222-7z 4 >5222-7z ...

How to find out memory leak?
Hi, all! I'm having memory leak problem with my program, but can't find where's the leak. I have a program written in pure-Ruby. It implemented a Chinese word segment algorithm. I call `segment(text)' to get the result. But each time I call `segment', the memory usage of my program (Ruby) increased several handred K- bytes -- the number of bytes of increment is roughly equal during each call. I can see from outside (using `ps' or `top') that the memory is leaking. But I can't find where the leaking goes. I tried to use ObjectSpace.each_object to find out what ...

pro/piping tutorials
Hello, does anyone know where I could find tutorials about pro/E module: pro/ piping? I would really appreciate any info on this. Thanks a lot! This is a multi-part message in MIME format. ------=_NextPart_000_0033_01C79EE3.B10443C0 Content-Type: text/plain; charset="iso-8859-2" Content-Transfer-Encoding: quoted-printable <dugave_11@yahoo.com> wrote in message = news:1180099763.016860.34260@o5g2000hsb.googlegroups.com... Hello, does anyone know where I could find tutorials about pro/E module: pro/ piping? I would really appreciate any info on ...

Help finding AIO w/Fax
I looking for an All In One with FAX (flatbed scanner) for less than $200. I have only been able to fine these two: Canon MP390 and HP Officejet 5510. I do not see any that are Epson, Dell, Ricoh or Samsung. I have ruled out Visioneer and Lexmark because of bad reviews. Does anyone know of any other they recommend. Thanks in advance for your help, Bill ...

Numeric Control does not generate Value Changed Event on Return
Hi all, &nbsp; I seem to have changed something on my LabVeiw 8.5.1 front panel which has stopped all of my Numeric Controls generating Value Changed Events when the Return or Enter key is pressed.&nbsp; The Value Change Event only fires if you TAB from the control or mouse click somewhere else. &nbsp; Can someone please tell me how to get the normal functionality back! &nbsp; Thanks &nbsp; P Did you have specific key navigation configured for the control(s)? Hi smercurio, No key navigation&nbsp;defined for any controls P &nbsp; you need to define a key navigati...

help
I received an Excel data file that contains a "+/-" symbol (html code \&plu= smn; &plusmn;), that can be copied and displayed in Word, Notepad, "Kompoze= r" html editor, unix vi, pico editors, and load to/retrieve from MySQL oper= ated on linux. But when I need to manipulate the data in perl, I am lost a= s how to recognize the symbol with RE. Could anyone help? Thanks in advance! joe On 2012-12-15 07:19, Joe wrote: > I received an Excel data file that contains a "+/-" symbol (html code \&plusmn; &plusmn;), that can be copied and...

Array of pointers in a struct
Hi all, I feel unclear about what my code is doing, although it works but I am not sure if there is any possible bug, please help me to verify it. This is a trie node (just similar to tree nodes) struct, I am storing an array of 27 pointers and a void pointer that can point to anything. typedef struct trieNode { struct trieNode *children[27]; // The children nodes void *obj; // The object stored } TrieNode; And this is the code I make a new trie node, initialize and return it. TrieNode *newTrieNode() { int i; // Allocate memory for the node TrieNode *node = (TrieNode *)malloc(sizeof(...