Blitz logo

Blitz Support :

From: daniel.egloff_at_[hidden]
Date: 2005-06-02 08:19:57


Dear Frank and Julian

I found a dimension independent way to implement a partial sort along
every dimension of a blitz++ array.
I basically use an blitz++ stl conformant iterator to walk through a
top region of the array.
This approach might also be usefull to extend the blitz reduction
functionality to any dimension.
The sorting is performed using the stl::make_heap algorithm. To make
it working I use a strided iterator
which is basically the one provided in the matrix template library
MTL.

As the solution might be interesting for other blitz++ users I post it
upstreams. Let me know what you
are thinking of the code snippet (ev. I still made some mistakes, its
not thoroughly tested) and if you
want to include it in some form in the blitz library.

template<typename T, int rank_N>
void partial_sort(Array<T, rank_N> &array, int sortDim, int num)
{
    TinyVector<int, rank_N> lower(0);
    TinyVector<int, rank_N> upper = array.ubound();
    upper[sortDim] = 0;

    RectDomain<rank_N> domain(lower, upper);
    Array<T, rank_N> top = array(domain);

    typename Array<T, rank_N>::iterator it = top.begin();
    typename Array<T, rank_N>::iterator et = top.end();

    int stride = array.stride(sortDim);
    for(; it != et; ++it)
    {
        strided_iterator<T*> start(&(*it), stride, 0);
        strided_iterator<T*> end = start + array.extent(sortDim);

        make_heap(start, end);
        for (typename strided_iterator<T*>::difference_type j = 0; j
< num; j++)
        {
            pop_heap(start, end - j);
        }
    }
}

(See attached file: blitzsort.cpp)

Freundliche Grüsse
Daniel Egloff
Zürcher Kantonalbank, VFEF
Josefstrasse 222, 8005 Zürich
Tel. +41 (0) 44 292 45 33, Fax +41 (0) 44 292 45 93
Briefadresse: Postfach, 8010 Zürich, http://www.zkb.ch

|---------+------------------------------->
| | Frank Schimmel |
| | <f.schimmel_at_fluid.tu|
| | wien.ac.at> |
| | |
| | 31.05.2005 12:59 |
| | Bitte antworten an |
| | Frank Schimmel |
| | |
|---------+------------------------------->
>----------------------------------------------------------------------------------------------------------------------------|
  | |
  | An: daniel.egloff_at_[hidden], blitz-support_at_[hidden] |
  | Kopie: |
  | Thema: Re: [Blitz-support] Search for dimension independent code |
>----------------------------------------------------------------------------------------------------------------------------|

>>>>> daniel egloff writes:

Hi Daniel!

> how can I make the following loop dimension independent, e.g. the
> dimension 4 should be variable, but compile time fixed.
> [...]

The number of dimensions of a blitz array must be a compile-time
constant. If that constant varies from compilation to compilation,
writing nested loops is not possible. (And besides: they don't perform
very well.) You could try to use an iterator instead (see below).

If that does not work, you can wrap you iteration in a function and
then you can have different versions of that function acting on you
arrays of different dimensionality using partial template
specialisation:

template<int dimension> class Helper;

template<>
class Helper<4> {
public:
  static const int dimension = 4;
  doStuff(blitz::Array<double,4> &array) {
    ...
  }
};

I'm not sure the syntax is all correct, though.
The kludge with the class is necessary, because partial template
specialisation exists for member functions only. :-/

You could also try to avoid the kludgy part using overloading to
emulate partial specialisation. That requires some deeper thinking,
tough. I'd reccomend the book "Modern C++ Design" by Andrei
Alexandrescu (http://www.moderncppdesign.com/book/main.html), which is
where I learned avout that trick. It's all implemented in the Loki C++
Library: http://sourceforge.net/projects/loki-lib/

> a related question: how can I make the following loop dimension
> independent:
> [code ...]
> // apply a function to every array element,
> [more code ...]

Have a look at the Array::iterator (typedef for ArrayIterator
instantiation).

for(blitz::Array::iterator i = array.begin(), e = array.end();
    i != e; ++i) {
  f(i.position(), *i);
}

The position() method gives the index of the element as a TinyVector.

> I think the second problem is related to functors. Whats about the
> first?
> Any help greatly appreciated.

Hope that does help
-Frank

PS: I hope your actual code is indented properly.
PPS: Please suggest to your legal department to shorten that
disclaimer.

___________________________________________________________________

Disclaimer:

Diese Mitteilung ist nur fuer die Empfaengerin / den Empfaenger
bestimmt.

Fuer den Fall, dass sie von nichtberechtigten Personen empfangen wird,
bitten wir diese hoeflich, die Mitteilung an die ZKB zurueckzusenden
und anschliessend die Mitteilung mit allen Anhaengen sowie allfaellige
Kopien zu vernichten bzw. zu loeschen. Der Gebrauch der Information
ist verboten.

This message is intended only for the named recipient and may contain
confidential or privileged information.

If you have received it in error, please advise the sender by return
e-mail and delete this message and any attachments. Any unauthorised
use or dissemination of this information is strictly prohibited.