Universidad de Costa Rica Escuela de Ciencias de la Computación e Informática

Una Clase Matriz Chirrisquitica^[*] Escrita en C++

Abstract   

Esta es implementación C++ de matrices con plantillas es muy simple.

This is a simple implementation of the Matrix C++ class with templates.

      Con alguna frecuencia mis alumnos me preguntan si puedo facilitarles una matriz funcionalmente completa que puedan usar para implementar sus algoritmos en el curso de CI-1303 de Estructuras de Datos y Análisis de Algoritmos. Aunque existen muchas implementaciones disponibles en Internet, es interesante para ellos contar con una clase simple que puedan usar como base para realizar sus trabajos. Después de todo, en eso consiste la reutilización de componentes de programación.

      Las matrices son contenedores que almacenan valores numéricos, para los que las operaciones aritméticas básicas están definidas. Aunque lo usual es que los valores almacenados en la matriz sean números de punto flotante, de tipo "float" o "double", también es relevante usar otros tipos de datos, como números complejos. Por eso, tiene sentido implementar las matrices con plantillas.

      Esta implementación no es particularmente eficiente, pues existen muchas bibliotecas que lo son, como por ejemplo GMM++ [RP-2004] o GSL [GDTG-2004], o cualquiera de las mencionadas en [VELD-2004]. Talvez la particularidad de estas matrices, además de que cuentan con las operaciones aritméticas, es que están almacenados usando un vector lineal, en lugar de usar un vector de vectores como ocurre en otras implementaciones. Los algoritmos para matrices se pueden programar de manera muy natural con ésta implementación.

Especificación   

      Los objetos contenedores requieren de todas las operaciones básicas que es posible definir para cualquier variable, por lo que en su especificación hay que constructores y destructores, operadores de copia y comparación, etc.

      Para usar matrices, las operaciones más importantes son las aritméticas, que permitan sumar, restar y multiplicar matrices. Otras operaciones muy importantes, que no están incluidas en esta implementación, son las operaciones de inversión que permite resolver sistemas de ecuaciones representados como matrices cuadradas. Tampoco aquí está implementado el algoritmo para encontrar el determinante de la matriz, pero sí se incluyen algunas funciones como "isDiagonal()" o "isNull()" que sirven para determinar si la matriz tiene una propiedad o atributo particular.

      En algunas bibliotecas de matrices el acceso al valor "M(i,j)" de una matriz se hace usando la notación C++ de acceso a matrices, por lo que el programador cliente en su programa escribe "M[i][j]". Para lograr ese efecto es necesario incluir una clase intermedia dentro de la clase "Matrix" en donde hay que sobrecargar 2 veces "operador[]". Aquí fue preferido el camino simple de usar "at()" o la sobrecarga directa de "operador(unsigned,unsigned)" en sus 2 versiones: const y modificadora. Por ejemplo, para la expresión:
      M(i,j) = N(j,i);
el valor de "M" debe usar la versión de "Matrix::operator(unsigned,unsigned)" que permite modificar el valor almacenado en "M", mientras que para acceder al valor de "N(j,i)M no hace falta modificar "N". La lista completa de campos y métodos, tanto públicos como privados, está en la documentación generada por Doxygen:
      Documentación → http://www.di-mare.com/adolfo/p/Matrix/classMx_1_1Matrix-members.html

Modelo de la clase   

      El modelo de una estructura de datos es un diagrama en el que se muestra la relación entre los campos que incluye la clase. La mente del programador generalmente trabaja del detalle concreto a la abstracción, por lo que es usual que comprenda mejor qué es una clase si puede examinar su diagrama.

/ \ | 0 1 2 3 | m_rows == 2 | 4 5 6 7 | m_cols == 4 \ / +---+ | 2 | M(i,j) ==> m_val[ (i * m_cols) + j ] +---+ | 4 | +---+ +---+---+---+---+---+---+---+---+ | *-|-->| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | +---+ +---+---+---+---+---+---+---+---+

Figura 1

      Para representar una matriz se necesita conocer sus dimensiones en filas y columnas que hay que almacenar en 2 campos, que en esta implementación se llaman "m_rows" y "m_cols". El campo "m_val" es un puntero hacia la memoria dinámica en donde están almacenados en forma consecutiva los valores de la matriz. Debido a que en C++ los índices de los vectores comienzan en cero "unsigned(0)" en una matriz "2x4" como la matriz "M[]" de la Figura 1 existen las entradas "M(0,0)" y "M(1,3)", pero los índices de "M(2,1)" y "M(1,4)" están fuera de rango. Si los valores de la matriz están almacenados por filas en la memoria de la computadora, el valor "M(i,j)" estaría almacenado "m_cols" posiciones hacia la derecha de donde está "M(i-1,j)", por lo que para llegar a la posición "j" hay que saltarse "m_cols" valores.

/ \ | a e | | b f | m_rows == 4 | c g | m_cols == 2 | d h | \ / +---+ | 4 | M(i,j) ==> m_val[ i + (j * m_rows) ] +---+ | 2 | +---+ +---+---+---+---+---+---+---+---+ | *-|-->| a | b | c | d | e | f | g | h | +---+ +---+---+---+---+---+---+---+---+

Figura 2

      En la Figura 2 está otra matriz, pero en este caso sus valores están almacenados por columnas, por lo que la fórmula de acceso a la entrada "M(i,j)" sirva para saltarse "m_rows" filas para pasar de una columna a la siguiente. En C++ los valores de las matrices se almacenan por filas, pero en otros lenguajes como Fortran se almacenan por columnas.

Modelo para la matriz rala   

      Una matriz rala es aquella en que la mayor parte de los valores de la matriz son iguales. Usualmente, en una matriz rala el valor que se repite es el cero. Al trabajar con matrices ralas surge la natural tendencia de los programadores de evitar duplicar el valor que se repite, almancenado únicamente aquellos valores que son diferentes.

____________________________________ / 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

Figura 3

      Como se puede notar al comparar el modelo de la matriz de la Figura 1 con el de rala matriz matriz que está en la Figura 3, se necesitan más campos en el Rep para representar una matriz rala. En esta implementación se han usado 3 vectores paralelos en los que se almacenan las coordenadas "(i,j)" de valor almacenado en "M(i,j)" que está en el vector "M.m_val[]". La implementación de la versión modificadora de "Sparse_Matrix::operator(unsigned,unsigned)" es especial, pues al acceder un valor "M(i,j)" que no estaba ya incluido en "M.m_val[]" es necesario agregarlo, para que el programador cliente de la clase pueda modificarlo. Esto puede resultar en que no todos los valores del vector "m_val[]" sean diferentes del valor "m_same", que es el que se repite en toda la matriz.

      Las 2 implementaciones de la matriz son funcionalmente equivalentes, aunque la matriz rala tiene algunas operaciones adicionales para permitirle al programador definir el valor que se repite más en la matriz.

Otras consideraciones de la implementación   

      Usar plantillas para programar no es sencillo. De hecho, Stroustrup [Str-98] recomienda escribir los algoritmos primero sin usar plantillas, para depurarlos, y luego resolver las dificultades inherentes a la parametrización, usando estrategias similares a las expuestas en [DiM-2004]. Para esta implementación se usó el compilador VC++ v7 de Microsoft, que tiene algunos compartimientos extraños que es difícil explicar (también funciona con el C++ de GNU).

      Al implementar los operadores aritméticos hubo que enfrentar varios problemas que no reporta el compilador, sino más bien el ligador de eslabonamiento ("linker"). Por ejemplo, para el operador de suma de matrices:
      Matrix<E> operator + (const Matrix<E>& A, const Matrix<E>& B);
ocurre que el ligador de eslabonamiento no puede encontrar la implementación, pese a que fue declarada como función amiga ("friend") dentro de la clase. La forma de sortear este obstáculo es usar el método void Matrix<E>::add() como recurso para implementar como función desarrollada en línea ("inline") dentro de la clase. Otro obstáculo similar aparece en la implementación del operador de indexación "M(i,j)"
      const Matrix<E>::value_type& Matrix<E>::operator() ( unsigned , unsigned ) const;
pues el compilador no es capaz de identificar el tipo redefinido "Matrix<E>::value_type&" y obliga a que en la implementación haya que usar directamente el identificador "E" que es el parámetro de la plantilla:
      template <class E> const E& Matrix<E>::operator() ( unsigned , unsigned ) const;
Definitivamente, estas anomalías dificultan el uso de plantillas y dificultan entender el código.

A--> / \ B--> | 0 1 2 3 | m_rows == 3 C--> | 4 5 6 7 | m_cols == 4 D--> | 8 9 A B | m_refs == 5 E--> \ / +---+ +---+ A | *-|---+---|-* | B +---+ | +---+ C v E +---+ +---+ +---+ | *-|--------> |<5>| m_refs <--|-* | +---+ | +---+ +---+ | | 2 | m_rows D | +---+ +---+ | | 4 | m_cols | *-|---+ +---+ +---+---+---+---+---+---+---+ +---+ m_val | *-|-->| 0 | 1 | 2 | 3 | 4 | 5 | 6 | +---+ +---+---+---+---+---+---+---+

Figura 4

      Otra consideración importante que es consecuencia de esta implementación es que la evaluación de expresiones matriciales es bastante ineficiente. Por ejemplo, al evaluar la expresión:
      A = B + C * (D - E) * F;
el resultado obtenido es el correcto, pero a costa de copiar varias veces resultados intermedios que surgen de la necesidad de copiar y destruir las matrices temporales que el compilador necesita crear para evaluar la expresión. Una forma de mejorar esta implementación es permitir que varias matrices compartan el mismo valor, conservando una copia de referencia para saber cuando destruir o no el valor almacenado. Esta estrategia acelera mucho la velocidad de la operación Matrix<E>::copy() que se puede implementar copiando un puntero e incrementando la cuenta de referencia que corresponde al valor almacenado de la matriz. En la Figura 4 se muestra cómo las 5 matrices "A", "B", "C", "D" y "E" comparten el mismo valor almacenado, y en el campo "m_refs" se indica que hay 5 matrices que comparten el valor. Si en algún momento hay que cambiar el valor de alguna de las matrices, por ejemplo si se ejecuta la operación
      A(0,0) = Matrix<E>::value_type(0);
es necesario decrementar a "4" el campo "m_refs", luego hacer una copia completa de la matriz de manera que en "A(0,0)" quede el valor cambiado. En realidad, una implementación que usara cuentas de referencia lo que logra es posponer la copia del valor almacenado de la matriz, pero es un mecanismo efectivo para evitar que operaciones simples, como por ejemplo:
      A = B + C;
resulten en copia y destrucción inmediata de resultados intermedios, pues aunque la operación de suma retorna una nueva matriz temporal, al copiarla para sustituir el valor de "A" lo que se hace es compartir el valora almacenado entre "A" y la matriz temporal; como inmediatamente hay que destruir la matriz temporal, lo que cabe es decrementa a "1" la cuenta de referencia, por lo que efectivamente su uso evita crear la matriz temporal para copiarla, y luego destruirla. Este pequeño truco disminuye en 2 órdenes de magnitud el costo de copiar una matriz, pues con cuentas de referencias la complejidad es O("1") mientras que sin ellas es O( rows() * cols() ). Para un estudiante esta observación es impresionante.

      En la implementación de la suma y la resta no hace falta calcular la posición de la entrada "M(i,j)", pues como las matrices que hay que sumar tienen sus valores en las mismas posiciones dentro del vector al que apunta "m_val", basta recorrer los vectores para obtener el resultado correcto. Este truco también sirve para implementar la operación de comparación.

      El valor "Matrix<E>::value_type(0)" se usa para obtener el neutro del valor contenido en la matriz, y "Matrix<E>::value_type(1)" el neutro para la multiplicación.

      En muchas ocasiones no conviene que un error en el uso de los índices resulte en la terminación del programa. Por ejemplo, a los estudiantes les ocurre muchas veces que algunos algoritmos les funcionan casi siempre, pero les dan problemas con casos extremos, en los que cancelar el programa más bien dificulta la depuración porque no se pueden obtener todos los resultados esperados, sino que hay que revisarlos con cuenta gotas, uno por uno. Para estos casos resulta mejor usar el mecanismo de excepciones del lenguaje para atajar las fallas que se produzcan. En la implementación de Matrix.h se ha incluido la macro ASSERT_NO_ABORT que permite usar compilación condicional para modificar el comportamiento de la macro assert() de manera que use excepciones en lugar de cancelar el programa. Para usar excepciones, antes de incluir el archivos de encabezado Matrix.h hay que definir la macro ASSERT_NO_ABORT de esta manera:

#define ASSERT_NO_ABORT
#include "Matrix.h"
{
    Matrix M; /* ... */
    try { Inversion(M); }
    catch ( std::out_of_range ex ) {
        std::cout << ex.what() << " en algoritmo Inversion()" ;
    }
}

      Si alguna de las condiciones assert() de Matrix.h falla, en lugar de cancelarse el programa, lo que ocurre es una excepción de tipo std::out_of_range será levantada, la que puede luego ser atajada en el bloque catch(){} respectivo.

Conclusión   

      Un estudiante que examine esta implementación tiene la oportunidad de ver no sólo cómo usar plantillas, sino también cómo especificar sus programas. El código fuente está aquí:

Archivo	Descripción	Texto	Doxygen
Matrix.h	Declaraciones y definiciones para la clase Matrix http://www.di-mare.com/adolfo/p/Matrix/Matrix.h	.txt	.html
Sparse_Matrix.h	Declaraciones y definiciones para la clase Sparse_Matrix http://www.di-mare.com/adolfo/p/Matrix/Sparse_Matrix.h	.txt	.html
Matrix_Lib.h	Algunas propiedades para la clase Matrix http://www.di-mare.com/adolfo/p/Matrix/Matrix_Lib.h	.txt	.html
test_Matrix.cpp	Ejemplos de uso de la clase Matrix http://www.di-mare.com/adolfo/p/Matrix/test_Matrix.cpp	.txt	.html
test_Matrix.vcproj	Compilación VC++ v7 http://www.di-mare.com/adolfo/p/Matrix/test_Matrix.vcproj	.txt	.vcproj
test_Matrix.dev	Compilación GNU C++ http://www.di-mare.com/adolfo/p/Matrix/test_Matrix.dev	.txt	.dev
Matrix.dxg	Configuración Doxygen [v 1.3.9.1] ftp://ftp.stack.nl/pub/users/dimitri/doxygen-1.3.9.1-setup.exe	.txt	.dxg

Documentación Doxygen → http://www.di-mare.com/adolfo/p/Matrix/index.html

Agradecimientos   

      Alejandro Di Mare aportó varias observaciones y sugerencias importantes para mejorar este trabajo. La documentación completa para los programas fue obtenida con Doxygen [VANH-2004] y está disponible en este sitio Internet:
      http://www.di-mare.com/adolfo/p/Matrix/index.html

Notas de pie de página   

Bibliografía   

Indice   

[-]	Abstract
[1]	Especificación
[2]	Modelo de la clase
[3]	Modelo para la matriz rala
[4]	Otras consideraciones de la implementación
[-]	Conclusión
[-]	Agradecimientos
[-]	Código fuente
	Notas de pie de página
	Bibliografía
	Indice
	Acerca del autor
	Acerca de este documento
	Principio Indice Final

Acerca del autor   

Adolfo Di Mare <adolfo@di-mare.com>

Acerca de este documento   

Referencia:	Di Mare, Adolfo: Una Clase Matriz Chirrisquitica Escrita en C++ , Reporte Técnico ECCI-2004-02, Escuela de Ciencias de la Computación e Informática, Universidad de Costa Rica, 2004.
Internet:	http://www.di-mare.com/adolfo/p/Matrix.htm
Autor:	Adolfo Di Mare <adolfo@di-mare.com>
Contacto:	Apdo 4249-1000, San José Costa Rica Tel: (506) 207-4020       Fax: (506) 438-0139
Revisión:	ECCI-UCR, Noviembre 2007
Visitantes: