// Matrix.h Copyright (C) 2004 adolfo@di-mare.com /** \file Matrix.h \brief Declaraciones y definiciones para la clase \c Matrix. \author Adolfo Di Mare \date 2004 - Why English names??? ==> http://www.di-mare.com/adolfo/binder/c01.htm#sc04 */ #ifndef Matrix_h #define Matrix_h // Macro definida para evitar que \c assert() cancele la ejecución. // Tira \c std::out_of_range si \c assert(X) falla. #ifdef ASSERT_NO_ABORT #include // std::out_of_range #include #include // cout #ifndef NDEBUG #ifdef assert #undef assert #endif #define assert(X) { \ std::string line, msg = "Assertion failed: " #X ", file " __FILE__ ; \ unsigned long n = __LINE__; \ if ( n==0 ) { line = "0"; } \ else do { \ char ch = (n%10) + '0'; \ line = ch + line; \ n = n/10; \ } while ( n!=0 ); \ msg += ", line " + line + "\n"; \ if ( !(X) ) { throw std::out_of_range( msg ); } \ } #endif #else #include // assert() #endif /// Definido por la biblioteca C++ estándar namespace std {} // Está acá para que Doxygen lo documente /// Matriz chirrisquitica de adolfo@di-mare.com namespace Mx { /** Esta es una clase matriz muy chirrisquitica que puede cambiar dinámicamente de tamaño. - La matriz tiene tamaño \c rows() x \c cols() - Se le puede cambiar el tamaño dinámicamente con el método \c reSize(). - La clase \c E debe incluir un neutro para la adición, cuyo valor debe poderse obtener invocando el convertidor \c Sparse_Matrix::value_type(). - Las operaciones aritméticas "+" "-" "*" deben estar definidas para \c Matrix::value_type y debe existir el valor \c Matrix::value_type() y también \c Matrix::value_type(1) (para matrices unitarias) - Los valores de la matriz pueden estar almacenados por filas o por columnas, según quede implementado el método \c Matrix::operator(unsigned, unsigned) \see http://www.oonumerics.org/oon/ */ template class Matrix { public: /// Tipo del objeto almacenado, similar al nombre usado en STL typedef E value_type; /// Tipo del objeto almacenado, similar al nombre usado en STL typedef value_type& reference; /// Tipo del objeto almacenado, similar al nombre usado en STL typedef const value_type& const_reference; /// Tipo del tamaño de un objeto, similar al nombre usado en STL typedef unsigned size_type; public: Matrix(unsigned m = 1, unsigned n = 1); Matrix(const Matrix& o); ///< Constructor de copia /// Matriz escalar de valor \c V. Matrix(const value_type V) : m_rows(0), m_cols(0), m_val(0) { reSize(1,1); (*this)(0,0) = V; } ~Matrix(); public: unsigned rows() const { return m_rows; } ///< Cantidad de filas de la matriz unsigned cols() const { return m_cols; } ///< Cantidad de columnas de la Matriz unsigned size() const { return m_rows * m_cols; } ///< Cantidad de valores almacenados en la matriz unsigned count() const { return size(); } ///< Cantidad de valores almacenados en la matriz /// Cantidad máxima posible de valores diferentes que pueden ser almacenados en la matriz size_type capacity() const { return size(); } public: Matrix& operator= (const Matrix &o) { return copy(o); } ///< Sinónimo de \c this->copy(o) Matrix& copy( const Matrix &o ); Matrix& move( Matrix &o ); Matrix& swap( Matrix &o ); public: /// ¿¿¿ (p == q) ??? friend bool operator == (const Matrix &p, const Matrix &q) { return p.equals(q); } /// ¿¿¿ (p != q) ??? friend bool operator != (const Matrix &p, const Matrix &q) { return !(p==q); } bool equals( const Matrix & o ) const; ///< ¿¿¿ (*this==o) ??? /// Retorna \c true si \c "o" comparte sus valores con \c "*this" bool same( const Matrix & o ) const { return this == &o; } protected: // disponibles para clases derivadas void add( const Matrix & ); void substract( const Matrix & ); void multiply( const Matrix &, const Matrix & ); public: friend Matrix operator + (const Matrix& A, const Matrix& B) { Matrix Res = A; Res.add(B); return Res; } ///< Retorna \c A+B friend Matrix operator - (const Matrix& A, const Matrix& B) { Matrix Res = A; Res.substract(B); return Res; } ///< Retorna \c A-B friend Matrix operator * (const Matrix& A, const Matrix& B) { Matrix Res; Res.multiply(A, B); return Res; } ///< Retorna \c A*B public: reference operator () (unsigned, unsigned); const_reference operator () (unsigned, unsigned) const; reference at (unsigned m, unsigned n) { return operator() (m,n); } ///< Retorna \c operator()(m,n). const_reference at (unsigned m, unsigned n) const { return operator() (m,n); } ///< Retorna \c operator()(m,n) \c "const". public: void reSize( unsigned, unsigned); void reShape(unsigned, unsigned); void transpose() { unsigned tmp = m_rows; m_rows = m_cols; m_cols = tmp; } ///< Transpone la matriz. template friend bool check_ok( const Matrix& M ); template friend class test_Matrix; ///< Datos de prueba para la clase private: value_type * m_val; ///< Vector de valores de la matriz unsigned m_rows; ///< Cantidad de filas de la matriz unsigned m_cols; ///< Cantidad de columnas de la matris }; // Matrix /** Verifica la invariante de la clase. - Es posible que la matriz tenga dimensiones nulas, lo que implica que todos los punteros a los vectors paralelos deben ser nulos. Este hecho se marca dándolo el valor \c 0 (cero) al campo \c m_val. - Las matrices quedan almacenadas en un vector de tamaño [M x N]. - En todos los algoritmos, "m" o "m_rows" es la cantidad de filas == \c rows() - En todos los algoritmos, "n" o "m_cols" es la cantidad de columnas == \c cols() \par Rep Modelo de la clase \code +---+ / \ | 2 | M(i,j) ==> m_val[ (i * m_cols) + j ] | 0 1 2 3 | m_rows == 2 +---+ (almacenamiento por filas) | 4 5 6 7 | m_cols == 4 | 4 | \ / +---+ +---+---+---+---+---+---+---+---+ | *-|-->| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | +---+ +---+---+---+---+---+---+---+---+ \endcode \par Rep Modelo de la clase \code +---+ | 4 | M(i,j) ==> m_val[ i + (j * m_rows) ] / a e \ +---+ (almacenamiento por columnas) | b f | m_rows == 4 | 2 | | c g | m_cols == 2 +---+ +---+---+---+---+---+---+---+---+ \ d h / | *-|-->| a | b | c | d | e | f | g | h | +---+ +---+---+---+---+---+---+---+---+ \endcode - http://www.di-mare.com/adolfo/binder/c03.htm#k1-Rep \remark Libera al programador de implementar el método \c Ok() - http://www.di-mare.com/adolfo/binder/c04.htm#sc11 */ template bool check_ok( const Matrix& M ) { if ( M.m_rows == 0 || M.m_cols == 0 ) { // si alguno es cero... if ( M.m_rows != 0 || M.m_cols != 0 ) { // ... los 2 deben serlo return false; /// - Invariante: (M.m_rows == 0) <==> (M.m_cols == 0) } else { if (M.m_val != 0) { return false; /// - Invariante: (M.m_rows == 0) <==> (M.m_val == 0) } return true; // La implementación sí permite (m_val == 0) } } unsigned MxN = M.m_rows * M.m_cols; for (unsigned k=0; kcheck_ok( m_val[k] ) } } return true; } /** \fn template inline Matrix::Matrix(const value_type V); \brief Constructor a partir de \c Matrix::value_type(V). - La matriz resultante es una matriz escalar, de dimensiones 1x1, y su valor es \c "V". */ /** Constructor de vector. - Obtiene suficiente memoria dinámica para almacenas los n * m valores de la matriz - Si \c "value_type" tiene un constructor de vector, lo usar para inicializar cada uno de los valores de la matriz; de lo contrario, los deja tal cual están en la memoria - Si \c "value_type" es uno de los tipos escalares básicos, como lo son \c int o \c float, los valores almacenados en la matriz quedan tal cual están y no son inicializados. \pre - m * n > 0 - (m > 0) && (n > 0) */ template inline Matrix::Matrix(unsigned m, unsigned n) { if (m == 0 || n == 0) { m_rows = m_cols = 0; m_val = 0; } else { m_rows = m; m_cols = n; m_val = new value_type [n*m]; } } template Matrix::Matrix(const Matrix& o) { if (o.m_val == 0) { // OJO: una matriz "vacía" produce errores en m_rows = m_cols = 0; // Matrix::operator(unsigned, unsigned) m_val = 0; return; } m_rows = o.m_rows; m_cols = o.m_cols; const unsigned MxN = o.m_rows * o.m_cols; m_val = new value_type [MxN]; // como las matrices son del mismo tamaño, // basta copiarlas entrada por entrada. value_type *pSelf = this->m_val; value_type *pO = o.m_val; value_type *pEND = &m_val[MxN]; for (; pSelf != pEND; ++pSelf, ++pO) { *pSelf = *pO; } } /// Destructor template inline Matrix::~Matrix() { if ( m_val != 0 ) { delete [] m_val; } } template bool Matrix::equals( const Matrix & o ) const { if (m_rows != o.m_rows || m_cols != o.m_cols) { return false; } else if ( m_val == 0 ) { return true; // las 2 matrices están vacías ==> son iguales } const unsigned MxN = o.m_rows * o.m_cols; value_type *pSelf = this->m_val; value_type *pO = o.m_val; value_type *pEND = &m_val[MxN]; for (; pSelf != pEND; ++pSelf, ++pO) { if (*pSelf != *pO) { return false; } } return true; } /** Copia desde \c "o". - Copia todo el valor de \c "o" sobre \c "*this", de forma que el nuevo valor de \c "*this" sea un duplicado exacto del valor de \c "o" - El valor anterior de \c "*this" se pierde - El objeto \c "o" mantiene su valor anterior - Luego de la copia, cuando el valor de \c "o" cambia, el de \c "*this" no cambiará, y viceversa, pues la copia es una copia profunda; no es superficial - Si \c "*this" es \c "o" entonces su valor no cambia - Después de esta operación siempre ocurre que *this == o \par Complejidad: O( rows() * cols() ) \returns *this \see http://www.di-mare.com/adolfo/binder/c04.htm#sc05 */ template Matrix& Matrix::copy( const Matrix &o ) { if (this == &o) { // evita auto-borrado return *this; } else if (o.m_val == 0) { if (m_val != 0) { delete [] m_val; m_rows = m_cols = 0; m_val = 0; } return *this; } // se asegura de que las dos matrices son del mismo tamaño const unsigned MxN = o.m_rows * o.m_cols; if (MxN != m_rows * m_cols) { // truco para evitar borrar la memoria dinámica delete [] m_val; m_val = new value_type [MxN]; } m_rows = o.m_rows; m_cols = o.m_cols; // como las matrices son del mismo tamaño, // basta copiarlas entrada por entrada. value_type *pSelf = this->m_val; value_type *pO = o.m_val; value_type *pEND = &m_val[MxN]; for (; pSelf != pEND; ++pSelf, ++pO) { *pSelf = *pO; } return *this; } /** Traslada el valor de \c "o" a \c "*this". - El valor anterior de \c "*this" se pierde - El nuevo valor de \c "*this" es el que \c "o" tuvo - \c "o" queda en el estado en que lo dejaría \c Erase() - Si \c "*this" es \c "o" entonces su valor no cambia - En general, después de esta operación casi nunca ocurre que (*this == o) \par Complejidad: O( rows() * cols() ) \returns \c *this \see http://www.di-mare.com/adolfo/binder/c04.htm#sc07 */ template Matrix& Matrix::move( Matrix &o ) { if (this == &o) { // evita auto-borrado return *this; } else if (o.m_val == 0) { if (m_val != 0) { delete [] m_val; m_rows = m_cols = 0; m_val = 0; } return *this; } else if ( m_val != 0 ) { delete [] m_val; } m_val = o.m_val; // me robo los valores de "o" m_rows = o.m_rows; m_cols = o.m_cols; o.m_rows = o.m_cols = 0; o.m_val = 0; return *this; } /** Intercambia los valores de \c "*this" y \c "o". - Debido a que algunos métodos retornan copias del valor de \c "*this" en lugar de una referencia, como ocurre con \c Matrix::Child(), a veces \c swap() no tiene el resultado esperado por el programador. - Por ejemplo, si se invoca T.Child(i). swap( T.Child(j) ) el resultado no es intercambiar los hijos, sino más bien intercambiar los valores de los sub-árboles temporales \c T.Child(i) y \c T.Child(j). La forma correcta de intercambiar hijos es usar \c Graft(). \par Complejidad: O( \c 1 ) \returns *this \see http://www.di-mare.com/adolfo/binder/c04.htm#sc08 */ template inline Matrix& Matrix::swap( Matrix &o ) { std::swap( this->m_val , o.m_val ); std::swap( this->m_rows , o.m_rows ); std::swap( this->m_cols , o.m_cols ); return *this; } /** Le cambia las dimensiones a la matriz. - En algunos casos los valores almacenados en la matriz no quedan todos iguales a \c Matrix::value_type(). - Si (m * n == 0) deja la matriz vacía. */ template void Matrix::reSize(unsigned m, unsigned n) { const unsigned MxN = m * n; if (MxN == 0) { if (m_val != 0) { // sólo borra si hay algo que borrar delete [] m_val; m_rows = m_cols = 0; m_val = 0; } } else if ( MxN == m_rows * m_cols ) { // truco para evitar borrar la memoria dinámica m_rows = m; m_cols = n; } else { if (m_val != 0) { delete [] m_val; } m_val = new value_type [MxN] ; m_rows = m; m_cols = n; } return; /* NOTA Esta es la antigua especificación de reSize(). Es incorrecta porque presume que el Rep de la matriz es un vector denso en donde están almacenados todos los valores de la matriz: +----------------------------------------------------------+ | reSize(): Le cambia las dimensiones a la matriz. | | ======== | | - Si ocurre que (m*n) == rows() * cols() los valores de | | la matriz se mantienen, aunque cambian sus dimensiones.| | - Si (m*n) != rows() * cols() los valores de la matriz | | quedarán inicializados de la misma forma en que los | | inicializaría CERO == Sparse_Matrix::value_type(). | | | | \pre (m > 0) && (n > 0) | +----------------------------------------------------------+ En la actual especificación, que ya está corregida, no queda implícita la presunción sobre cómo está organizada internamente la matriz. Por eso, esta nueva especificación sí sirve para una matriz implementada con un vector denso de valores, o para la matriz implementada como una matriz rala. Estos pequeños detalles en algunas ocasiones surgen cuando el programador de una clase introduce mejoras o modificaciones, pues muchas veces es muy difícil o prácticamente imposible predecir todos los pormenores y detalles de una especificación o de una implementación. */ } /** Le ajusta las dimensiones a la matriz. - Si ocurre que (m*n) == rows()*cols() hace lo mismo que haría \c reSize(m,n). - En caso contrario, no hace nada. */ template inline void Matrix::reShape(unsigned m, unsigned n) { const unsigned MxN = m * n; if (MxN == m_rows * m_cols) { m_rows = m; m_cols = n; } } /// Retorna una referencia al elemento [i,j] de la matriz ( \c const ). /// - \c M(i,j) significa lo que en arreglos se denota con \c M[i][j]. /// - val = M(i,j); // M(i,j) es un "rvalue" (const) template inline const E& Matrix::operator() (unsigned i, unsigned j) const { assert( "Matrix::operator()()" && (i < rows()) ); assert( "Matrix::operator()()" && (j < cols()) ); return m_val[ (i * m_cols) + j ] ; // almacena los valores por filas // return m_val[ i + (j * m_rows) ] ; // almacena los valores por columnas } /// Retorna una referencia al elemento [i,j] de la matriz. /// - \c M(i,j) significa lo que en arreglos se denota con \c M[i][j]. /// - M(i,j) = val; // M(i,j) es un "lvalue" (modificable) template inline E& Matrix::operator() (unsigned i, unsigned j) { assert( "Matrix::operator()()" && (i < rows()) ); assert( "Matrix::operator()()" && (j < cols()) ); return m_val[ (i * m_cols) + j ] ; // almacena los valores por filas // return m_val[ i + (j * m_rows) ] ; // almacena los valores por columnas } /** Le suma a \c "*this" la matriz \c "O". \pre - \c "*this" y \c "O" deben tener las mismas dimensiones - rows() == O.rows() && cols() == O.cols() \par Complejidad: O( rows() * cols() ) \remarks - Esta es la implementación de Matrix operator+( Matrix&, Matrix ) - El compilador tiene problemas en compilar un función amiga ("friend") definida con plantillas si esa función amiga no está definida (implementada) dentro de la declaración de la clase. Para solventar esta deficiencia existe este método que realiza el trabajo, aunque es poco cómodo de usar. */ template void Matrix::add( const Matrix & O ) { // verifica que las dos matrices sean del mismo tamaño assert( (rows() == O.rows()) && (cols() == O.cols()) ); // Como las matrices son del mismo tamaño basta sumarlas entrada por entrada. value_type *pThis = m_val; value_type *pO = & O.m_val[0]; value_type *pEND = &m_val[m_cols * m_rows]; for ( ; pThis != pEND; ++pThis, ++pO) { *pThis += *pO; } return; } /** Le resta a \c "*this" la matriz \c "O". \pre - \c "*this" y \c "O" deben tener las mismas dimensiones - rows() == O.rows() && cols() == O.cols() \par Complejidad: O( rows() * cols() ) \remarks - Esta es la implementación de Matrix operator-( Matrix&, Matrix ) - El compilador tiene problemas en compilar un función amiga ("friend") definida con plantillas si esa función amiga no está definida (implementada) dentro de la declaración de la clase. Para solventar esta deficiencia existe este método que realiza el trabajo, aunque es poco cómodo de usar. */ template void Matrix::substract( const Matrix & O ) { // verifica que las dos matrices sean del mismo tamaño assert( (rows() == O.rows()) && (cols() == O.cols()) ); // Como las matrices son del mismo tamaño basta sumarlas entrada por entrada. value_type *pThis = m_val; value_type *pO = & O.m_val[0]; value_type *pEND = &m_val[m_cols * m_rows]; for ( ; pThis != pEND; ++pThis, ++pO) { *pThis -= *pO; } return; } /** Calcula la multiplicación A * B y la almacena en \c "*this". - Las dimensiones de \c "*this" se ajustan de manera que: - rows() == A.rows() && cols() == B.cols() \pre - \c "A" y \c "B" deben tener dimensiones compatibles - A.cols() == B.rows() - La multiplicación se hace [Fila x Columna], lo que implica que la cantidad de valores en la filas de \c "A" debe ser igual a la cantidad de columnas de \c "B" \par Complejidad: O( A.cols() * B.cols() * A.cols() ) \remarks - Esta es la implementación de Matrix operator*( Matrix&, Matrix ) - El compilador tiene problemas en compilar un función amiga ("friend") definida con plantillas si esa función amiga no está definida (implementada) dentro de la declaración de la clase. Para solventar esta deficiencia existe este método que realiza el trabajo, aunque es poco cómodo de usar. */ template void Matrix::multiply( const Matrix & A, const Matrix & B ) { // Verifica que las matrices se puedan multiplicar assert( (A.cols() == B.rows()) && " => Matrix::multiply()" ); reSize( A.rows(), B.cols() ); value_type sum; for (unsigned i=0; i(i,j) = sum; // produce un error de compilación // this->operator()(i,j) = sum; // también funciona (*this)(i,j) = sum; // también funciona } } return; } /// Graba en el flujo \c COUT el valor de \c M[][]. template std::ostream& operator<<(std::ostream& COUT, const Matrix& M) { COUT << '[' << M.rows() << 'x' << M.cols() << ']' << std::endl; for (unsigned i=0; i < M.rows(); ++i) { for (unsigned j=0; j < M.cols(); ++j) { COUT << " " << M(i,j); } COUT << std::endl; } return COUT; } /// Obtiene del flujo \c CIN el valor para \c M[][]. template std::istream& operator>>(std::istream& CIN, Matrix& M) { assert( "This code has not been tested" ); unsigned rows,cols; CIN >> rows >> cols; M.reSize(rows,cols); for (unsigned i=0; i> M(i,j); } } return CIN; } } // namespace Mx #include "Matrix_Lib.h" #endif // Matrix_h // EOF: Matrix.h