Looking at Type inferencing compiler

[David Bateman]

Jens,

I've had a quick look at your type inferencing compiler, more with the
idea of seeing how much work would be involved to get this to a point
where it can be included in octave, and be useful. I have several
stumbling blocks, but the first one is that you use a class oc_matrix
locally in your compiler and define the operator * in it, and then use
it. This class doesn't have reference counting, etc and is a very simple
matrix class.

Basically, I don't think its a good idea to recreate such additional
classes when octave itself has such classes available. The equivalent of
your class oc_matrix<OC_INT,10,10> would be
int32NDArray(dim_vector(10,10)). So as a test I've tried writing a cover
oct-file that uses code coming from your compiler for one example, using
int32NDArray for another and the standard Matrix array for the finally
version... One complication is that there is currently no * operator for
the octave integer matrices and so I added this.

I however notice large slowup in the code relative to your version. For
example

# Jens original cover code with 10x10 matrix
tic; cover(10000,1); toc
Elapsed time is 0.060677 seconds.
# cover using a int32NDArray and * operator based on xGEMM
tic; cover(10000,2); toc
Elapsed time is 0.398460 seconds.
# cover using the Matrix class and netlib blas
tic; cover(10000,3); toc
Elapsed time is 1.368492 seconds.
# cover using the Matrix class and ATLAS. SLOWER!!!!! matrix small at 10x10
tic; cover(10000,3); toc
Elapsed time is 2.240105 seconds.
# m-file version of cover
tic; cover(10000); toc
Elapsed time is 0.277359 seconds.

So the only example above that is faster than the m-file is yours!!! Do
you or anyone else know where I'm paying the speed penalty for the
octave in-built classes? In any case I suspect that the improvement you
showed in the example in your thesis was largely due to the different
matrix class that the type inferencing, though further benchmarks will
be needed to prove that.

I attach my cover.cc with the example code for the three implementation
together with my modified cover.m.

Regards
David

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

// C++ code followes

#include <octave/oct.h>
#include <iostream>

using namespace std;

template< class T, int I, int J> class oc_matrix
{
 public:
    T rep[I][J];
    int rows, cols;
    
    inline oc_matrix( T ( &_m )[I][J] )
        : rows(I), cols(J){ for( int i=0; i<I; ++i ) for( int j=0; j<J; ++j ) 
rep[i][j] = _m[i][j]; }

//   inline explicit oc_matrix( T ( &_m )[I][J] )
//      : rows(I), cols(J){ for( int i=0; i<I; ++i ) for( int j=0; j<J; ++j ) 
rep[i][j] = _m[i][j]; }

    inline explicit oc_matrix( T  &m )
        : rows(1), cols(1){ for( int i=0; i<I; ++i ) for( int j=0; j<J; ++j ) 
rep[i][j] = T(); rep[0][0] = m; }

    inline explicit oc_matrix( int i, int j )
        : rows(i), cols(j){}
    
    inline operator octave_value(){ 
        Matrix m(rows,cols); for (int i = 0; i<rows; ++i) for ( int j = 0 ; 
j<cols; ++j ) m(i,j) = rep[i][j];
        return m;
    };

};


template<class T, int I, int J> inline 
oc_matrix<T,I,J> operator*( oc_matrix<T,I,J>& a ,  oc_matrix<T,I,J>& b ){
    oc_matrix<T,I,J> r (I,J);
    if (a.rows == 1 and a.cols == 1) {
        for (int i=0; i<I; ++i)
            for (int j=0; j<I; ++j)
                r.rep[i][j] = a.rep[0][0] * b.rep[i][j];
        return r;
    }
    if (b.rows == 1 and b.rows == 1 ){
        for (int i=0; i<I; ++i)
            for (int j=0; j<I; ++j)
                r.rep[i][j] = a.rep[i][j] * b.rep[0][0];
        return r;
    }

    for (int i=0; i<I; ++i)
        for (int j=0; j<J; ++j){
            r.rep[i][j] = T();
            for (int z=0; z< (I); ++z)
                r.rep[i][j] += a.rep[i][z] * b.rep[z][j];
        }
    return r;

};    

template<class T, int I, int J> 
std::ostream& operator<< (std::ostream& s, const oc_matrix<T,I,J> m){
    for (int i=0; i<m.rows; ++i)
        for (int j=0; j<m.cols; ++j)
            s << m.rep[i][j] << " ";
    return s;
}

inline
octave_value operator<= (const octave_value& a, const double &i){
    if (a.is_scalar_type()){          
        double a_ = a.double_value();  
        return  octave_value(a_ <= i);  
    }                   
    return a <= octave_value(i);
}

#define OC_INT octave_idx_type

extern void
cover3 ( octave_value &n0 , octave_value& c0 )
{
//THE FUNCTION BODY
double c1 = 1;
double tmp0[10][10] = { 
{0,0,0,0,1,1,0,0,1,1},{0,1,1,1,1,0,1,0,0,0},{0,1,0,1,1,1,0,0,0,0},{0,1,1,1,1,0,1,0,0,1},{1,1,1,1,1,0,0,1,1,1},{1,0,1,0,0,0,1,0,0,0},{0,1,0,1,0,1,0,1,1,1},{0,0,0,0,1,0,1,1,1,0},{1,0,0,0,1,0,1,1,0,0},{1,0,0,1,1,0,1,0,0,1}
 };
Matrix a0 (10, 10);
double *a0_vec = a0.fortran_vec();
memcpy (a0_vec, tmp0, 100 *sizeof(double));
OC_INT i0 = 0;
bool oc_loop_control0 = true;
Matrix c2 (10, 10, 0);
double *c2_vec = c2.fortran_vec();
for (octave_idx_type i = 0; i < 100; i += 11)
  c2_vec[i] = 1;
do { 
c2 = c2 * a0;
i0 = i0 + 1;
oc_loop_control0 =  octave_value( n0 <= i0 ).is_true() ;
} while ( ! oc_loop_control0 );

c0 = c2;

}


// Define a matrix multiply operator for ints

template <class T>
intNDArray<T>
operator * (const intNDArray<T>& m, const intNDArray<T>& a) 
{
  intNDArray<T> retval;

  dim_vector mdv = m.dims();
  dim_vector adv = a.dims();

  if (mdv.length() != 2 || adv.length() != 2)
    gripe_nonconformant ("operator *", mdv, adv);
  else
    {
      octave_idx_type nr = m.rows ();
      octave_idx_type nc = m.cols ();

      octave_idx_type a_nr = a.rows ();
      octave_idx_type a_nc = a.cols ();

      if (nc != a_nr)
        gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
      else
        {
          if (nr == 0 || nc == 0 || a_nc == 0)
            retval.resize (dim_vector(nr, a_nc), 0.0);
          else
            {
              retval.resize(dim_vector(nr, a_nc));
              T *vec = retval.fortran_vec();
              const T *a_vec = a.fortran_vec();
              const T *m_vec = m.fortran_vec();

              for (octave_idx_type j = 0; j < a_nc; ++j)
                {
                  for (octave_idx_type i = 0; i < nr; ++i)
                    vec[i+j*nr] = T();
                  for (octave_idx_type l = 0; l < nc; ++l)
                    {
                      if (a_vec[l+j*nc] != T())
                        {
                          T temp = a_vec[l+j*nc];
                          for (octave_idx_type i = 0; i < nr; ++i)
                            vec[i+j*nr] += temp * m_vec[i+l*nc];
                        }
                    }
                }
            }
        }
    }
  return retval;
}

template
intNDArray<octave_int32> operator * (const intNDArray<octave_int32>& m, 
                                      const intNDArray<octave_int32>& a);

extern void
cover2 ( octave_value &n0 , octave_value& c0 )
{
//THE FUNCTION BODY
octave_int32 c1 = 1;
octave_int32 tmp0[10][10] = { 
{0,0,0,0,1,1,0,0,1,1},{0,1,1,1,1,0,1,0,0,0},{0,1,0,1,1,1,0,0,0,0},{0,1,1,1,1,0,1,0,0,1},{1,1,1,1,1,0,0,1,1,1},{1,0,1,0,0,0,1,0,0,0},{0,1,0,1,0,1,0,1,1,1},{0,0,0,0,1,0,1,1,1,0},{1,0,0,0,1,0,1,1,0,0},{1,0,0,1,1,0,1,0,0,1}
 };
int32NDArray a0 (dim_vector(10, 10));
octave_int32 *a0_vec = a0.fortran_vec();
memcpy (a0_vec, tmp0, 100 *sizeof(octave_int32));
OC_INT i0 = 0;
bool oc_loop_control0 = true;
int32NDArray c2 (dim_vector(10,10), 0);
octave_int32 *c2_vec = c2.fortran_vec();
for (octave_idx_type i = 0; i < 100; i += 11)
  c2_vec[i] = 1;
do { 
c2 = c2 * a0;
i0 = i0 + 1;
oc_loop_control0 =  octave_value( n0 <= i0 ).is_true() ;
} while ( ! oc_loop_control0 );

c0 = c2;

}

extern void
cover1 ( octave_value &n0 , octave_value& c0 )
{
//THE FUNCTION BODY
OC_INT c1 = 1;
OC_INT tmp0[10][10] = { 
{0,0,0,0,1,1,0,0,1,1},{0,1,1,1,1,0,1,0,0,0},{0,1,0,1,1,1,0,0,0,0},{0,1,1,1,1,0,1,0,0,1},{1,1,1,1,1,0,0,1,1,1},{1,0,1,0,0,0,1,0,0,0},{0,1,0,1,0,1,0,1,1,1},{0,0,0,0,1,0,1,1,1,0},{1,0,0,0,1,0,1,1,0,0},{1,0,0,1,1,0,1,0,0,1}
 };
oc_matrix<OC_INT,10,10> a0 = tmp0;
OC_INT i0 = 0;
bool oc_loop_control0 = true;
oc_matrix<OC_INT,10,10> c2( c1 );
do { 
c2 = c2 * a0;
i0 = i0 + 1;
oc_loop_control0 =  octave_value( n0 <= i0 ).is_true() ;
} while ( ! oc_loop_control0 );

c0 = c2;

}

DEFUN_DLD(cover, args, , "cover") {
     if(args.length() > 0 && args(0).is_real_scalar()) {
          octave_value a;
          octave_value n=args(0).double_value();
          if (args.length() > 1)
            {
              double k = args(1).double_value();
              if (k > 2.)
                cover3 (n, a);
              else if (k > 1.)
                cover2 (n, a);
              else
                cover1 (n, a);
            }
          else
            cover1 (n, a);
          return octave_value(a);
     }
     return octave_value(1);
}

function c = cover(n)
  c = 1;
  a = [1 1 0 0; 0 1 1 0; 0 0 1 1; 0 1 0 1];
  i = 0;
  do 
     c = c * a;
     i = i + 1;
  until (  n <= i );
endfunction