Logo Search packages:      
Sourcecode: blender version File versions

matrix.c

/*
 * $Id: matrix.c,v 1.13 2004/10/14 23:34:19 jesterking Exp $
 *
 * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version. The Blender
 * Foundation also sells licenses for use in proprietary software under
 * the Blender License.  See http://www.blender.org/BL/ for information
 * about this.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software Foundation,
 * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
 * All rights reserved.
 *
 * Contributor(s): Michel Selten & Joseph Gilbert
 *
 * ***** END GPL/BL DUAL LICENSE BLOCK *****
 */

#include "matrix.h"

//doc strings
char Matrix_Zero_doc[] = "() - set all values in the matrix to 0";
char Matrix_Identity_doc[] =
      "() - set the square matrix to it's identity matrix";
char Matrix_Transpose_doc[] = "() - set the matrix to it's transpose";
char Matrix_Determinant_doc[] = "() - return the determinant of the matrix";
char Matrix_Invert_doc[] =
      "() - set the matrix to it's inverse if an inverse is possible";
char Matrix_TranslationPart_doc[] =
      "() - return a vector encompassing the translation of the matrix";
char Matrix_RotationPart_doc[] =
      "() - return a vector encompassing the rotation of the matrix";
char Matrix_Resize4x4_doc[] = "() - resize the matrix to a 4x4 square matrix";
char Matrix_toEuler_doc[] = "() - convert matrix to a euler angle rotation";
char Matrix_toQuat_doc[] = "() - convert matrix to a quaternion rotation";

//methods table
struct PyMethodDef Matrix_methods[] = {
      {"zero", ( PyCFunction ) Matrix_Zero, METH_NOARGS,
       Matrix_Zero_doc},
      {"identity", ( PyCFunction ) Matrix_Identity, METH_NOARGS,
       Matrix_Identity_doc},
      {"transpose", ( PyCFunction ) Matrix_Transpose, METH_NOARGS,
       Matrix_Transpose_doc},
      {"determinant", ( PyCFunction ) Matrix_Determinant, METH_NOARGS,
       Matrix_Determinant_doc},
      {"invert", ( PyCFunction ) Matrix_Invert, METH_NOARGS,
       Matrix_Invert_doc},
      {"translationPart", ( PyCFunction ) Matrix_TranslationPart,
       METH_NOARGS,
       Matrix_TranslationPart_doc},
      {"rotationPart", ( PyCFunction ) Matrix_RotationPart, METH_NOARGS,
       Matrix_RotationPart_doc},
      {"resize4x4", ( PyCFunction ) Matrix_Resize4x4, METH_NOARGS,
       Matrix_Resize4x4_doc},
      {"toEuler", ( PyCFunction ) Matrix_toEuler, METH_NOARGS,
       Matrix_toEuler_doc},
      {"toQuat", ( PyCFunction ) Matrix_toQuat, METH_NOARGS,
       Matrix_toQuat_doc},
      {NULL, NULL, 0, NULL}
};

/*****************************/
//    Matrix Python Object   
/*****************************/

PyObject *Matrix_toQuat( MatrixObject * self )
{
      float *quat, *mat;

      if( self->colSize < 3 ) {
            return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "inappropriate matrix size\n" );
      } else if( self->colSize > 2 ) {    //3 or 4 col
            if( self->rowSize < 3 ) //3 or 4 row
                  return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                          "inappropriate matrix size\n" );

            mat = PyMem_Malloc( 3 * 3 * sizeof( float ) );
            if( mat == NULL ) {
                  return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                          "problem allocating matrix\n\n" ) );
            }
            mat[0] = self->matrix[0][0];
            mat[1] = self->matrix[0][1];
            mat[2] = self->matrix[0][2];
            mat[3] = self->matrix[1][0];
            mat[4] = self->matrix[1][1];
            mat[5] = self->matrix[1][2];
            mat[6] = self->matrix[2][0];
            mat[7] = self->matrix[2][1];
            mat[8] = self->matrix[2][2];
      }
      quat = PyMem_Malloc( 4 * sizeof( float ) );
      if( quat == NULL ) {
            PyMem_Free( mat );
            return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                    "problem allocating quat\n\n" ) );
      }
      Mat3ToQuat( ( float ( * )[3] ) mat, quat );

      return ( PyObject * ) newQuaternionObject( quat );
}


PyObject *Matrix_toEuler( MatrixObject * self )
{
      float *eul, *mat;
      int x;

      if( self->colSize < 3 ) {
            return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "inappropriate matrix size\n" );
      } else if( self->colSize > 2 ) {    //3 or 4 col
            if( self->rowSize < 3 ) //3 or 4 row
                  return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                          "inappropriate matrix size\n" );

            mat = PyMem_Malloc( 3 * 3 * sizeof( float ) );
            if( mat == NULL ) {
                  return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                          "problem allocating mat\n\n" ) );
            }
            mat[0] = self->matrix[0][0];
            mat[1] = self->matrix[0][1];
            mat[2] = self->matrix[0][2];
            mat[3] = self->matrix[1][0];
            mat[4] = self->matrix[1][1];
            mat[5] = self->matrix[1][2];
            mat[6] = self->matrix[2][0];
            mat[7] = self->matrix[2][1];
            mat[8] = self->matrix[2][2];
      }
      eul = PyMem_Malloc( 3 * sizeof( float ) );
      if( eul == NULL ) {
            PyMem_Free( mat );
            return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                    "problem allocating eul\n\n" ) );
      }
      Mat3ToEul( ( float ( * )[3] ) mat, eul );

      for( x = 0; x < 3; x++ ) {
            eul[x] *= ( float ) ( 180 / Py_PI );
      }

      return ( PyObject * ) newEulerObject( eul );
}

PyObject *Matrix_Resize4x4( MatrixObject * self )
{
      float *mat;
      int x, row, col;

      if( self->colSize == 4 && self->rowSize == 4 )
            return EXPP_incr_ret( Py_None );

      mat = PyMem_Malloc( 4 * 4 * sizeof( float ) );
      if( mat == NULL ) {
            return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                    "problem allocating mat\n\n" ) );
      }
      for( x = 0; x < 16; x++ ) {
            mat[x] = 0.0f;
      }

      if( self->colSize == 2 ) {    //2x2, 2x3, 2x4
            mat[0] = self->matrix[0][0];
            mat[1] = self->matrix[0][1];
            mat[4] = self->matrix[1][0];
            mat[5] = self->matrix[1][1];
            if( self->rowSize > 2 ) {
                  mat[8] = self->matrix[2][0];
                  mat[9] = self->matrix[2][1];
            }
            if( self->rowSize > 3 ) {
                  mat[12] = self->matrix[3][0];
                  mat[13] = self->matrix[3][1];
            }
            mat[10] = 1.0f;
            mat[15] = 1.0f;
      } else if( self->colSize == 3 ) {   //3x2, 3x3, 3x4
            mat[0] = self->matrix[0][0];
            mat[1] = self->matrix[0][1];
            mat[2] = self->matrix[0][2];
            mat[4] = self->matrix[1][0];
            mat[5] = self->matrix[1][1];
            mat[6] = self->matrix[1][2];
            if( self->rowSize > 2 ) {
                  mat[8] = self->matrix[2][0];
                  mat[9] = self->matrix[2][1];
                  mat[10] = self->matrix[2][2];
            }
            if( self->rowSize > 3 ) {
                  mat[12] = self->matrix[3][0];
                  mat[13] = self->matrix[3][1];
                  mat[14] = self->matrix[3][2];
            }
            if( self->rowSize == 2 )
                  mat[10] = 1.0f;
            mat[15] = 1.0f;
      } else if( self->colSize == 4 ) {   //2x4, 3x4
            mat[0] = self->matrix[0][0];
            mat[1] = self->matrix[0][1];
            mat[2] = self->matrix[0][2];
            mat[3] = self->matrix[0][3];
            mat[4] = self->matrix[1][0];
            mat[5] = self->matrix[1][1];
            mat[6] = self->matrix[1][2];
            mat[7] = self->matrix[1][3];
            if( self->rowSize > 2 ) {
                  mat[8] = self->matrix[2][0];
                  mat[9] = self->matrix[2][1];
                  mat[10] = self->matrix[2][2];
                  mat[11] = self->matrix[2][3];
            }
            if( self->rowSize == 2 )
                  mat[10] = 1.0f;
            mat[15] = 1.0f;
      }

      PyMem_Free( self->matrix );
      PyMem_Free( self->contigPtr );
      self->contigPtr = PyMem_Malloc( 4 * 4 * sizeof( float ) );
      if( self->contigPtr == NULL ) {
            return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                    "problem allocating array space\n\n" ) );
      }
      self->matrix = PyMem_Malloc( 4 * sizeof( float * ) );
      if( self->matrix == NULL ) {
            PyMem_Free( self->contigPtr );
            return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                    "problem allocating pointer space\n\n" ) );
      }
      for( x = 0; x < 4; x++ ) {
            self->matrix[x] = self->contigPtr + ( x * 4 );
      }

      for( row = 0; row < 4; row++ ) {
            for( col = 0; col < 4; col++ ) {
                  self->matrix[row][col] = mat[( row * 4 ) + col];
            }
      }
      PyMem_Free( mat );

      self->colSize = 4;
      self->rowSize = 4;

      return EXPP_incr_ret( Py_None );
}

PyObject *Matrix_TranslationPart( MatrixObject * self )
{
      float *vec = NULL;
      PyObject *retval;

      if( self->colSize < 3 ) {
            return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "inappropriate matrix size\n" );
      } else if( self->colSize > 2 ) {    //3 or 4 columns
            if( self->rowSize < 4 ) //all 4 rows
                  return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                          "inappropriate matrix size\n" );

            vec = PyMem_Malloc( 3 * sizeof( float ) );
            if( vec == NULL ) {
                  return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                          "problem allocating vec\n\n" ) );
            }
            vec[0] = self->matrix[3][0];
            vec[1] = self->matrix[3][1];
            vec[2] = self->matrix[3][2];
      }

      retval =  ( PyObject * ) newVectorObject( vec, 3 );
      PyMem_Free( vec );
      return retval;
}

PyObject *Matrix_RotationPart( MatrixObject * self )
{
      float *mat;

      if( self->colSize < 3 ) {
            return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "inappropriate matrix size\n" );
      } else if( self->colSize > 2 ) {    //3 or 4 col
            if( self->rowSize < 3 ) //3 or 4 row
                  return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                          "inappropriate matrix size\n" );

            mat = PyMem_Malloc( 3 * 3 * sizeof( float ) );
            if( mat == NULL ) {
                  return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                          "problem allocating mat\n\n" ) );
            }
            mat[0] = self->matrix[0][0];
            mat[1] = self->matrix[0][1];
            mat[2] = self->matrix[0][2];
            mat[3] = self->matrix[1][0];
            mat[4] = self->matrix[1][1];
            mat[5] = self->matrix[1][2];
            mat[6] = self->matrix[2][0];
            mat[7] = self->matrix[2][1];
            mat[8] = self->matrix[2][2];
      }

      return ( PyObject * ) newMatrixObject( mat, 3, 3 );
}

PyObject *Matrix_Invert( MatrixObject * self )
{
      float det;
      int x, y, z;
      float *mat = NULL;
      float t;

      if( self->rowSize != self->colSize )
            return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "only square matrices are supported\n" );

      //calculate the determinant
      if( self->rowSize == 2 ) {
            det = Det2x2( self->matrix[0][0], self->matrix[0][1],
                        self->matrix[1][0], self->matrix[1][1] );
      } else if( self->rowSize == 3 ) {
            det = Det3x3( self->matrix[0][0], self->matrix[0][1],
                        self->matrix[0][2], self->matrix[1][0],
                        self->matrix[1][1], self->matrix[1][2],
                        self->matrix[2][0], self->matrix[2][1],
                        self->matrix[2][2] );
      } else if( self->rowSize == 4 ) {
            det = Det4x4( *self->matrix );
      } else {
            return EXPP_ReturnPyObjError( PyExc_StandardError,
                                    "error calculating determinant for inverse()\n" );
      }

      if( det != 0 ) {

            //calculate the classical adjoint
            if( self->rowSize == 2 ) {
                  mat = PyMem_Malloc( self->rowSize * self->colSize *
                                  sizeof( float ) );
                  if( mat == NULL ) {
                        return ( EXPP_ReturnPyObjError
                               ( PyExc_MemoryError,
                                 "problem allocating mat\n\n" ) );
                  }
                  mat[0] = self->matrix[1][1];
                  mat[1] = -self->matrix[1][0];
                  mat[2] = -self->matrix[0][1];
                  mat[3] = self->matrix[0][0];
            } else if( self->rowSize == 3 ) {
                  mat = PyMem_Malloc( self->rowSize * self->colSize *
                                  sizeof( float ) );
                  if( mat == NULL ) {
                        return ( EXPP_ReturnPyObjError
                               ( PyExc_MemoryError,
                                 "problem allocating mat\n\n" ) );
                  }
                  Mat3Adj( ( float ( * )[3] ) mat, *self->matrix );
            } else if( self->rowSize == 4 ) {
                  mat = PyMem_Malloc( self->rowSize * self->colSize *
                                  sizeof( float ) );
                  if( mat == NULL ) {
                        return ( EXPP_ReturnPyObjError
                               ( PyExc_MemoryError,
                                 "problem allocating mat\n\n" ) );
                  }
                  Mat4Adj( ( float ( * )[4] ) mat, *self->matrix );
            }
            //divide by determinate
            for( x = 0; x < ( self->rowSize * self->colSize ); x++ ) {
                  mat[x] /= det;
            }

            //set values
            z = 0;
            for( x = 0; x < self->rowSize; x++ ) {
                  for( y = 0; y < self->colSize; y++ ) {
                        self->matrix[x][y] = mat[z];
                        z++;
                  }
            }

            //transpose
            if( self->rowSize == 2 ) {
                  t = self->matrix[1][0];
                  self->matrix[1][0] = self->matrix[0][1];
                  self->matrix[0][1] = t;
                  
/* 
   Note: is the code below correct?  
   transposing mat and not copying into self->matrix? 
   s. swaney 11-oct-2004
*/
            } else if( self->rowSize == 3 ) {
                  Mat3Transp( ( float ( * )[3] ) mat );
            } else if( self->rowSize == 4 ) {
                  Mat4Transp( ( float ( * )[4] ) mat );
            }
      } else {
            printf( "matrix does not have an inverse - none attempted\n" );
      }
      PyMem_Free( mat );
      return EXPP_incr_ret( Py_None );
}


PyObject *Matrix_Determinant( MatrixObject * self )
{
      float det;

      if( self->rowSize != self->colSize )
            return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "only square matrices are supported\n" );

      if( self->rowSize == 2 ) {
            det = Det2x2( self->matrix[0][0], self->matrix[0][1],
                        self->matrix[1][0], self->matrix[1][1] );
      } else if( self->rowSize == 3 ) {
            det = Det3x3( self->matrix[0][0], self->matrix[0][1],
                        self->matrix[0][2], self->matrix[1][0],
                        self->matrix[1][1], self->matrix[1][2],
                        self->matrix[2][0], self->matrix[2][1],
                        self->matrix[2][2] );
      } else if( self->rowSize == 4 ) {
            det = Det4x4( *self->matrix );
      } else {
            return EXPP_ReturnPyObjError( PyExc_StandardError,
                                    "error in determinant()\n" );
      }
      return PyFloat_FromDouble( det );
}

PyObject *Matrix_Transpose( MatrixObject * self )
{
      float t;

      if( self->rowSize != self->colSize )
            return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "only square matrices are supported\n" );

      if( self->rowSize == 2 ) {
            t = self->matrix[1][0];
            self->matrix[1][0] = self->matrix[0][1];
            self->matrix[0][1] = t;
      } else if( self->rowSize == 3 ) {
            Mat3Transp( *self->matrix );
      } else if( self->rowSize == 4 ) {
            Mat4Transp( *self->matrix );
      } else
            return ( EXPP_ReturnPyObjError( PyExc_TypeError,
                                    "unable to transpose matrix\n" ) );

      return EXPP_incr_ret( Py_None );
}

PyObject *Matrix_Zero( MatrixObject * self )
{
      int row, col;

      for( row = 0; row < self->rowSize; row++ ) {
            for( col = 0; col < self->colSize; col++ ) {
                  self->matrix[row][col] = 0.0f;
            }
      }
      return EXPP_incr_ret( Py_None );
}

PyObject *Matrix_Identity( MatrixObject * self )
{
      if( self->rowSize != self->colSize )
            return ( EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "only square matrices supported\n" ) );

      if( self->rowSize == 2 ) {
            self->matrix[0][0] = 1.0f;
            self->matrix[0][1] = 0.0f;
            self->matrix[1][0] = 0.0f;
            self->matrix[1][1] = 1.0f;
      } else if( self->rowSize == 3 ) {
            Mat3One( *self->matrix );
      } else if( self->rowSize == 4 ) {
            Mat4One( *self->matrix );
      } else
            return ( EXPP_ReturnPyObjError( PyExc_TypeError,
                                    "unable to create identity matrix\n" ) );

      return EXPP_incr_ret( Py_None );
}

static void Matrix_dealloc( MatrixObject * self )
{
      PyMem_Free( self->contigPtr );
      PyMem_Free( self->matrix );

      PyObject_DEL( self );
}

static PyObject *Matrix_getattr( MatrixObject * self, char *name )
{
      if( strcmp( name, "rowSize" ) == 0 ) {
            return PyInt_FromLong( ( long ) self->rowSize );
      } else if( strcmp( name, "colSize" ) == 0 ) {
            return PyInt_FromLong( ( long ) self->colSize );
      }

      return Py_FindMethod( Matrix_methods, ( PyObject * ) self, name );
}

static int Matrix_setattr( MatrixObject * self, char *name, PyObject * v )
{
      /* This is not supported. */
      return ( -1 );
}

static PyObject *Matrix_repr( MatrixObject * self )
{
      PyObject *repr, *str;
      int x, y;
      char ftoa[24];

      repr = PyString_FromString( "" );
      if( !repr )
            return ( EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "Attribute error in PyMatrix (repr)\n" ) );

      for( x = 0; x < self->rowSize; x++ ) {
            str = PyString_FromString( "[" );
            PyString_ConcatAndDel( &repr, str );

            for( y = 0; y < ( self->colSize - 1 ); y++ ) {
                  sprintf( ftoa, "%.4f, ", self->matrix[x][y] );
                  str = PyString_FromString( ftoa );
                  PyString_ConcatAndDel( &repr, str );
            }
            sprintf( ftoa, "%.4f]\n", self->matrix[x][y] );
            str = PyString_FromString( ftoa );
            PyString_ConcatAndDel( &repr, str );
      }
      return repr;
}

//no support for matrix[x][y] so have to return by sequence index
//will return a row from the matrix to support previous API
//compatability
static PyObject *Matrix_item( MatrixObject * self, int i )
{
      float *vec = NULL;
      PyObject *retval;
      int x;

      if( i < 0 || i >= self->rowSize )
            return EXPP_ReturnPyObjError( PyExc_IndexError,
                                    "matrix row index out of range\n" );

      vec = PyMem_Malloc( self->colSize * sizeof( float ) );
      if( vec == NULL ) {
            return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                    "problem allocating vec\n\n" ) );
      }
      for( x = 0; x < self->colSize; x++ ) {
            vec[x] = self->matrix[i][x];
      }

      retval =( PyObject * ) newVectorObject( vec, self->colSize );
      PyMem_Free( vec );
      return retval;
}

static PyObject *Matrix_slice( MatrixObject * self, int begin, int end )
{
      PyObject *list;
      int count, maxsize, x, y;

      maxsize = self->colSize * self->rowSize;
      if( begin < 0 )
            begin = 0;
      if( end > maxsize )
            end = maxsize;
      if( begin > end )
            begin = end;

      list = PyList_New( end - begin );

      for( count = begin; count < end; count++ ) {
            x = ( int ) floor( ( double ) ( count / self->colSize ) );
            y = count % self->colSize;
            PyList_SetItem( list, count - begin,
                        PyFloat_FromDouble( self->matrix[x][y] ) );
      }

      return list;
}

static int Matrix_ass_item( MatrixObject * self, int i, PyObject * ob )
{
      int maxsize, x, y;

      maxsize = self->colSize * self->rowSize;
      if( i < 0 || i >= maxsize )
            return EXPP_ReturnIntError( PyExc_IndexError,
                                  "array assignment index out of range\n" );
      if( !PyInt_Check( ob ) && !PyFloat_Check( ob ) )
            return EXPP_ReturnIntError( PyExc_IndexError,
                                  "matrix member must be a number\n" );

      x = ( int ) floor( ( double ) ( i / self->colSize ) );
      y = i % self->colSize;
      self->matrix[x][y] = ( float ) PyFloat_AsDouble( ob );

      return 0;
}

static int Matrix_ass_slice( MatrixObject * self, int begin, int end,
                       PyObject * seq )
{
      int count, maxsize, x, y, z;

      maxsize = self->colSize * self->rowSize;
      if( begin < 0 )
            begin = 0;
      if( end > maxsize )
            end = maxsize;
      if( begin > end )
            begin = end;

      if( !PySequence_Check( seq ) )
            return EXPP_ReturnIntError( PyExc_TypeError,
                                  "illegal argument type for built-in operation\n" );
      if( PySequence_Length( seq ) != ( end - begin ) )
            return EXPP_ReturnIntError( PyExc_TypeError,
                                  "size mismatch in slice assignment\n" );

      z = 0;
      for( count = begin; count < end; count++ ) {
            PyObject *ob = PySequence_GetItem( seq, z );
            z++;
            if( !PyInt_Check( ob ) && !PyFloat_Check( ob ) )
                  return EXPP_ReturnIntError( PyExc_IndexError,
                                        "list member must be a number\n" );

            x = ( int ) floor( ( double ) ( count / self->colSize ) );
            y = count % self->colSize;
            if( !PyArg_Parse( ob, "f", &self->matrix[x][y] ) ) {
                  Py_DECREF( ob );
                  return -1;
            }
      }
      return 0;
}

static int Matrix_len( MatrixObject * self )
{
      return ( self->colSize * self->rowSize );
}

PyObject *Matrix_add( PyObject * m1, PyObject * m2 )
{
      float *mat;
      int matSize, rowSize, colSize, x, y;

      if( ( !Matrix_CheckPyObject( m1 ) )
          || ( !Matrix_CheckPyObject( m2 ) ) )
            return EXPP_ReturnPyObjError( PyExc_TypeError,
                                    "unsupported type for this operation\n" );

      if( ( ( MatrixObject * ) m1 )->flag > 0
          || ( ( MatrixObject * ) m2 )->flag > 0 )
            return EXPP_ReturnPyObjError( PyExc_ArithmeticError,
                                    "cannot add scalar to a matrix\n" );

      if( ( ( MatrixObject * ) m1 )->rowSize !=
          ( ( MatrixObject * ) m2 )->rowSize
          || ( ( MatrixObject * ) m1 )->colSize !=
          ( ( MatrixObject * ) m2 )->colSize )
            return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "matrices must be the same same for this operation\n" );

      rowSize = ( ( ( MatrixObject * ) m1 )->rowSize );
      colSize = ( ( ( MatrixObject * ) m1 )->colSize );
      matSize = rowSize * colSize;

      mat = PyMem_Malloc( matSize * sizeof( float ) );
      if( mat == NULL ) {
            return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                    "problem allocating mat\n\n" ) );
      }
      for( x = 0; x < rowSize; x++ ) {
            for( y = 0; y < colSize; y++ ) {
                  mat[( ( x * rowSize ) + y )] =
                        ( ( MatrixObject * ) m1 )->matrix[x][y] +
                        ( ( MatrixObject * ) m2 )->matrix[x][y];
            }
      }

      return newMatrixObject( mat, rowSize, colSize );
}

PyObject *Matrix_sub( PyObject * m1, PyObject * m2 )
{
      float *mat;
      int matSize, rowSize, colSize, x, y;

      if( ( !Matrix_CheckPyObject( m1 ) )
          || ( !Matrix_CheckPyObject( m2 ) ) )
            return EXPP_ReturnPyObjError( PyExc_TypeError,
                                    "unsupported type for this operation\n" );

      if( ( ( MatrixObject * ) m1 )->flag > 0
          || ( ( MatrixObject * ) m2 )->flag > 0 )
            return EXPP_ReturnPyObjError( PyExc_ArithmeticError,
                                    "cannot subtract a scalar from a matrix\n" );

      if( ( ( MatrixObject * ) m1 )->rowSize !=
          ( ( MatrixObject * ) m2 )->rowSize
          || ( ( MatrixObject * ) m1 )->colSize !=
          ( ( MatrixObject * ) m2 )->colSize )
            return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "matrices must be the same same for this operation\n" );

      rowSize = ( ( ( MatrixObject * ) m1 )->rowSize );
      colSize = ( ( ( MatrixObject * ) m1 )->colSize );
      matSize = rowSize * colSize;

      mat = PyMem_Malloc( matSize * sizeof( float ) );
      if( mat == NULL ) {
            return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                    "problem allocating mat\n\n" ) );
      }
      for( x = 0; x < rowSize; x++ ) {
            for( y = 0; y < colSize; y++ ) {
                  mat[( ( x * rowSize ) + y )] =
                        ( ( MatrixObject * ) m1 )->matrix[x][y] -
                        ( ( MatrixObject * ) m2 )->matrix[x][y];
            }
      }

      return newMatrixObject( mat, rowSize, colSize );
}

PyObject *Matrix_mul( PyObject * m1, PyObject * m2 )
{
      PyObject *retval;
      int matSizeV, rowSizeV, colSizeV, rowSizeW, colSizeW, matSizeW, x, y,z;
      MatrixObject *matV;
      MatrixObject *matW;
      float *mat = NULL;
      float dot = 0;
      

      if( ( !Matrix_CheckPyObject( m1 ) )
          || ( !Matrix_CheckPyObject( m2 ) ) )
            return EXPP_ReturnPyObjError( PyExc_TypeError,
                                    "unsupported type for this operation\n" );

      //get some vars
      rowSizeV = ( ( ( MatrixObject * ) m1 )->rowSize );
      colSizeV = ( ( ( MatrixObject * ) m1 )->colSize );
      matSizeV = rowSizeV * colSizeV;
      rowSizeW = ( ( ( MatrixObject * ) m2 )->rowSize );
      colSizeW = ( ( ( MatrixObject * ) m2 )->colSize );
      matSizeW = rowSizeW * colSizeW;
      matV = ( ( MatrixObject * ) m1 );
      matW = ( ( MatrixObject * ) m2 );

      //coerced int or float for scalar multiplication
      if( matW->flag > 1 || matW->flag > 2 ) {

            if( rowSizeV != rowSizeW && colSizeV != colSizeW )
                  return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                          "Matrix dimension error during scalar multiplication\n" );

            mat = PyMem_Malloc( matSizeV * sizeof( float ) );
            if( mat == NULL ) {
                  return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                          "problem allocating mat\n\n" ) );
            }
            for( x = 0; x < rowSizeV; x++ ) {
                  for( y = 0; y < colSizeV; y++ ) {
                        mat[( ( x * rowSizeV ) + y )] =
                              matV->matrix[x][y] *
                              matW->matrix[x][y];
                  }
            }
            retval = ( PyObject* ) newMatrixObject( mat, rowSizeV, colSizeV );
            PyMem_Free( mat );
            return retval;
      } else if( matW->flag == 0 && matV->flag == 0 ) {     //true matrix multiplication
            if( colSizeV != rowSizeW ) {
                  return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                          "Matrix multiplication undefined...\n" );
            }

            mat = PyMem_Malloc( ( rowSizeV * colSizeW ) *
                            sizeof( float ) );
            if( mat == NULL ) {
                  return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                          "problem allocating mat\n\n" ) );
            }
            for( x = 0; x < rowSizeV; x++ ) {
                  for( y = 0; y < colSizeW; y++ ) {
                        for( z = 0; z < colSizeV; z++ ) {
                              dot += ( matV->matrix[x][z] *
                                     matW->matrix[z][y] );
                        }
                        mat[( ( x * rowSizeV ) + y )] = dot;
                        dot = 0;
                  }
            }
            retval = ( PyObject* )  newMatrixObject( mat, rowSizeV, colSizeW );
            PyMem_Free( mat );
            return retval;
      } else
            return EXPP_ReturnPyObjError( PyExc_AttributeError,
                                    "Error in matrix_mul...\n" );
}

//coercion of unknown types to type MatrixObject for numeric protocols
int Matrix_coerce( PyObject ** m1, PyObject ** m2 )
{
      long *tempI;
      double *tempF;
      float *mat;
      int x, matSize;

      matSize =
            ( ( ( MatrixObject * ) * m1 )->rowSize ) *
            ( ( ( MatrixObject * ) * m1 )->rowSize );
      if( Matrix_CheckPyObject( *m1 ) ) {
            if( Matrix_CheckPyObject( *m2 ) ) { //matrix & matrix
                  Py_INCREF( *m1 );
                  Py_INCREF( *m2 );
                  return 0;
            } else {
                  if( VectorObject_Check( *m2 ) ) {   //matrix & vector?
                        printf( "use MatMultVec() for column vector multiplication\n" );
                        Py_INCREF( *m1 );
                        return 0;
                  } else if( PyNumber_Check( *m2 ) ) {      //& scalar?
                        if( PyInt_Check( *m2 ) ) {    //it's a int
                              tempI = PyMem_Malloc( 1 *
                                                sizeof( long ) );
                              if( tempI == NULL ) {
                                    return ( EXPP_ReturnIntError
                                           ( PyExc_MemoryError,
                                             "problem allocating tempI\n\n" ) );
                              }
                              *tempI = PyInt_AsLong( *m2 );
                              mat = PyMem_Malloc( matSize *
                                              sizeof( float ) );
                              if( mat == NULL ) {
                                    PyMem_Free( tempI );
                                    return ( EXPP_ReturnIntError
                                           ( PyExc_MemoryError,
                                             "problem allocating mat\n\n" ) );
                              }
                              for( x = 0; x < matSize; x++ ) {
                                    mat[x] = ( float ) *tempI;
                              }
                              PyMem_Free( tempI );
                              *m2 = newMatrixObject( mat,
                                                 ( ( ( MatrixObject * ) * m1 )->rowSize ), ( ( ( MatrixObject * ) * m1 )->colSize ) );
                              ( ( MatrixObject * ) * m2 )->flag = 1;    //int coercion
                              PyMem_Free( mat );
                              Py_INCREF( *m1 );
                              return 0;
                        } else if( PyFloat_Check( *m2 ) ) { //it's a float
                              tempF = PyMem_Malloc( 1 *
                                                sizeof
                                                ( double ) );
                              if( tempF == NULL ) {
                                    return ( EXPP_ReturnIntError
                                           ( PyExc_MemoryError,
                                             "problem allocating tempF\n\n" ) );
                              }
                              *tempF = PyFloat_AsDouble( *m2 );
                              mat = PyMem_Malloc( matSize *
                                              sizeof( float ) );
                              if( mat == NULL ) {
                                    PyMem_Free( tempF );
                                    return ( EXPP_ReturnIntError
                                           ( PyExc_MemoryError,
                                             "problem allocating mat\n\n" ) );
                              }
                              for( x = 0; x < matSize; x++ ) {
                                    mat[x] = ( float ) *tempF;
                              }
                              PyMem_Free( tempF );
                              *m2 = newMatrixObject( mat,
                                                 ( ( ( MatrixObject * ) * m1 )->rowSize ), ( ( ( MatrixObject * ) * m1 )->colSize ) );
                              ( ( MatrixObject * ) * m2 )->flag = 2;    //float coercion
                              PyMem_Free( mat );
                              Py_INCREF( *m1 );
                              return 0;
                        }
                  }
                  //unknom2n type or numeric cast failure
                  printf( "attempting matrix operation m2ith unsupported type...\n" );
                  Py_INCREF( *m1 );
                  return 0;   //operation m2ill type check
            }
      } else {
            //1st not Matrix
            printf( "numeric protocol failure...\n" );
            return -1;  //this should not occur - fail
      }
      return -1;
}

//******************************************************************
//                                      Matrix definition
//******************************************************************
static PySequenceMethods Matrix_SeqMethods = {
      ( inquiry ) Matrix_len, /* sq_length */
      ( binaryfunc ) 0, /* sq_concat */
      ( intargfunc ) 0, /* sq_repeat */
      ( intargfunc ) Matrix_item,   /* sq_item */
      ( intintargfunc ) Matrix_slice,     /* sq_slice */
      ( intobjargproc ) Matrix_ass_item,  /* sq_ass_item */
      ( intintobjargproc ) Matrix_ass_slice,    /* sq_ass_slice */
};

static PyNumberMethods Matrix_NumMethods = {
      ( binaryfunc ) Matrix_add,    /* __add__ */
      ( binaryfunc ) Matrix_sub,    /* __sub__ */
      ( binaryfunc ) Matrix_mul,    /* __mul__ */
      ( binaryfunc ) 0, /* __div__ */
      ( binaryfunc ) 0, /* __mod__ */
      ( binaryfunc ) 0, /* __divmod__ */
      ( ternaryfunc ) 0,      /* __pow__ */
      ( unaryfunc ) 0,  /* __neg__ */
      ( unaryfunc ) 0,  /* __pos__ */
      ( unaryfunc ) 0,  /* __abs__ */
      ( inquiry ) 0,          /* __nonzero__ */
      ( unaryfunc ) 0,  /* __invert__ */
      ( binaryfunc ) 0, /* __lshift__ */
      ( binaryfunc ) 0, /* __rshift__ */
      ( binaryfunc ) 0, /* __and__ */
      ( binaryfunc ) 0, /* __xor__ */
      ( binaryfunc ) 0, /* __or__ */
      ( coercion ) Matrix_coerce,   /* __coerce__ */
      ( unaryfunc ) 0,  /* __int__ */
      ( unaryfunc ) 0,  /* __long__ */
      ( unaryfunc ) 0,  /* __float__ */
      ( unaryfunc ) 0,  /* __oct__ */
      ( unaryfunc ) 0,  /* __hex__ */
};

PyTypeObject matrix_Type = {
      PyObject_HEAD_INIT( NULL )    /*   required python macro   */
      0,    /*ob_size */
      "Matrix",         /*tp_name */
      sizeof( MatrixObject ), /*tp_basicsize */
      0,                /*tp_itemsize */
      ( destructor ) Matrix_dealloc,      /*tp_dealloc */
      ( printfunc ) 0,  /*tp_print */
      ( getattrfunc ) Matrix_getattr,     /*tp_getattr */
      ( setattrfunc ) Matrix_setattr,     /*tp_setattr */
      0,                /*tp_compare */
      ( reprfunc ) Matrix_repr,     /*tp_repr */
      &Matrix_NumMethods,     /*tp_as_number */
      &Matrix_SeqMethods,     /*tp_as_sequence */
};

//******************************************************************
//Function:                              newMatrixObject
//******************************************************************
PyObject *newMatrixObject( float *mat, int rowSize, int colSize )
{
      MatrixObject *self;
      int row, col, x;

      if( rowSize < 2 || rowSize > 4 || colSize < 2 || colSize > 4 )
            return ( EXPP_ReturnPyObjError( PyExc_RuntimeError,
                                    "row and column sizes must be between 2 and 4\n" ) );

      self = PyObject_NEW( MatrixObject, &matrix_Type );

      //generate contigous memory space
      self->contigPtr = PyMem_Malloc( rowSize * colSize * sizeof( float ) );
      if( self->contigPtr == NULL ) {
            return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                    "problem allocating array space\n\n" ) );
      }
      //create pointer array
      self->matrix = PyMem_Malloc( rowSize * sizeof( float * ) );
      if( self->matrix == NULL ) {
            PyMem_Free( self->contigPtr );
            return ( EXPP_ReturnPyObjError( PyExc_MemoryError,
                                    "problem allocating pointer space\n\n" ) );
      }
      //pointer array points to contigous memory
      for( x = 0; x < rowSize; x++ ) {
            self->matrix[x] = self->contigPtr + ( x * colSize );
      }

      if( mat ) {       //if a float array passed
            for( row = 0; row < rowSize; row++ ) {
                  for( col = 0; col < colSize; col++ ) {
                        self->matrix[row][col] =
                              mat[( row * colSize ) + col];
                  }
            }
      } else {          //or if NULL passed
            for( row = 0; row < rowSize; row++ ) {
                  for( col = 0; col < colSize; col++ ) {
                        self->matrix[row][col] = 0.0f;
                  }
            }
      }

      //set size vars of matrix
      self->rowSize = rowSize;
      self->colSize = colSize;

      //set coercion flag
      self->flag = 0;

      return ( ( PyObject * ) self );
}

Generated by  Doxygen 1.6.0   Back to index