// Matrix_Sparse.h Copyright (C) 2004 adolfo@di-mare.com #ifdef Spanish_dox /** \file Matrix_Sparse.h \brief Declaraciones y definiciones para la clase \c Matrix_Sparse<>. \author Adolfo Di Mare \date 2009 - Why English names??? ==> http://www.di-mare.com/adolfo/binder/c01.htm#sc04 */ #endif #ifdef English_dox /** \file Matrix_BASE.h \brief Declarations and definitiones for class \c Matrix_Sparse<>. \author Adolfo Di Mare \date 2009 - Why English names??? ==> http://www.di-mare.com/adolfo/binder/c01.htm#sc04 */ #endif #ifndef Matrix_Sparse_h #define Matrix_Sparse_h #include "Matrix_BASE.h" #include // assert() namespace Mx { #ifdef Spanish_dox #endif #ifdef English_dox #endif #ifdef Spanish_dox /// Matriz chirrisquitica almacenada como una matriz rala. #endif #ifdef English_dox /// Chirrisquitica matrix stored as a sparse matrix. #endif /// \copydetails Matrix_BASE template class Matrix_Sparse : public Matrix_BASE { public: Matrix_Sparse(unsigned m = 1, unsigned n = 1); Matrix_Sparse(const Matrix_Sparse& o); /// \copydoc Matrix_BASE::Matrix_BASE(const E&) Matrix_Sparse(const E& V) : m_capacity(0) { reSize(1,1); (*this)(0,0) = V; } ~Matrix_Sparse(); ///< Destructor. public: unsigned rows() const { return m_rows; } unsigned cols() const { return m_cols; } unsigned size() const { return size_Matrix( *this ); } unsigned count() const { return count_Matrix( *this ); } unsigned capacity() const { return m_capacity; } public: Matrix_Sparse& operator= (const Matrix_Sparse& M) { return copy(M); } Matrix_Sparse& copy( const Matrix_Sparse& M ); Matrix_Sparse& move( Matrix_Sparse& M ); Matrix_Sparse& swap( Matrix_Sparse& M ); void clear() { return clear_Matrix( *this ); } public: bool equals( const Matrix_Sparse & M ) const { return equals_Matrix( *this , M ); } public: E& operator()(unsigned i, unsigned j); const E& operator()(unsigned i, unsigned j) const; E& at(unsigned i, unsigned j) // throw (out_of_range); { return at_Matrix( *this , i,j ); } const E& at(unsigned i, unsigned j) const { return at_Matrix( *this , i,j ); } void reShape(unsigned m, unsigned n); void reSize( unsigned m, unsigned n); void transpose(); void setDefault(const E& same); const E& getDefault() { return m_same; } public: template friend bool check_ok( const Matrix_Sparse& M ); template friend class test_Matrix; ///< BUnit Test. private: /// Ajusta la matriz para que pueda almacenar \c n valores diferentes a \c getDefault(). void reserve(unsigned _Count); void reSize(unsigned newsize); private: unsigned* m_I; ///< Indice "i" de \c M(i,j) 0 <= i < m_capacity unsigned* m_J; ///< Indice "j" de \c M(i,j) 0 <= i < m_capacity E * m_val; ///< Valor para \c M(i,j) 0 <= i < m_capacity unsigned m_size; ///< Cantidad de valores insertados en los 3 vectores paralelos unsigned m_capacity; ///< Tamaño de los 3 vectores paralelos unsigned m_rows; ///< Cantidad de filas de la matriz unsigned m_cols; ///< Cantidad de columnas de la matris E m_same; ///< Valor almacenado en la mayor parte de la \c Matrix_Sparse private: template friend void add_Matrix( Matrix_Sparse& Res, const Matrix_Sparse& M ); }; // Matrix_Sparse /** Verifica la invariante de la clase. - El campo \c m_same indica cuál es el valor que se repite más en toda la matriz. - Usualmente \c same es el neutro aditivo \c value_type(). - No existe un constructor explícito para darle a \c m_same su valor inicial, que es siempre inicializado en \c value_type(). Para cambiarlo es necesario invocar el método \c setgetDefault(). - Los vectores \c m_I[], \c m_J[] y \c m_val[] son vectores paralelos, todos de longitud \c Matrix_Sparse::m_capacity. - La cantidad máxima de valores diferente que pueden ser almacenados en la matriz es \c Matrix_Sparse::m_capacity. - El largo de estos vectores aumenta cuando se hace referencia a un valor M(i,j) mediante la versión que no es \c const del \c operator()(i,j). O sea, que es ese método el encargado de agregar valores en \c m_val[], pues el operador homónimo \c const operator()(i,j) nunca agrega nada y, como es \c const, en general retorna una referencia constante a \c m_same. - 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_capacity. - 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 ____________________________________ / m_capacity \ +---+---+---+---+---+---+-..-+---+---+ 0 1 2 3 4 5 6 m_I--->| 1 | 3 | 3 |'-'| ... ... |'-'| 0 / - - - - - - - \ +---+---+---+---+ ... ... +---+ 1 | - a - - - - - | m_J--->| 1 | 2 | 1 |'-'| ... ... |'-'| 2 | - - - - - - - | +---+---+---+---+ ... ... +---+ 3 | - c b - - - - | m_val-->|'a'|'b'|'c'|'-'| ... ... |'-'| 4 \ - - - - - - - / +---+---+---+---+---+---+-..-+---+---+ 0 1 2 | m_size--------+ == 3 m_same == '-' m_rows == 5 m_cols == 7 \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_Sparse& M ) { if ( M.m_capacity == 0 ) { // m_capacity es el "marcador" que determina el valor de los demás if ( M.m_I != 0 ) { return false; /// - check_ok(): (m_capacity == 0) <==> (m_I == 0) } if ( M.m_J != 0 ) { return false; /// - check_ok(): (m_capacity == 0) <==> (m_J == 0) } if ( M.m_val != 0 ) { return false; /// - check_ok(): (m_capacity == 0) <==> (m_val == 0) } } else { if ( M.m_rows == 0 ) { return false; /// - check_ok(): (m_rows == 0) ==> (m_capacity == 0) } if ( M.m_cols == 0 ) { return false; /// - check_ok(): (m_cols == 0) ==> (m_capacity == 0) } if ( M.m_I == 0 ) { return false; /// - check_ok(): (m_capacity != 0) <==> (m_I != 0) } if ( M.m_J == 0 ) { return false; /// - check_ok(): (m_capacity != 0) <==> (m_J != 0) } if ( M.m_val == 0 ) { return false; /// - check_ok(): (m_capacity != 0) <==> (m_val != 0) } } if (M.m_rows == 0 || M.m_cols == 0) { if (M.m_rows == 0 && M.m_cols == 0) { // OK ==> los 2 son nulos } else { return false; /// - check_ok(): (m_rows == 0) <==> (m_cols == 0) } } if ( M.m_size > M.m_capacity ) { return false; /// - check_ok(): ( m_size <= m_capacity ) } if ( ! check_ok (M.m_same) ) { return false; /// - check_ok(): check_ok (m_same) } for (unsigned k=0; k( (0<=m_I[k]) && (m_I[k] < m_rows) ) k = [0..m_size-1] } if ( ! ( (0<=M.m_J[k]) && (M.m_J[k]( (0<=m_J[k]) && (m_J[k] < m_cols) ) k = [0..m_size-1] } if ( ! check_ok( M.m_val[k] ) ) { return false; /// - check_ok(): check_ok( m_val[k] ) } } return true; } /// Define el escalar que por defecto está en todas las entradas de la \c Matrix_Sparse. /// - Si same != getDefault() la matriz queda vacía. /// - De lo contrario, nada ocurre. template inline void Matrix_Sparse::setDefault(const E& same) { if ( m_same != same ) { m_same = same; m_size = 0; } } /** 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_Sparse::Matrix_Sparse(unsigned m, unsigned n) : m_same() { if (m == 0 || n == 0) { m_rows = m_cols = 0; m_val = 0; m_capacity = 0; } else { m_rows = m; m_cols = n; m_capacity = 16; // pues 16 == 2*2*2*2 ... m_I = new unsigned [m_capacity]; m_J = new unsigned [m_capacity]; m_val = new E[m_capacity]; } m_size = 0; // m_same = E(); // Usa el número "cero" como neutro tipo "E" } template inline Matrix_Sparse::Matrix_Sparse(const Matrix_Sparse& o) { m_capacity = 0; // m_capacity==0 indica que NO se está usando memoria dinámica copy( o ); // assert( this != &o ); return; // ... pues ya "o" existe y se está usando para incializar a "*this" /* NOTA Cuando un objeto usa memoria dinámica, copiarlo es un proceso diferente a inicializarlo a partir de otro valor de su misma clase. En otras palabras, el trabajo que hace el constructor CLS::CLS( const CLS& ) es muy diferente al que hace el copiador CLS& operator=( const CLS & ). El truco de programación usado en en esta implementación consiste en "marcar" el valor "m_capacity" para saber si se ha obtenido o no memoria dinámica para los vectores paralelos. La implementación de copy() está hecha de manera que si "m_capacity" es un puntero nulo, en copy() se inicializan correctamente todos los del Rep de Matrix_Sparse. Eso explica por qué la implementación de copy() es un tanto más elaborada que lo que a primera vista pareciera que es necesario. */ } template inline Matrix_Sparse::~Matrix_Sparse() { if (m_capacity != 0) { // sólo devuelve la memoria dinámica que sí ha adquirido delete [] m_I; delete [] m_J; // Retorna la memoria dinámica en uso delete [] m_val; } } /** 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_Sparse& Matrix_Sparse::copy( const Matrix_Sparse &o ) { if (this == &o) { // evita auto-borrado return *this; } if (o.m_capacity == 0) { if (m_capacity != 0) { // Como copy() puede ser invocado cuando el Rep delete [] m_val; // NO ha sido inicializado, en esta implementación delete [] m_J; // hay que tener cuidado de no usar los otros campos delete [] m_I; // del Rep antes de darles valor. } // - "m_capacity" SI debe tener su valor correcto. m_rows = m_cols = 0; m_val = 0; m_same = o.m_same; m_size = 0; // assert( o.m_size == 0 ); m_capacity = 0; // assert( o.m_capacity == 0 ); return *this; } // se asegura de que las dos matrices sean del mismo tamaño if (m_capacity != o.m_capacity) { // truco para evitar borrar la memoria dinámica if (m_capacity != 0) { // sólo devuelve la memoria dinámica que sí ha adquirido delete [] m_val; delete [] m_J; delete [] m_I; } m_capacity = o.m_capacity; m_I = new unsigned [m_capacity]; m_J = new unsigned [m_capacity]; m_val = new typename Matrix_BASE::value_type[m_capacity]; } m_same = o.m_same; m_size = o.m_size; // como las matrices son del mismo tamaño, // basta copiarlas entrada por entrada. m_rows = o.m_rows; m_cols = o.m_cols; for (unsigned k=0; k (*this == o) \par Complejidad: O( rows() * cols() ) \returns \c *this \see http://www.di-mare.com/adolfo/binder/c04.htm#sc07 */ template Matrix_Sparse& Matrix_Sparse::move( Matrix_Sparse &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; m_rows = o.m_rows; m_cols = o.m_cols; o.m_val = 0; o.m_rows = o.m_cols = 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_Sparse::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 Matrix_Sparse& Matrix_Sparse::swap( Matrix_Sparse &o ) { std::swap( this->m_I , o.m_I ); std::swap( this->m_J , o.m_J ); std::swap( this->m_val , o.m_val ); std::swap( this->m_size , o.m_size ); std::swap( this->m_capacity , o.m_capacity ); std::swap( this->m_rows , o.m_rows ); std::swap( this->m_cols , o.m_cols ); std::swap( this->m_same , o.m_same ); 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_Sparse::value_type(). \pre (m > 0) && (n > 0) */ template void Matrix_Sparse::reSize(unsigned m, unsigned n) { // NO hace falta cambiar m_capacity porque lo único que está cambiando // el la dimensión de la matriz if (m != m_rows || n != m_cols) { m_size = 0; // desecha todos los valores anteriores m_rows = m; // cambia las dimensiones de la matriz m_cols = n; } if (m==0 || n==0) { m_rows = 0; m_cols = 0; m_size = 0; if (m_capacity != 0) { m_capacity = 0; // todo esto mantiene la invariante de la clase delete [] m_I; m_I = 0; delete [] m_J; m_J = 0; delete [] m_val; m_val = 0; } } 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 == Matrix_Sparse::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. */ /* NOTA Esta es la imnplementación anterior, que debe reacomodar los ínices (i,j) de aquellos valores almacenados en los vectores paralalos para simular que la matriz rala en realidad está almacenada como una matriz densa, dentro de un vector. - Al modificar la especificación de para reSize() ya no hace falta hacer todo este trabajo. NO hace falta cambiar m_capacity porque lo único que está cambiando el la dimensión de la matriz */ /* Nota para esta implementación: Este método funciona de una forma similar a reSize() para la matriz implementada con un vector. Para lograr que las clases Matrix y Matrix_Sparse tengan la misma funcionalidad, esta implementación reacomoda los índices de la matriz de manera que las m_rows*m_cols entradas queden reacomodadas como si la implementación usara un vector de dimension m*n. Sin embargo, lo lógico sería eliminar todos los valores de la matriz trocando el ciclo "for(;k;){}" por un eficiente "m_size=0;" - Esto muestra que una especificación que aparenta ser "amplia", como lo es la de Matrix::reSize(n,m), en realidad se queda "corta" si se toma en cuenta que la misma especificación debe servir para varias implementaciones diferentes. - Además, muestra que algunas operaciones sólo tienen sentido para una implementación específica, como ocurre con las operaciones Matrix_Sparse>::getDefault() y su "seteador" Matrix_Sparse>::setgetDefault(). */ #if 0 if (m * n == m_rows * m_cols) { m_size = 0; // desecha todos los valores anteriores } else { unsigned row = 0, col = 0; for (unsigned k=0; k m_val[ (i * m_cols) + j ] ; // si almacena los valores por filas // M(i,j) <==> m_val[ i + (j * m_rows) ] ; // si almacena los valores por columnas unsigned lineal = (m_I[k] * m_cols) + m_J[k]; // por filas // unsigned lineal = m_I[k] + (m_j[k] * m_rows); // por columnas m_I[k] = lineal / n; // lineal % m; m_J[k] = lineal % n; // lineal / m; } } m_rows = m; // cambia las dimensiones de la matriz m_cols = n; #endif } /** 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_Sparse::reShape(unsigned m, unsigned n) { if ( m * n == m_rows * m_cols ) { reSize(n,m); } } /// 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_Sparse::operator()(unsigned i, unsigned j) const { assert( "Matrix_Sparse::operator()()" && (i < rows()) ); assert( "Matrix_Sparse::operator()()" && (j < cols()) ); for (unsigned k=0; kM(i,j) = val; // M(i,j) es un "lvalue" (modificable) template inline E& Matrix_Sparse::operator()(unsigned i, unsigned j) { assert( "Matrix_Sparse::operator()()" && (i < rows()) ); assert( "Matrix_Sparse::operator()()" && (j < cols()) ); // Busca al elemento (i,j) en los vectores paralelos for (unsigned k=0; k Agotado m_val[] Matrix_Sparse::operator()()" ); if (m_size == m_capacity) { unsigned new_capacity = m_capacity; if (m_capacity % 16 == 0) { new_capacity += 16; } else { new_capacity /= 16; // Hace que el nuevo tamaño sea divisible por 16 new_capacity *= 16; new_capacity += 16; } reSize( new_capacity ); // Agrega 16 porque 16 == 2*2*2*2 ... } // Agrega un nuevo valor modificable en los vectores paralelos m_I[m_size] = i; m_J[m_size] = j; m_val[m_size] = m_same; // m_last_change = & m_val[m_size]; return m_val[m_size++]; } /// Le cambia el tamaño máximo posible a la matriz. /// - Le aumenta la cantidad de valores diferentes a \c getDefault(). /// - No hace nada cuando size() < newsize. template void Matrix_Sparse::reSize(unsigned newsize) { if ( newsize < m_size ) { // hay más valores en uso que "newsize" // assert((m_size < newsize) && "=>Matrix_Sparse::reSize()" ); return; } unsigned * new_I = new unsigned [ newsize ]; // Obtiene los vectores nuevos unsigned * new_J = new unsigned [ newsize ]; E * new_val = new E [ newsize ]; if (m_capacity > 0) { // Copia los valores en uso for (unsigned k=0; k void add_Matrix( Matrix_Sparse& Res, const Matrix_Sparse & M ) { // verifica que las dos matrices sean del mismo tamaño assert( "Matrix_Sparse::add()" && ( Res.rows() == M.rows() ) ); assert( "Matrix_Sparse::add()" && ( Res.cols() == M.cols() ) ); // Como las matrices son del mismo tamaño basta sumarlas entrada por entrada. typename Matrix_BASE::value_type *pRes = Res.m_val; typename Matrix_BASE::value_type *pM = & M.m_val[0]; typename Matrix_BASE::value_type *pEND = & Res.m_val[M.m_cols * M.m_rows]; for ( ; pRes != pEND; ++pRes, ++pM ) { *pRes += *pM; } return; } #if 0 /** 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_Sparse operator-( Matrix_Sparse&, Matrix_Sparse ) - 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 esta función amiga que realiza el trabajo, aunque es poco cómoda de usar. */ /// \copydoc Matrix_BASE::substract() template void substract_Matrix( Matrix_Sparse& Res, const Matrix_Sparse & M ) { // verifica que las dos matrices sean del mismo tamaño assert( "Matrix_Sparse::substract()" && ( Res.rows() == M.rows() ) ); assert( "Matrix_Sparse::substract()" && ( Res.cols() == M.cols() ) ); // Como las matrices son del mismo tamaño basta restarlas entrada por entrada. typename Matrix_BASE::value_type *pRes = Res.m_val; typename Matrix_BASE::value_type *pM = & M.m_val[0]; typename Matrix_BASE::value_type *pEND = & Res.m_val[M.m_cols * M.m_rows]; for ( ; pRes != pEND; ++pRes, ++pM ) { *pRes += *pM; } return; } #endif } // namespace Mx #endif // Matrix_Sparse_h // EOF: Matrix_Sparse.h