Matrix_Sparse.h

Go to the documentation of this file.

// 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 <adolfo@di-mare.com>
    \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 <adolfo@di-mare.com>
    \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 <cassert> // 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 E>
class Matrix_Sparse : public Matrix_BASE<E> {
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 <class T> friend bool  check_ok( const Matrix_Sparse<T>& M );
    template <class T> 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) <code>0 <= i < m_capacity</code>
    unsigned* m_J;   ///< Indice "j" de \c M(i,j) <code>0 <= i < m_capacity</code>
    E       * m_val; ///< Valor para \c M(i,j) <code>0 <= i < m_capacity</code>
    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 <class T> friend
    void add_Matrix( Matrix_Sparse<T>& Res, const Matrix_Sparse<T>& 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 <em>Rep</em> 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 <class T>
bool check_ok( const Matrix_Sparse<T>& 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(): <code>(m_capacity == 0) <==> (m_I == 0)</code>
        }
        if ( M.m_J != 0 ) {
            return false; /// - check_ok(): <code>(m_capacity == 0) <==> (m_J == 0)</code>
        }
        if ( M.m_val != 0 ) {
            return false; /// - check_ok(): <code>(m_capacity == 0) <==> (m_val == 0)</code>
        }
    }
    else {
        if ( M.m_rows == 0 ) {
            return false; /// - check_ok(): <code>(m_rows == 0) ==> (m_capacity == 0)</code>
        }
        if ( M.m_cols == 0 ) {
            return false; /// - check_ok(): <code>(m_cols == 0) ==> (m_capacity == 0)</code>
        }
        if ( M.m_I == 0 ) {
            return false; /// - check_ok(): <code>(m_capacity != 0) <==> (m_I != 0)</code>
        }
        if ( M.m_J == 0 ) {
            return false; /// - check_ok(): <code>(m_capacity != 0) <==> (m_J != 0)</code>
        }
        if ( M.m_val == 0 ) {
            return false; /// - check_ok(): <code>(m_capacity != 0) <==> (m_val != 0)</code>
        }
    }

    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(): <code>(m_rows == 0) <==> (m_cols == 0)</code>
        }
    }

    if ( M.m_size > M.m_capacity ) {
        return false; /// - check_ok(): <code>( m_size <= m_capacity )</code>
    }
    if ( ! check_ok (M.m_same) ) {
        return false; /// - check_ok(): <code>check_ok (m_same)</code>
    }
    for (unsigned k=0; k<M.m_size; ++k) {
        if ( ! ( (0<=M.m_I[k]) && (M.m_I[k]<M.m_rows) ) ) {
            return false; /// - check_ok(): <code>( (0<=m_I[k]) && (m_I[k] < m_rows) ) k = [0..m_size-1]</code>
        }
        if ( ! ( (0<=M.m_J[k]) && (M.m_J[k]<M.m_cols) ) ) {
            return false; /// - check_ok(): <code>( (0<=m_J[k]) && (m_J[k] < m_cols) ) k = [0..m_size-1]</code>
        }
        if ( ! check_ok( M.m_val[k] ) ) {
            return false; /// - check_ok(): <code>check_ok( m_val[k] )</code>
        }
    }
    return true;
}

/// Define el escalar que por defecto está en todas las entradas de la \c Matrix_Sparse.
/// - Si <code>same != getDefault()</code> la matriz queda vacía.
/// - De lo contrario, nada ocurre.
template<class E>
inline void Matrix_Sparse<E>::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
      <code> n * m </code> 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
     -  <code> m * n > 0  </code>
     -  <code> (m > 0) && (n > 0) </code>
*/
template <class E>
inline Matrix_Sparse<E>::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 <class E>
inline Matrix_Sparse<E>::Matrix_Sparse(const Matrix_Sparse<E>& 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 <class E>
inline Matrix_Sparse<E>::~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 <code> *this == o </code>.

    \par Complejidad:
         O( <code>  rows() * cols() </code> )

    \returns *this
    \see http://www.di-mare.com/adolfo/binder/c04.htm#sc05
*/
template <class E>
Matrix_Sparse<E>& Matrix_Sparse<E>::copy( const Matrix_Sparse<E> &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<E>::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<m_size; ++k) {
        m_I[k]   = o.m_I[k];
        m_J[k]   = o.m_J[k];
        m_val[k] = o.m_val[k];
    }
    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 <code> (*this == o) </code>

    \par Complejidad:
         O( <code>  rows() * cols() </code> )

    \returns \c *this

    \see http://www.di-mare.com/adolfo/binder/c04.htm#sc07
*/
template <class E>
Matrix_Sparse<E>& Matrix_Sparse<E>::move( Matrix_Sparse<E> &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 <code> T.Child(i). swap( T.Child(j) ) </code>
      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 <class E>
Matrix_Sparse<E>& Matrix_Sparse<E>::swap( Matrix_Sparse<E> &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<E>::value_type().
    \pre <code> (m > 0) && (n > 0) </code>
*/
template <class E>
void Matrix_Sparse<E>::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<E>::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<E> y Matrix_Sparse<E> 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<E>::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<E>>::getDefault() y su
      "seteador" Matrix_Sparse<E>>::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_size; ++k) {
            // M(i,j) <==> 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 <code> (m*n) == rows()*cols() </code> hace
      lo mismo que haría \c reSize(m,n).
    - En caso contrario, no hace nada.
*/
template <class E>
inline void Matrix_Sparse<E>::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].
/// - <code>val = M(i,j); // M(i,j) es un "rvalue" (const)</code>
template <class E>
inline const E& Matrix_Sparse<E>::operator()(unsigned i, unsigned j) const {
    assert( "Matrix_Sparse<E>::operator()()" && (i < rows()) );
    assert( "Matrix_Sparse<E>::operator()()" && (j < cols()) );

    for (unsigned k=0; k<m_size; k++) {
        if ( m_I[k] == i ) {
            if ( m_J[k] == j ) {
                return m_val[k];
            }
        }
    }
    return m_same;
/*  NOTA
    Como este método es "const", de antemano se sabe que el programador no puede
    usarlo para modificar el valor de la matriz. Por eso, aunque el valor
    retornado sea una referencia a valor común por defecto m_same, de antemano
    el compilador asegura que ese valor no será modificado.

    Sin embargo, cuando el programador usa el método homónimo operator()(i,j)
    que no es "const", es posible que el valor retornado sí sea modificado.
    En ese caso, ya no es correcto retornar una referencia al valor común "m_same".
    - Por eso, cuando se usa el hace referencia en el otro operator()(i,j) es
      necesario agregar una entrada en los vectores paralelos en aquellos casos
      en que no existe un valor diferente a "m_same" para (i,j).
    - Esto quiere decir que sí es posible que al utilizar la versión modificable
      de operator()(i,j) quede en el vector "m_val[]" un valor que es igual a
      "m_same". En el caso peor, podría ocurrir que todos los valores almacenados
      en el vector "m_val[]" sean todos iguales a "m_same".
    - Una forma de corregir esta anomalía es "revisar después" si existe un valor
      en el vector "m_val[]" que es igual a "m_same". Una forma eficiente de
      hacerlo es mantener el el Rep un puntero "m_last_change" que apunte al
      último valor "m_val[]" que la versión modificable de operator()(i,j) retornó.
      En la siguiente invocación a operator()(i,j), se puede verificar si hubo un
      ese valor anterior fue cambiado de manera que tiene almacenado "m_same".
*/
#if 0
    // Este código podría ser agregado al principio de operator()(i,j)
    if (m_last_change != 0) {
        if ( *m_last_change == m_same ) { // el último uso de op()(i,j) dejó "m_same"
            if ( m_size == 1 ) { // Elimina los vectores viejos
                delete [] m_I;    m_I = 0;
                delete [] m_J;    m_J = 0;
                delete [] m_val;  m_val = 0;
                m_capacity = 0;
                m_size     = 0;
            }
            else { // intercambia este valor con el último
                m_size--;
                unsigned k = (m_last_change - m_val); // índice de "* m_last_change"
                *m_last_change = m_val[m_size];
                m_I[k] = m_I[m_size];
                m_J[k] = m_J[m_size];
            }
        }
        m_last_change = 0;
    }
#endif
}
/// 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].
/// - <code>M(i,j) = val; // M(i,j) es un "lvalue" (modificable)</code>
template <class E>
inline E& Matrix_Sparse<E>::operator()(unsigned i, unsigned j) {
    assert( "Matrix_Sparse<E>::operator()()" && (i < rows()) );
    assert( "Matrix_Sparse<E>::operator()()" && (j < cols()) );

    // Busca al elemento (i,j) en los vectores paralelos
    for (unsigned k=0; k<m_size; k++) {
        if ( m_I[k] == i ) {
            if ( m_J[k] == j ) {
            //  m_last_change = & m_val[k];
                return m_val[k];
            }
        }
    }

    assert( (m_size <= m_capacity) && " => Agotado m_val[] Matrix_Sparse<E>::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 <code>size() < newsize</code>.
template <class E>
void Matrix_Sparse<E>::reSize(unsigned newsize) {
    if ( newsize < m_size ) { // hay más valores en uso que "newsize"
        // assert((m_size < newsize) && "=>Matrix_Sparse<E>::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<m_size; ++k) {
            new_I[k]   = m_I[k];
            new_J[k]   = m_J[k];
            new_val[k] = m_val[k];
        }
        delete [] m_I;  // Elimina los vectores viejos
        delete [] m_J;
        delete [] m_val;
    }
    m_I = new_I;      // Instala los vectores nuevos
    m_J = new_J;
    m_val = new_val;
    m_capacity = newsize;
}

/// \copydoc Matrix_BASE::add_Matrix()
template <class T>
void add_Matrix( Matrix_Sparse<T>& Res, const Matrix_Sparse<T> & M ) {
    // verifica que las dos matrices sean del mismo tamaño
    assert( "Matrix_Sparse<E>::add()" && ( Res.rows() == M.rows() ) );
    assert( "Matrix_Sparse<E>::add()" && ( Res.cols() == M.cols() ) );

    // Como las matrices son del mismo tamaño basta sumarlas entrada por entrada.

    typename Matrix_BASE<T>::value_type *pRes =   Res.m_val;
    typename Matrix_BASE<T>::value_type *pM   = & M.m_val[0];
    typename Matrix_BASE<T>::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
    - <code> rows() == O.rows() && cols() == O.cols() </code>

    \par Complejidad:
         O( <code> rows() * cols() </code> )

    \remarks
    - Esta es la implementación de
      Matrix_Sparse<E> operator-( Matrix_Sparse<E>&, Matrix_Sparse<E> )
    - El compilador tiene problemas en compilar un función amiga (<em>"friend"</em>)
      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 <class T>
void substract_Matrix( Matrix_Sparse<T>& Res, const Matrix_Sparse<T> & M ) {
    // verifica que las dos matrices sean del mismo tamaño
    assert( "Matrix_Sparse<E>::substract()" && ( Res.rows() == M.rows() ) );
    assert( "Matrix_Sparse<E>::substract()" && ( Res.cols() == M.cols() ) );

    // Como las matrices son del mismo tamaño basta restarlas entrada por entrada.

    typename Matrix_BASE<T>::value_type *pRes =   Res.m_val;
    typename Matrix_BASE<T>::value_type *pM   = & M.m_val[0];
    typename Matrix_BASE<T>::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