I am looking for a C++ library that implements handling of 2D (or nD)
matrixes which are static at runtime (but not at compile time).
My problem is that with standart C or C++ array definitions (i.e:
<type> **A), the access to one element of A via A[x][y] is too slow.
As at runtime the size of the matrix is known, all the extension features
offered by the <type> **A declaration are not needed.
I thought about library reserving memory for the matrix as if it was a
unidimensional vector of size #lines * #columns, but letting me access
to the elements with the standard 2D array notation (for readability of
the code).
Using a 1D vector should accelerate the access to the elements because
you only load one memory reference to the array and make a few easy
computations, instead of loading 2 diffrent addresses.
By the way, the library should be written to accept any type of elements.

If such a library does not exist, can one of you tell me how I should
construct my class declarations so that this could work ?
I tried to define operator [][], but this doesn't work.

Thanks for every help (via E-mail or News).

       Markus

---
-------------------------------------------------------------------------
Markus Schutz
E-Mail: schutz@ltssg3.epfl.ch
-------------------------------------------------------------------------

In article <2q2le7$3m6@info.epfl.ch>, <schutz@resumix.portal.com> writes:

> I am looking for a C++ library that implements handling of 2D (or nD)
> matrixes which are static at runtime (but not at compile time).
> My problem is that with standart C or C++ array definitions (i.e:
> <type> **A), the access to one element of A via A[x][y] is too slow.

The reason that this is slow is due to the intermediate object that you have to
create.  Try something like this:

inline TYPE operator() (int a, int b) {return p[a][b];}

> As at runtime the size of the matrix is known, all the extension features
> offered by the <type> **A declaration are not needed.
> I thought about library reserving memory for the matrix as if it was a
> unidimensional vector of size #lines * #columns, but letting me access
> to the elements with the standard 2D array notation (for readability of
> the code).

You could provide the same interface as above if you thought it would be
faster, but why do you think it would be.  For 2D you have *(p+a)+b, while in
1D you have p+(a*rows+b).

> Using a 1D vector should accelerate the access to the elements because
> you only load one memory reference to the array and make a few easy
> computations, instead of loading 2 diffrent addresses.
> By the way, the library should be written to accept any type of elements.
>

Oh, I see.  But your address loads are just pointers and pointer arithmetic as
I have shown above.  I still don't see where it is faster unless you are
talking about keeping the pipeline moving or something.

> If such a library does not exist, can one of you tell me how I should
> construct my class declarations so that this could work ?
> I tried to define operator [][], but this doesn't work.

Wait a minute, I thought you said that you tried this and it was too slow.  Now
you are saying you couldn't make it work.  It will be slower than the
operator() because of the intermediate object and therefore should be avoided.
 You want it to go faster not slower, right?

>
> Thanks for every help (via E-mail or News).
>
>        Markus
>
> ---
> -------------------------------------------------------------------------
> Markus Schutz
> E-Mail: schutz@ltssg3.epfl.ch
> -------------------------------------------------------------------------
>
>
>

In article <1994May2.134148.27544@tin.monsanto.com>,
 <raspar@ccmail.monsanto.com> wrote:
>
>In article <2q2le7$3m6@info.epfl.ch>, <schutz@resumix.portal.com> writes:
>
>> I am looking for a C++ library that implements handling of 2D (or nD)
>> matrixes which are static at runtime (but not at compile time).
>> My problem is that with standart C or C++ array definitions (i.e:
>> <type> **A), the access to one element of A via A[x][y] is too slow.
>
>The reason that this is slow is due to the intermediate object that you have to
>create.  Try something like this:
>
>inline TYPE operator() (int a, int b) {return p[a][b];}
>

Ughhh.

How about something simple like this

#include <stddef.h>
#include <iostream.h>

template <class T> class Matrix {
public:
  ~Matrix() {
    for (int i = 0; i < xRange; i++) {
      delete [] data[i];
    }
    delete [] data;
  }
  Matrix(size_t x = 1,size_t y = 1) :
    xRange(x), yRange(y), data(new double* [x]) {
    for (int i = 0; i < x ; i++) {
      data[i] = new double[y];
    }
  }
  T * operator [](int value) { return data[value]; }
  const T * operator [](int value) const { return data[value]; }
private:
  size_t xRange;
  size_t yRange;
  T ** data;
};

main()
{
  Matrix<double> XX(11,11);

  XX[10][5] = 10.1;

  cout << "XX[10][5] = " << (XX[10][5]) << endl;
}

Clearly this example won't let you do a 3 dimensional, so the simple way to
fix that is to have matrices implemented as vectors of vectors.

>Wait a minute, I thought you said that you tried this and it was too slow.  Now
>you are saying you couldn't make it work.  It will be slower than the
>operator() because of the intermediate object and therefore should be avoided.

This really depends on the compiler, and the architecture.  I think unless
you are thrashing the cache, an array of arrays (the normal 2D array in C)
will probably be quite faster, since this is something which compilers have
had to deal with for a long time.

Regards,
 Rohan

--
----------------------------------------------------------------------------
rjl@iassf.easams.com.au | "In the beginning there was cave man, and he had a
Rohan Lenard            |  club, which he used to seduce woman.  Hence the
+61-2-367-4555          |  saying `not tonight dear I've got a headache'"

In article <2q4237$r8a@f111.iassf.easams.com.au> rjl@f111.iassf.easams.com.au (Rohan LENARD) writes:

In article <1994May2.134148.27544@tin.monsanto.com>,
<raspar@ccmail.monsanto.com> wrote:
>
>In article <2q2le7$3m6@info.epfl.ch>, <schutz@resumix.portal.com> writes:
>
>> I am looking for a C++ library that implements handling of 2D (or nD)
>> matrixes which are static at runtime (but not at compile time).
>> My problem is that with standart C or C++ array definitions (i.e:
>> <type> **A), the access to one element of A via A[x][y] is too slow.
>

If you want the matrix to be allocated statically (rather than using
new), though this does not address the problem of access speed, it is
generally a _lot_ faster to place objects statically than to call new.
It is is possible for C++ programs which use new a lot become bound by
the heap usage overhead, so this may help.

template <class T, int SIZE>
class Vector {
protected:
  T rep[SIZE];
public:
  Vector();
  Vector(const T&);
  Vector(const Vector<T,SIZE>&);

  inline T& operator[](int i) { return rep[i]; }

// etc

};

and square matrix:

template <class T, int SIZE>
class Square_Matrix {
protected:
  Vector<T,SIZE> rep[SIZE];
public:
  Square_Matrix();
  Square_Matrix(const T&);
  Square_Matrix(const Vector<T,SIZE>&);
  Square_Matrix(const Matrix<T,SIZE>&);

  inline T& operator[](int i) { return rep[i]; }

// etc

};

The operator[] for square matrices returns a reference to the vector
representing the row, and the operator[] for vectors can then be used
on this object.  This allows the same effect as an operator[][]:

{
  Square_Matrix<double, 10> A;

  A[3][4] = 3.0;
}

this is equivalent to

{
  Square_Matrix<double, 10> A;

  A.operator[](3).operator[](4) = 3.0;
}

Square matrices should be ok using this approach, for instance a
multiply operation:

template <class T, int SIZE>
Square_Matrix<T,SIZE> operator*(const Square_Matrix<T,SIZE>& A,
    const Square_Matrix<T,SIZE>& B) {
  // code to multiply A*B
}

Non-square matrices are more difficult:

template <class T, int ROW, int COL>
class Matrix {
protected:
  Vector<T,COL> rep[ROW];
public:
  Matrix();
  Matrix(const T&);
  Matrix(const Vector<T,COL>&);
  Matrix(const Matrix<T,ROW,COL>&);
};

non-square matrix operations now look like:

template <class T, int ROW, int COL, int COL2>
Matrix<T,ROW,COL2> operator*(const Matrix<T,ROW,COL>& A,
        const Matrix<T,COL,COL2>& B) {
  // code to multiply A[ROW,COL] * B[COL,COL2] = RES[ROW,COL2]
}

Non-square matrices can handle square matrices too, and just as
efficiently as for square matrices.

Example use:

{
  Matrix<double,10,3> A;
  Matrix<double,3,5> B;
  Matrix<double,10,5> C;
  C = A * B;  // instantiates and calls [1]
}

[1]

template<double,10,3,5>
Matrix<double,10,5> operator*(const Matrix<double,10,3>&,
         const Matrix<double,3,5>&);

Of course you could get a bit of code explosion with the template
instantiations for these operations, as you get an instantiation for
each combination of sizes you use for each operation, but if
efficiency is important it may be worth it.

This definition of matrices allows the compiler to place the object
statically so if you rely on local matrices within functions a lot the
improvement in efficiency can be good as the local matrices go on the
stack.  The cost of allocating on the stack can be orders of magnitude
faster than calling new.

Adam
--
______________________________________________________________________________
A.Back@exeter.ac.uk (Internet)      (JANET) A.Back@uk.ac.exeter

In article <2q4237$r8a@f111.iassf.easams.com.au> rjl@f111.iassf.easams.com.au (Rohan LENARD) writes:

<stuff deleted>

> How about something simple like this
>
> #include <stddef.h>
> #include <iostream.h>
>
> template <class T> class Matrix {
> public:
>   ~Matrix() {
>     for (int i = 0; i < xRange; i++) {
>       delete [] data[i];
>     }
>     delete [] data;
>   }
>   Matrix(size_t x = 1,size_t y = 1) :
>     xRange(x), yRange(y), data(new double* [x]) {
>     for (int i = 0; i < x ; i++) {
>       data[i] = new double[y];
>     }
>   }
>   T * operator [](int value) { return data[value]; }
>   const T * operator [](int value) const { return data[value]; }

Are the above two lines legal?  My compiler (xlC) calls the above two
lines ambiguous when called on a non-const object.  I looked through
the ARM, and it wasn't obvious to me that the above are required to be
non-ambiguous.  In fact, the ARM explicitly states that a const member
function may be called on a const or non-const object.

> private:
>   size_t xRange;
>   size_t yRange;
>   T ** data;
> };

<stuff deleted>

Thanks

     -Scott Berg

--
J. Scott Berg       Real mail: Varian Physics; Stanford  CA  94305-4060
email: ALHY@slac.stanford.edu
phone: (415) 926-4732 (w)  (415) 326-2631 (h)

alhy@courant.SLAC.Stanford.EDU (J. Scott Berg) writes:

>In article <2q4237$r8a@f111.iassf.easams.com.au> rjl@f111.iassf.easams.com.au (Rohan LENARD) writes:

><stuff deleted>

>> How about something simple like this
>>
>> #include <stddef.h>
>> #include <iostream.h>
>>
>> template <class T> class Matrix {
>> public:
>>   ~Matrix() {
>>     for (int i = 0; i < xRange; i++) {
>>       delete [] data[i];
>>     }
>>     delete [] data;
>>   }
>>   Matrix(size_t x = 1,size_t y = 1) :
>>     xRange(x), yRange(y), data(new double* [x]) {
>>     for (int i = 0; i < x ; i++) {
>>       data[i] = new double[y];
>>     }
>>   }
>>   T * operator [](int value) { return data[value]; }
>>   const T * operator [](int value) const { return data[value]; }

>Are the above two lines legal?  My compiler (xlC) calls the above two
>lines ambiguous when called on a non-const object.  I looked through
>the ARM, and it wasn't obvious to me that the above are required to be
>non-ambiguous.  In fact, the ARM explicitly states that a const member
>function may be called on a const or non-const object.

>> private:
>>   size_t xRange;
>>   size_t yRange;
>>   T ** data;
>> };

><stuff deleted>

The ARM requires that C++ compilers properly differentiate
between const and non-const member functions with the same
name and arguments.

Unforunately, many current C++ compilers do not conform.

In article <2q67og$msv@fsgi01.fnal.gov>, b91926@fsgi01.fnal.gov (David Sachs) writes:
|> alhy@courant.SLAC.Stanford.EDU (J. Scott Berg) writes:
|>
|> >In article <2q4237$r8a@f111.iassf.easams.com.au> rjl@f111.iassf.easams.com.au (Rohan LENARD) writes:
|>
|> ><stuff deleted>
|>
|> >> How about something simple like this
|> >>
|> >> #include <stddef.h>
|> >> #include <iostream.h>
|> >>
|> >> template <class T> class Matrix {
|> >> public:
|> >>   ~Matrix() {
|> >>     for (int i = 0; i < xRange; i++) {
|> >>       delete [] data[i];
|> >>     }
|> >>     delete [] data;
|> >>   }
|> >>   Matrix(size_t x = 1,size_t y = 1) :
|> >>     xRange(x), yRange(y), data(new double* [x]) {
|> >>     for (int i = 0; i < x ; i++) {
|> >>       data[i] = new double[y];
|> >>     }
|> >>   }
|> >>   T * operator [](int value) { return data[value]; }
|> >>   const T * operator [](int value) const { return data[value]; }
|>
|> >Are the above two lines legal?  My compiler (xlC) calls the above two
|> >lines ambiguous when called on a non-const object.  I looked through
|> >the ARM, and it wasn't obvious to me that the above are required to be
|> >non-ambiguous.  In fact, the ARM explicitly states that a const member
|> >function may be called on a const or non-const object.
|>
|> >> private:
|> >>   size_t xRange;
|> >>   size_t yRange;
|> >>   T ** data;
|> >> };
|>
|> ><stuff deleted>
|>
|> The ARM requires that C++ compilers properly differentiate
|> between const and non-const member functions with the same
|> name and arguments.
|>
|> Unforunately, many current C++ compilers do not conform.

xlC was the name of IBM's C++ compiler on AIX in a former life.
The version complained about here must be very old. We have
used this approach in code for a long time.  Here is a version
of the "Matrix" code above, fixed for syntax errors, that compiles
with the current IBM compiler, now called C Set ++ (ugh):

 #include <stddef.h>
 #include <iostream.h>

 template <class T> class Matrix {
 public:
   ~Matrix() {
     for (int i = 0; i < xRange; i++) {
       delete [] data[i];
     }
     delete [] data;
   }
   Matrix(size_t x = 1,size_t y = 1) :
     xRange(x), yRange(y), data(new T* [x]) {
     for (int i = 0; i < x ; i++) {
       data[i] = new T[y];
     }
   }
   T * operator [](int value) { return data[value]; }
   const T * operator [](int value) const { return data[value]; }
 private:
   size_t xRange;
   size_t yRange;
   T ** data;
 };

int main() {
  Matrix<int> m(4,4);
  *(m[2]) = 1;
  const Matrix<int>& mr = m;
  cout << *(m[2]) << endl;
  return 0;
}


--
John.

John J. Barton        jjb@watson.ibm.com            (914)784-6645
H1-C13 IBM Watson Research Center P.O. Box 704 Hawthorne NY 10598