Paralelizaci´on en sistemas multicore con Maple · acceder desde programas escritos en lenguaje C,...

Paralelizacion

en sistemas multicorecon Maple

Alfonso Lopez Murcia

Indice

1. Introduccion 11.1. El punto de partida: master IMACI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2. Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.3. Estado del arte en computacion simbolica en paralelo . . . . . . . . . . . . . . . . . . . . . . . . . 11.4. Estructura de la memoria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2. Herramientas Maple para calculo paralelo 42.1. C OpenMaple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2. Maple Grid Computing Toolbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3. Algoritmo de Berkowitz y Subresultantes 103.1. Algoritmo Mejorado de Berkowitz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3.1.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113.2. Subresultantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.3. Otras maneras de definir Subresultantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

4. Implementacion del Algoritmo de Berkowitz en Maple Grid Computing Toolbox 164.1. Versiones del codigo Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164.2. Resultados experimentales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174.3. Analisis de los tiempos obtenidos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5. Trabajo futuro 21

A. Anexo A: Mathematica vs Maple 23A.1. Motivacion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23A.2. Calculo de determinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23A.3. La resultante de dos polinomios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24A.4. Analisis y conclusiones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

B. Anexo B: OpenMaple API 27B.1. Funciones especıficas de OpenMaple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27B.2. Funciones de conversion de Maple a estructuras de datos en C . . . . . . . . . . . . . . . . . . . . 27B.3. Funciones de conversion de C a objetos Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28B.4. Funciones para mostrar objetos Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28B.5. Funciones de comprobacion de tipo de dato Maple . . . . . . . . . . . . . . . . . . . . . . . . . . 28B.6. Funciones del sistema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29B.7. Funciones de manipulacion de tablas y listas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29B.8. Funciones de asignacion y seleccion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29B.9. Funciones de gestion de memoria y objetos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30B.10.Funciones de evaluacion y manejo de errores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31B.11.Funciones de manipulacion de tablas rectangulares (rtables) . . . . . . . . . . . . . . . . . . . . . 31

C. Anexo C: Algoritmo de Berkowitz en OpenMaple 33C.1. Codigo C del Algoritmo de Berkowitz en secuencial . . . . . . . . . . . . . . . . . . . . . . . . . . 33C.2. Compilacion y ejecucion del codigo C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

D. Anexo D: Implementacion OpenMaple de Berkowitz vs BerkowitzAlgorithm de Maple 41D.1. Tiempo para BerkowitzAlgorithm de Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41D.2. Tiempo para la implementacion OpenMaple de Berkowitz . . . . . . . . . . . . . . . . . . . . . . 41D.3. Matriz polinomial empleada en la toma de tiempos . . . . . . . . . . . . . . . . . . . . . . . . . . 41

i

E. Anexo E: Grid Computing Toolbox for Maple 43E.1. Primer programa paralelo. Reparto de matrices Toeplitz entre los procesadores . . . . . . . . . . 43E.2. Segundo programa paralelo. Reparto de matrices y producto casi paralelo . . . . . . . . . . . . . 47E.3. Tercer programa paralelo. Reparto de matrices Toeplitz y mejora del producto de matrices . . . 50E.4. Cuarto programa paralelo. Anillo de nodos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

Bibliografıa 62

ii

Agradecimientos.A Gema, Domingo, Vicente y mi familia.

Gracias por vuestra inestimable ayuda.

1. Introduccion

1.1. El punto de partida: master IMACI

El desarrollo de los ordenadores ha propiciado en los ultimos cuarenta anos el estudio de diversos problemasmatematicos ligados al calculo con cantidades masivas de objetos geometricos (puntos, rectas, triangulos, etc.) oalgebraicos (matrices, polinomios, etc.) relativamente simples. El modo de abordar los mismos viene determinadopor tres condiciones: el planteamiento constructivo (algorıtmico) de las soluciones, las implementaciones y elanalisis de su eficiencia. A su vez, el desarrollo de tecnicas especıficas para este calculo algebraico y geometricoha dado lugar a la formulacion de nuevos e interesantes problemas matematicos. Las aplicaciones tecnologicasson multiples.

Hay una gran variedad de herramientas para calculo algebraico. Algunas de ellas, como Maple, Matlab,Mathematica y librerıas matriciales se han visto en la asignatura Herramientas Informaticas para ComputacionCientıfica (HICC) dentro del master de Ingenierıa y Matematica Aplicadas en Ciencia e Ingenierıa (IMACI). Lasposibilidades de estas herramientas y la disponibilidad en la actualidad de procesadores con varios nucleos hanmotivado la determinacion de implementar un algoritmo paralelo, en concreto el Algoritmo de Berkowitz. Portanto tenemos el planteamiento constructivo inicial, seguidamente se implementara el Algoritmo de Berkowitzen paralelo y se evaluara su rendimiento.

Como necesitamos utilizar software para calculo simbolico, tras un analisis (ver anexo A) se decidio usarMaple junto con otro software para el calculo en paralelo en un procesador con cuatro nucleos.

1.2. Maple

Maple es un software matematico comercial de proposito general capaz de realizar calculos simbolicos,algebraicos y de algebra computacional. De las diversas funciones disponibles en Maple hay dos relacionadascon el calculo del polinomio caracterıstico y una con el algoritmo de Berkowitz:

charpoly dentro del paquete linalg.

CharacteristicPolynomial dentro del paquete LinearAlgebra.

BerkowitzAlgorihm dentro del paquete LinearAlgebra[Generic].

El algoritmo de Berkowitz se distingue porque calcula el polinomio caracterıstico sin realizar ninguna division.Este aspecto se trata en detalle en la seccion 3 (Algoritmo de Berkowitz y Subresultantes).

Sin embargo Maple 11 no utiliza todos los nucleos de un procesador multicore, pese a afirmar lo contrariola herramienta. Por este motivo se opto por implementar el algoritmo paralelo de Berkowitz y contrastar lostiempos de las funciones anteriores con la implementacion motivo de esta tesis de master. Como Maple nopermite el paralelismo es necesario un software adicional.

1.3. Estado del arte en computacion simbolica en paralelo

Por computacion simbolica entendemos la transformacion de expresiones matematicas a su forma simbolica.Los sistemas y programas que proporcionan un lenguaje para realizar computacion simbolica se denominan

Sistemas de Algebra Computacional (CAS, del ingles Computer Algebra System). Ademas existen librerıas paralos lenguajes de programacion actuales.

Actualmente los CAS se dividen en dos categorıas:

1. Sistemas de proposito general, como:

Axiom (libre). Web http://www.axiom-developer.org/

Macsyma (comercial). Web http://www.symbolics.com/

GNU Maxima (libre). Web http://maxima.sourceforge.net/

MuPAD (comercial). Web http://www.sciface.com/

Maple (comercial). Web http://www.maplesoft.com/

1

Mathematica (comercial). Web http://www.wolfram.com/

2. Sistemas especializados, algunos ejemplos son:

CoCoA (para calculos polinomicos con enfasis especial en las necesidades del algebra conmutativa ybases de Grobner). Web http://cocoa.dima.unige.it/

MAGMA (para calculos algebracos, teorıa de numeros, geometrıa y combinatoria). Webhttp://magma.maths.usyd.edu.au/

SINGULAR (para calculos polinomicos con enfasis especial en las necesidades del algebra conmutati-va, de la geometrıa algebraica y de la teorıa de singularidades). Web http://www.singular.uni-kl.de/

JACAL (para la simplificacion y manipulacion de ecuaciones y expresiones algebraicas con una ovarias variables). Web http://swissnet.ai.mit.edu/∼jaffer/JACAL.html

PARI-GP (computacion de curvas). Web http://pari.math.u-bordeaux.fr/

GAP (teorıa de grupos). Web http://www.gap-system.org/

La mayorıa de los CAS no estan preparados para el calculo en paralelo. Cuando se utiliza un CAS en unentorno paralelo normalmente la infraestructura subyacente es un grid, pues permite tratar todos los recur-sos como un unico superordenador de forma transparente. El software que facilita la conectividad entre lasplataformas heterogeneas del grid ofreciendo un conjunto de servicios que hacen posible el funcionamiento deaplicaciones distribuidas se denomina middleware. Un ejemplo de middleware para CAS es NetSolve/GridSolve.Hay otros comerciales como MathGridLink para Mathematica. En [12] puede verse un analisis sobre el grid.

Como no vamos a utilizar un cluster de ordenadores sino un unico equipo con cuatro nucleos veamos lasalternativas posibles que permiten extender Maple a un entorno paralelo:

Distributed Maple.

PVmaple.

El paquete Threads de Maple.

OpenMaple y MPI.

Grid Computing Toolbox.

Distributed Maple es una extension de Maple para para distribuir computaciones entre un grupo de esta-ciones de trabajo agrupadas como una maquina paralela virtual. Distributed Maple permite a las aplicacionesMaple interoperar entre ordenadores heterogeneos mediante el lenguaje de programacion Java. El softwareesta disponible en http://www.risc.uni-linz.ac.at/software/distmaple. Este software no esta actualizado y enteorıa funciona con versiones antiguas de Maple. Se ha intentado configurar Distributed Maple con Maple 11pero desafortunadamente no ha sido posible la integracion. Por el motivo anterior se descarto Distributed Maple.

Parallel Virtual Maple (PVMaple) utiliza el entorno de paso de mensajes PVM para la intercomunicacionde aplicaciones Maple entre equipos heterogeneos. Los principios de diseno son similares a Distributed Maple.El software esta disponible en http://www.info.uvt.ro/∼petcu/pvmaple. La ultima version disponible es paraMaple 6, por tanto no es viable con Maple 11.

En [20] puede verse un estudio de Distributed Maple y PVMaple.El paquete Threads de Maple 11 permite realizar un mayor numero de computaciones en Maple, incremen-

tando el rendimiento en sistemas multicore, como es nuestro caso. Sin embargo en Maple 11 no funciona elpaquete Threads. Por este motivo no se ha avanzado en esta alternativa. Recientemente Maplesoft anuncio laversion 12 de Maple donde asegura que funciona esta caracterıstica. Queda para un futuro una evaluacion deeste paquete en entornos multicore.

OpenMaple es un Application Program Interface (API) para efectuar llamadas a funciones y datos Mapledesde programas escritos en lenguaje C, Java o Visual. Tambien permite enriquecer el catalogo de funcionesMaple con funciones escritas en otro lenguaje. Su posible uso con MPI se detalla en la proxima seccion (ver2.1).

Grid Computing Toolbox es un software comercial para dotar de paralelismo a Maple. En la seccion 2.2 sepuede ver en detalle.

2

1.4. Estructura de la memoria

En la seccion 2 se describen las herramientas de Maple para el calculo en paralelo. La seccion 3 introduceel concepto de polinomio Subresultante y el Algoritmo de Berkowitz. En la seccion 4, se comentan los distintosalgoritmos en paralelo realizados en la memoria y conclusiones. Finalmente, se presentan lıneas de investigaciona desarrollar en la seccion 5.

Respecto los anexos, en el anexo A se compara Maple y Mathematica. El anexo B introduce el API deOpenMaple y los anexos C y D presentan el estudio de la implementacion del Algoritmo de Berkowitz conOpenMaple. Finalmente, el anexo E describe las diferentes implementaciones en paralelo del algoritmo con GridComputing Toolbox.

3

2. Herramientas Maple para calculo paralelo

2.1. C OpenMaple

El software matematico Maple dispone de un conjunto de funciones denominado OpenMaple que permiteacceder desde programas escritos en lenguaje C, Java o Visual Basic a los algoritmos y datos de Maple. En elcaso del lenguaje C el API disponible facilita iniciar una sesion Maple, ejecutar comandos Maple, manipularestructuras de datos Maple y controlar la salida. Para Java y Visual Basic, sus respectivas APIs facilitan lasmismas acciones. En la documentacion de Maple viene detallada toda la API de OpenMaple. Vease tambien elanexo B para un resumen de las funciones mas notables.

Message Passing Interface (MPI) es un estandar que define la sintaxis y la semantica de las funciones con-tenidas en una librerıa de paso de mensajes disenada para ser usada en programas que exploten la existenciade multiples procesadores. MPI no precisa de memoria compartida, por lo que es muy importante en la progra-macion para sistemas distribuidos. Los elementos principales en el paso de mensajes son el proceso que envıa,el proceso que recibe y el mensaje.

El binomio OpenMaple y MPI es una alternativa interesante para implementar el algoritmo de Berkowitz.El primer paso ha sido conocer OpenMaple. Veamos unos ejemplos sencillos para comprender la filosofıa de estaAPI. En primer lugar vamos a ver una expresion Maple sencilla

r := int(sin(x), x);eval(r, x = 0);

El codigo en lenguaje C puede verse en el cuadro 1. La funcion textCallBack es utilizada para mostrar unresultado de Maple. Forma parte de la estructura MCallBackVector creada en la lınea 12 que se pasa comoargumento en la funcion StartMaple de la lınea 15.

Como se aprecia en la lınea 15, cuando se crea una sesion OpenMaple se asocia a un descriptor, en este casodenominado kv. Este descriptor se emplea en todas las funciones OpenMaple posteriores del programa.

Las expresiones Maple se evaluan en las lıneas 22 y 24. Como estamos dentro de un programa en lenguajeC necesitamos definir x = 0 en Maple, por este motivo es necesaria la lınea 23.

Pasemos a un segundo ejemplo. Supongamos que nos interesa determinar el polinomio caracterıstico de unamatriz empleando las funciones de Maple.

with(LinearAlgebra):M := Matrix([[2, 1, 4], [3, 2, 1], [0, 0, 5]]);CharacteristicPolynomial(M,x);with(LinearAlgebra[Generic]):(Z[‘0‘], Z[‘1‘], Z[‘+‘], Z[‘-‘], Z[‘*‘], Z[‘=‘]) := (0, 1, ‘+‘, ‘-‘, ‘*‘, ‘=‘):C := BerkowitzAlgorithm[Z](M):< x3|x2|x|1 > . C;

La salida Maple en ambos casos es −5 + x3 − 9 ∗ x2 + 21 ∗ x. El codigo OpenMaple puede verse en el cuadro2. Con EvalMapleStatement es muy facil pasar expresiones a Maple.

Este segundo ejemplo puede inducir a pensar que es sencillo utilizar OpenMaple. Veamos como una expresiontan sencilla como

for i to r-1 do S[i] := A[i,r] end do;

se vuelve compleja en nuestro programa OpenMaple del cuadro 3.En OpenMaple para poder trabajar con matrices (Matrix), vectores (Vector) y listas (Array) se utilizan

tablas rectangulares (rtables). Si en rt tenemos una matriz, para poder acceder a una posicion empleamos eltipo de dato RTableData especificando los ındices.

Finalmente se implemento en OpenMaple el algoritmo de Berkowitz (ver anexo C.1). Sin embargo se detec-taron algunos inconvenientes como:

La conversion frecuente de tipos.

4

1 #include <stdio.h>

2 #include <stdlib.h>

3 #include ”maplec.h”4 static void M DECL textCallBack(void *data, int tag, char *output)5 {6 printf(” %s\n”,output);7 }8 /* Main */9 int main(int argc, char *argv[])10 char err[2048]; /* command input and error string buffers */11 MKernelVector kv; /* Maple kernel handle */12 MCallBackVectorDesc cb = { textCallBack, 0, 0, 0, 0, 0, 0, 0 };13 ALGEB r; /* Maple data-structure */14 /* Initialize Maple */15 if((kv=StartMaple(argc,argv,&cb,NULL,NULL,err)) == NULL)16 {17 printf(“Fatal error,%s\n”,err);18 return(1);19 }20 /* Example 1: compute an integral */21 printf(“ Evaluate an integral: ”);22 r = EvalMapleStatement(kv,”int(sin(x), x);”);22 /* Example 2: assign x a value and reevaluate the integral */23 MapleAssign(kv,ToMapleName(kv,”x”,TRUE),ToMapleInteger(kv,0));24 r = MapleEval(kv,r);25 MapleALGEB Printf(kv,”\n Evaluated at x=0, the integral is:%a”,r);26 StopMaple(kv);27 return(0);28 }

Cuadro 1: Fichero ejemplo om1.c

1 EvalMapleStatement(kv,”with(LinearAlgebra) : ”);2 EvalMapleStatement(kv,”M := Matrix([[2, 1, 4], [3, 2, 1], [0, 0, 5]]);”);3 printf(”\n CharacteristicPolynomial(M,x) = ”);4 r = EvalMapleStatement(kv,”CharacteristicPolynomial(M, x);”);5 MapleALGEB Printf(kv,”%a”,r);6 EvalMapleStatement(kv,”with(LinearAlgebra[Generic]) :”);7 EvalMapleStatement(kv,”(Z[‘0‘], Z[‘1‘], Z[‘ + ‘], Z[‘ − ‘], Z[‘ ∗ ‘], Z[‘ = ‘]) := (0, 1, ‘ + ‘, ‘ − ‘, ‘ ∗ ‘, ‘ = ‘) :”);8 EvalMapleStatement(kv,”C := BerkowitzAlgorithm[Z](M) :”);9 printf(”\n BerkowitzAlgorithm[Z](M) = ”);10 EvalMapleStatement(kv,”< x3|x2|x|1 > .C;”);

Cuadro 2: Ejemplo segundo de OpenMaple

La complejidad y extension del codigo en lenguaje C.

El rendimiento, medido en tiempo de ejecucion no mejora frente a Maple. Ver anexo D.

Es necesario proteger algunas variables (como la matriz polinomial) del recolector de basura automaticode Java para evitar el borrado de datos que causan un fallo en la ejecucion.

Si a la lista anterior anadimos la programacion en MPI el codigo C final es enorme y muy complejo. Portodo lo anterior se descarto el uso de OpenMaple con MPI para paralelizar el algoritmo de Berkowitz.

5

1 M INT i, r, index[2];2 ALGEB rt; /* Rtable like matrix */3 RTableData rtd; /* Rtable cell data */4 for(i=1; i < r; i++)5 {6 index1[0] = i; index1[1] = r;7 rtd = RTableSelect(kv,rt,index1); /* A[i,r] */8 MapleTableAssign(kv,S,ToMapleInteger(kv,i),rtd.dag); /* S[i]:=A[i,r] */9 }

Cuadro 3: Ejemplo tercero de OpenMaple

Sin embargo OpenMaple es una herramienta interesante. Pues podemos hacer nuevas funciones en lenguajeC, incluirla en una biblioteca y usarlas dentro de Maple mediante el paquete ExternalCalling.

2.2. Maple Grid Computing Toolbox

Grid Computing Toolbox es un software comercial para calculo paralelo distribuido en Maple. Esta bibliotecafacilita un conjunto de funciones al estilo MPI para distribuir computaciones en un grupo de maquinas o CPUsde una misma maquina. Cada maquina o CPU se denomina nodo en el grid.

Cada nodo en el grid ejecuta un proceso grid servidor que debe lanzarse explıcitamente en cada nodo. Cadaproceso grid servidor descubre los otros procesor servidores automaticamente formando de forma transparenteel grid. Este descubrimiento automatico se puede deshabilitar cuando se usa en un entorno controlado, comoPortable Batch System (PBS), un software de control de trabajos (job scheduling). Por tanto tenemos dos fasesbien definidas:

1. Crear el grid.

2. Ejecutar un programa Maple en el grid.

Crear el grid Para crear el grid hay tres alternativas:

1. En modo grafico, desde Maple 11 podemos utilizar el worksheet PersonalGridServer para lanzar variosnodos, ver su estado y parar el servidor.

2. En modo texto, ejecutando un fichero. En MS Windows es gridserver.bat y en GNU Linux es el scriptgridserver.sh.

6

3. Desde un fichero Maple mediante el comando StartServer. Tambien existe la orden StopServer paraparar el servidor grid en una maquina.

Veamos un ejemplo de como lanzar el servidor desde un fichero Maple.

> with(Grid):

> port:=2000: cpus:=4: broadcast:="127.0.0.255": bcport:=4400: logfile:="logs/grid.log":

> Server[StartServer](port, cpus, broadcast,bcport,logfile");

1196

> (TotalNodes, AvailableNodes) := Status();

4,3

El comando Status se utiliza para verificar el numero de nodos disponibles antes de enviar un trabajo algrid. Como se aprecia, devuelve dos valores. El primero es el numero de nodos que componen el grid, en estecaso cuatro. El segundo valor especifica la cantidad de nodos disponibles para computacion, en este caso sontres.

Ejecutar un programa Maple en el grid Una vez creado el grid es el momento de enviar trabajos. Siestamos comenzando con la herramienta lo mas didactico es practicar desde el entorno grafico los cuatro ejemplosincluidos en el worksheet Grid.

Sin embargo lo normal es lanzar los trabajos de una de estas dos formas:

1. En MS Windows con el fichero launcher.bat y en GNU Linux con el script launcher.sh seguido de dosargumentos. El primer argumento corresponde al fichero con los comandos Maple y el segundo argumentoes el fichero con los nodos participantes.

2. Desde un fichero Maple mediante el comando Launch.

Veamos un ejemplo de como enviar un trabajo desde un fichero Maple al grid creado anteriormente.

> code := "printf(\"Hello from node %a of %a \\n\", Grid:-Util:-MyNode(),

Grid:-Util:-NumNodes());":

> result := Grid[Launch](numnodes, code, printf, proc() false end);

Node 2: Hello from node 2 of 4



> Server[StopServer](port,cpus);

0

En este ejemplo vemos dos comandos del toolbox. El comando NumNodes devuelve el numero de nodosdisponibles para la computacion. El comando MyNode devuelve el numero de nodo dentro del grid. Si hay Nnodos en el grid este numero varıa de 0 a N-1. En el ejemplo anterior hay tres nodos, cada uno numerado 0, 1

7

y 2 correspondientes a tres nucleos de un procesador Intel Quad Core. Como se aprecia, no aparecen los textosordenados por numero de procesador. Si se efectuase otra ejecucion aparecerıa en otro orden.

Grid Computing Toolbox sigue la misma filosofıa de MPI en cuanto a un nodo maestro (nodo 0) y paso demensajes entre nodos. Veamos un ejemplo de envıo de mensajes analizando el fichero message.mpl:

1 me := Grid:-Util:-MyNode():

2 n := Grid:-Util:-NumNodes():

3 dest := me+1:

4 if dest=n then dest := 0 fi:

6 if me=0 then # master

7 msg := [seq(i^2, i=1..n)]: # create a message

8 Grid:-Util:-Send(dest, msg); # send it on to node number 1

9 m := Grid:-Util:-Receive(); # wait for a message to come back

10 else # other nodes

11 msg := Grid:-Util:-Receive(); # wait for a message to arrive

12 Grid:-Util:-Send(dest, msg); # send it on to the next node

13 fi;

El nodo maestro (nodo 0) crea el mensaje inicial en la lınea 7. Si hay tres nodos disponibles en el gridentonces msg es [1, 4, 9]. Este mensaje se envıa en la lınea 8 al nodo 1 y el nodo 0 se queda a la espera de recibirdatos. El resto de nodos (nodo 1 y nodo 2) reciben el mensaje y lo envıan al siguiente. Es decir, el nodo 1 envıael mensaje al nodo 2. El nodo 2, sin embargo cumple la condicion de la lınea 4 y por tanto envıa el mensajeal nodo 0, cerrando el cırculo. Como el nodo 0 estaba pendiente de recibir el mensaje, termina la ejecucion delcodigo.

En el ejemplo anterior no aparece el comando Launch. En consecuencia es de suponer que el grid estabacreado con el script gridserver.sh y se ha ejecutado con el script launcher.sh de este modo:

gridserver.sh -n 4 -b 4400 -a 127.0.0.255 -p 2000 -f logs/gridlog.txt startmultiple &

launcher.sh -h localhost -n 4 message.mpl

Launching on 4 nodes

[1, 4, 9, 16]

El comando Send solamente bloquea al nodo mientras dura la transmision, por contra, el comando Receive

bloquea el nodo hasta la recepcion del mensaje. Es muy importante establecer el orden de envıos y recepcionespara evitar bloqueos (deadlocks).

Otra forma de distribuir computacion en el grid sin usar los comandos Send y Receive es mediante Seq yMap. El comando Seq es similar al comando seq de Maple pero distribuyendo el calculo de la secuencia en elgrid. El resultado esta disponible unicamente en el nodo 0. Una limitacion de Seq es que no permite estructurasde mas de una dimension como matrices.

r := Grid:-Util:-Seq(i^16, i=1..100000);

if Grid:-Util:-MyNode() = 0 then

printf("%q", r)

fi;

El comando Map es similar al comando map de Maple pero distribuyendo la computacion en el grid. Elresultado esta disponible unicamente en el nodo 0. Una limitacion de Map es que no permite estructuras demas de una dimension como matrices.

r := Grid:-Util:-Map(proc(a, b) a^2+b end, [seq(i, i=1..10)], 10);

if Grid:-Util:-MyNode() = 0 then

printf("%a", r)

fi;

Los comandos Seq y Map dejan de ser utiles en calculos iterativos donde todos los nodos deben tener losdatos, es decir, como el resultado esta en el nodo 0 es preciso enviarlo al resto de nodos para proseguir con laiteracion siguiente. En este caso recurrimos a Send y Receive directamente descartando Seq y Map.

8

Finalmente, debe destacarse el uso del comando Time para determinar el tiempo empleando en la resolucionde un algoritmo. Mientras el comando time de Maple mide el uso de la CPU el comando Time mide el tiempoen milisegundos. Esta segunda medida es mejor para trabajos en paralelo.

9

3. Algoritmo de Berkowitz y Subresultantes

Esta seccion tiene como objetivo describir el algoritmo implementado en paralelo y una aplicacion del mismo.El algoritmo de Berkowitz calcula el polinomio caracterıstico de una matriz, usando tan solo operaciones basicas(suma, restas y multiplicaciones). Esto implica que si una matriz tiene entradas en un dominio, aplicandoeste algoritmo se evita las divisiones exactas en su cuerpo de fracciones a la hora de calcula su polinomiocaracterıstico.

Como aplicacion, se presenta el calculo de subresultantes. En [3] y [15] se propone el calculo de subresultantesde polinomios con coeficientes con parametros, utilizando el Algoritmo de Berkowitz. En esta memoria se hancalculado las subresultantes de los polinomios del anexo 4.2, aplicando varios metodos descritos en la seccion3.2, tanto en secuencial como en paralelo.

3.1. Algoritmo Mejorado de Berkowitz

Como antes se ha mencionado, el algoritmo de Berkowitz calcula el polinomio caracterıstico de una matriz.Ademas, usa solo operaciones basicas (suma, restas y multiplicaciones) en el dominio base y tiene una buenacomplejidad secuencial. En [1] se puede encontrar un estudio detallado tanto de su complejidad secuencial comoparalela.

A continuacion se describen los aspectos mas notables del algoritmo y su mejora (para mas detalle, vease[1], [2] y [7]).

Proposicion 3.1. (Formula de Samuelson)Sea D un dominio. Sean r, n dos numeros enteros verificando 1 ≤ r < n y Ar (resp. Ar+1 ) la submatriz

principal directora de orden r (resp. r + 1 ) de una matriz dada A ∈ Mn(D).Se considera la siguiente particion de Ar+1 :

Ar+1 =

(Ar Sr

Rr ar+1,r+1

)

donde Sr ∈ Dr×1 y Rr ∈ D

1×r . Sean

Pr(x) = det(xIr − Ar) y Pr+1(x) = det(xIr+1 − Ar+1)

los polinomios caracterısticos de las matrices Ar y Ar+1 respectivamente con

Pr(x) = p0xr + p1x

r−1 + . . . + pr.

Entonces:

Pr+1(x) = (x − ar+1,r+1)Pr(x)

−r∑

k=1

[(RrAk−1r Sr)p0 + (RrA

k−2r Sr)p1 + . . . + (RrArSr)pk−2 + (RrSr)pk−1]x

r−k

En [7], S. J. Berkowitz presenta la Formula de Samuelson mediante el producto de ciertas matrices triangu-lares inferiores Toeplitz (n + 1) × n. Dado un polinomio

P (x) =

d∑

k=0

akxn−k,

con coeficientes en D, y su vector de coeficientes

−→P =

a0

a1

...ad

,

10

Toep (P ) denota la matriz (d+1)× d triangular inferior Toeplitz asociada a la lista (ak) 0≤k≤d de coeficientesde P (x):

Toep (P ) =

a0 0 · · · 0a1 a0 · · · 0...

.... . .

...ad−1 ad−2 · · · a0

ad ad−1 · · · a1

.

Con esta notacion, la Formula de Samuelson se puede escribir en forma de igualdad matricial:

−−→Pr+1 = Toep (Qr+1) ×

−→Pr

donde Qr+1 es el polinomio

Qr+1 = xr+1 − ar+1,r+1xr −

r−1∑

i=0

RrAirSrx

r−1−i

La aplicacion de esta relacion sobre los vectores−→Pr , . . . ,

−→P2 para r ≥ 2 junto con la igualdad:

−→P1 =

(1

−a11

)

= Toep (Q1)

permite describir el vector de coeficientes del polinomio caracterıstico de A como:

−→Pn = Toep (Qn) × Toep (Qn−1) × · · · × Toep (Q1).

Ademas, en [7], S. J. Berkowitz describe la implementacion del algoritmo en paralelo. En [1] se presenta unaversion secuencia mejorada del Metodo de Berkowitz que se describe a continuacion.

Algoritmo Mejorado de Berkowitz (IBA)

Entrada: Una matriz cuadrada de orden n, A ∈ Mn(D).Salida: El polinomio caracterıstico Pn de A .

(IBA.1) InicioEl vector Vect toma el valor:

Vect :=

(1

−a11

)

.

(IBA.2) para r desde 1 a n − 1 ,(IBA.2.1) Calculo de las entradas {RrA

k−1r Sr}

rk=1 de Toep (Qr+1);

(IBA.2.2) Actualizar Vect en Vect := Toep (Qr+1) × Vect ;

(IBA.3) Solucion

Pn sera el unico polinomio tal que−→Pn = Vect.

3.1.1. Ejemplo

A continuacion se muestra un ejemplo aplicando el Algoritmo Mejorado de Berkowitz para el calculo delpolinomio caracterıstico de una matriz con entradas en Z[y].

Se considera la matriz A definida como

A =

y2 + 2 y −y2 + 2 y −y2 + 2 −y2 + y + 1

2 y2 − y −y2 + 2 y − 1 −y2 + y + 1 −y2

2 y2 − y − 1 −y2 − y − 1 2 y2 − y + 1 −y2 + 1

−1 + 2 y y2 + 2 y y2 y2 + y − 1

.

11

Tras haber hecho los calculos necesarios, los polinomios Q2(x, y), Q3(x, y), Q4(x, y) son:

Q2 = x2 +(y2 − 2y + 1

)x +

(2y2 − y

) (y2 − 2y

),

Q3 = x3 −(2y2 − y + 1

)x2 +

(y4 − 4y2 − y3 + 4y + 3

)x + 4y − y6 − 1 − 8y4 +

+13y2 + 7y5 − 11y3,

Q4 = x4 −(y2 + y − 1

)x3 +

(2y4 + 4y3 − 4y2 − y + 1

)x2 +

(3y6 − 7y5 + 3y4 − 5y2

−4y + 2)x − y(y7 + 9y6 − 43y5 − 6y4 + 88y3 + 4y2 − 23y − 2

).

Tras multiplicar las matrices Toeplitz definidas por los coeficientes de estos polinomios, se obtiene:

4∏

i=1

Toep(Qi) =

1

−3y2 − 4y + 1

6y4 + 11y3 − 8y2 + 3y + 3

−6y6 − 11y5 − 15y4 + 32y3 − 32y2 − 4y + 6

−3y7 + 50y6 − 4y5 − 49y4 + 40y3 − 3y2 − 3y + 4

,

y con ello los coeficientes del polinomio caracterıstico de la matriz A.

3.2. Subresultantes

En algebra computacional, la Teorıa de Subresultantes proporciona un metodo potente para resolver dediferente naturaleza, que implican manipulacion de polinomios. Por ejemplo, podemos calcular el maximo comundivisor de dos polinomios con una variable con coeficientes enteros evitando el crecimiento exponencial si usamosel Algoritmo de Euclides.

Si se extiende la Teorıa de Subresultantes al caso multivariante, dados dos polinomios con coeficientes enun dominio, el calculo de sus Polinomios Subresultantes nos permite hallar el maximo comun divisor de dichospolinomios sin aplicar el Algoritmo de Euclides y sin la necesidad de trabajar en el cuerpo de fracciones de D.

El concepto de determinante polinomial asociado a una matriz con entradas en D proporciona una de lasmaneras mas usuales de introducir la nocion de Polinomio Subresultante.

Definicion 3.2. Sea ∆ una matriz m × n con m ≤ n. El determinante polinomial asociado a ∆, detpol(∆),viene definido por:

detpol(∆) =

n−m∑

k=0

det(∆k)xn−m−k

donde ∆k es la submatriz cuadrada de ∆ formada con las primeras m − 1 columnas de ∆ y la (k + m)–esimacolumna de ∆.

Observese que si ∆ es cuadrada, entonces detpol(∆) det(∆).

Definicion 3.3. Sean A(x), B(x) dos polinomios en D[x] y n, m ∈ N con gr(A) = n y gr(B) = m:

A = a0xn + a1x

n−1 + · · · + an B = b0xm + b1x

m−1 + · · · + bm.

Dado i ∈ {0, . . . , ınf(n, m) − 1}, la matriz de Sylvester de ındice i asociada a (A, B) se define de la siguientemanera:

Sylvi(A, B)

n+m−i︷︸︸︷

a0 . . . an

. . .. . .

a0 . . . an

b0 . . . bm

. . .. . .

b0 . . . bm

m − i

n − i

12

Definicion 3.4. La Resultante de A(x) y B(x), Resultante(A, B), se define como el determinante de la matrizde Sylvester,

Resultante(A, B) = det

n+m︷︸︸︷

a0 . . . an

. . .. . .

a0 . . . an

b0 . . . bm

. . .. . .

b0 . . . bm

m

n

Definicion 3.5. En las mismas condiciones que en la Definicion 3.3:

El Polinomio Subresultante de ındice i asociado a (A, B) se define como el determinante polinomial de lamatriz de Sylvester de ındice i asociada a (A, B):

Sresi(A, B) = detpol(Sylvi(A, B)).

La Subresultante Principal de ındice i asociada a (A, B) es el coeficiente del termino de grado i delPolinomio Subresultante Sresi(A, B):

psci(A, B) = coefi(Sresi(A, B)),

es decir, el determinante de la submatriz de Sresi(A, B) formada por las primeras n + m − 2i columnas.

Definicion 3.6. La Subresultante Sresi(A, B) es defectiva si su grado es estrictamente menor que i,

gr(Sresi(A, B)) < i,

o, equivalentemente, si psci(A, B) = 0. En caso contrario Sresi(A, B) es regular.

Definicion 3.7.

Se denomina Sucesion de Subresultantes asociada a A(x) y B(x) (con gr(A) ≥ gr(B)), a la coleccion depolinomios:

{Sresi(A, B)}i∈{0,...,gr(B)−1}

Se denomina Sucesion de Subresultantes Principales asociada a A(x) y B(x) (con gr(A) ≥ gr(B)), a lacoleccion de coeficientes:

{psci(A, B)}i∈{0,...,gr(B)−1}

La siguiente proposicion relaciona la Sucesion de Subresultantes con el maximo comun divisor de A(x) yB(x) con coeficientes en D.

Proposicion 3.8. Siguiendo con la misma notacion que en la Definicion 3.7, se verifica:

gr(mcd(A, B)) = i ⇐⇒

psc0(A, B) = . . . = psci−1(A, B) = 0

psci(A, B) 6= 0(1)

Ademas si se verifica cualquiera de las dos afirmaciones entonces Sresi(A, B) es un maximo comun divisorde A(x) y B(x).

Demostracion. Vease por ejemplo [14].Respecto a su calculo, el clasico Teorema de las Subresultantes (vease [14] o [18]) proporciona un algoritmo

recursivo basado en el calculo de pseudorestos y divisiones exactas para la determinacion eficiente de la Sucesionde Subresultantes asociada a los polinomios A(x) y B(x): conocidas los Polinomios Subresultantes Sj+1(x) y

13

Sj(x) (donde Sj+1(x) ha de ser regular y Sj(x) no necesariamente), el teorema indica cual es la expresion delos Polinomios Subresultantes de ındices j − 1, . . . , gr(Sj) − 1, en funcion de Sj+1(x) y Sj(x).

En [6] los autores introducen un algoritmo que mejora el clasico Teorema de las Subresultantes en el casodefectivo. La mejora se consigue evitando el aumento de de los coeficientes intermedios que sucede en el algoritmode las Subresultantes en el caso defectivo.

Ambos algoritmos eluden el uso el uso de determinantes e implican divisiones exactas y restos. Sin embargo,si los coeficientes son a la vez polinomios (por ejemplo, cuando dependen de parametros), estos algoritmospueden resultar costosos y lentos (ver [3]).

3.3. Otras maneras de definir Subresultantes

Una forma novedosa de definir i–esimo polinomio Subresultante se expone en [16].

Sresi(A, B) = (−1)i(m−i+1)

∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣

a0 a1 a2 . . . . . . an

. . .. . .

. . .. . .

a0 a1 a2 . . . . . . an

1 −x

. . .. . .

1 −x

b0 b1 b2 . . . . . . . . . bm

. . .. . .

. . .. . .

b0 b1 b2 . . . . . . . . . bm

∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣

m − i

i

n − i

(2)

Tambien podemos pensar en calcular Subresultantes a partir de Matrices de Bezout.

Definicion 3.9. La matriz de Bezout asociada a dos polinomios, A(x) y B(x), es la matriz simetrica

Bez(A, B) =

c0,0 . . . c0,n−1

......

cn−1,0 . . . cn−1,n−1

donde las entradas ci,j vienen definidas por la expresion

A(x)B(y) − A(y)B(x)

x − y=

n−1∑

i,j=0

ci,jxiyj .

Definicion 3.10. Dados A(x) y B(x), deg(A) = n ≥ deg(B) = m, la matriz de Bezout hıbrida asociada a A(x)y B(x), Hbez(A, B), es una matriz cuadrada de orden n,

Hbez(A, B) =

h0,0 . . . h0,n−1

......

hn−1,0 . . . hn−1,n−1

,

cuyas entradas vienen definidas por:

para 1 ≤ i ≤ m, 1 ≤ j ≤ n, la entrada (i, j)–esima es el coeficiente de xn−j en el polinomio

Km−i+1 = (a0xm−i + . . . + am−i)(bm−i+1x

n−m+i−1 + . . . + bmxn−m)

−(am−i+1xn−m+i−1 + . . . + an)(b0x

m−i + . . . + bm−i);

para m+1 ≤ i ≤ n, 1 ≤ j ≤ n, la entrada (i, j)–esima es el coeficiente de xn−j en el polinomio xn−iB(x).

14

Entonces se verifican las identidades siguientes (ver [4] y [22] para mas detalles)

(−1)(n−i)(n−i−1)/2pn−m0 Sresi(A, B)(−1)i

∣∣∣∣∣∣∣∣∣∣∣∣∣∣

cn−1,n−1 cn−1,n−2 . . . . . . . . . cn−1,0

...... . . . . . . . . .

...ci,n−1 ci,n−2 . . . . . . . . . ci,0

1 −x

. . .. . .

1 −x

∣∣∣∣∣∣∣∣∣∣∣∣∣∣

n − i

i

, (3)

Sresi(A, B)(−1)i(n−i+1)

∣∣∣∣∣∣∣∣∣∣∣∣∣∣

hn−1,0 hn−1,1 . . . . . . . . . hn−1,n−1

...... . . . . . . . . .

...hi,0 hi,1 . . . . . . . . . hi,n−1

1 −x

. . .. . .

1 −x

∣∣∣∣∣∣∣∣∣∣∣∣∣∣

n − i

i

. (4)

Todo esto nos lleva a considerar tres matrices diferentes para calcular la secuencia de Subresultantes a travesdel calculo de determinantes, bien usando la Ecuacion (2), la Ecuacion (3) o la Ecuacion (4).

Como antes se ha mencionado, si los coeficientes de los polinomios considerados son a su vez polinomios, sedebe evitar un algoritmo que implique el calculo de divisiones exactas, por el coste que esto implica.

Debido a que el termino independiente del polinomio caracterıstico de una matriz es justamente el deter-minantes de esta (salvando el signo), el calculo de los determinantes (2), (3) y (4) aplicando el Algoritmo deBerkowitz evita precisamente el computo de divisiones exactas.

Por ello, el objetivo de la memoria consiste en:

1. Implementar el algoritmo de Berkowitz en paralelo con Maple.

2. Una vez implementado, aplicar Ecuaciones (2) y (4) para hallar la secuencia de Subresultantes de lospolinomios que figuran en el anexo 4.2,

a) tanto en secuencial,

con la funcion LinearAlgebra[Determinant], implementada en Maple,

y con la funcion LinearAlgebra[Generic][Determinant](method=BerkowitzAlgorithm), im-plementada en Maple;

b) como en paralelo, aplicando nuestra implementacion del algoritmo de Berkowitz.

15

4. Implementacion del Algoritmo de Berkowitz en Maple Grid Com-

puting Toolbox

4.1. Versiones del codigo Maple

En un principio se abordo la implementacion desde la filosofıa MPI. Es decir, un procesador maestro paraestablecer y repartir el trabajo junto con un grupo de procesadores esclavos encargados de realizar el trabajo.Los calculos se pasan entre procesadores mediante intercambio de mensajes.

Maestro(Nodo 0)

(Nodo 1) (Nodo 2) (Nodo 3)

Esclavo Esclavo Esclavo

Send/Receive

Figura 1: Comunicacion maestro y esclavos.

Un primer enfoque ha sido implementar una bolsa de tareas. La idea subyacente es un procesador maestrocon la lista de ındices de las matrices Toep(Qi) a calcular. Los procesadores esclavos solicitan un ındice, recibenpor ejemplo, el valor j, calculan la matriz Toep(Qj), la envıan al procesador maestro y vuelven a solicitar otroındice. Sin embargo, esta estrategia de asignar de forma dinamica el trabajo a los procesadores esclavos noha funcionado bien pues se producen aleatoriamente abrazos mortales (deadlocks) abortando la ejecucion delprograma.

El segundo enfoque es una bolsa de tareas estatica. Es decir, el procesador maestro determina el trabajode cada procesador esclavo y envıa una lista de ındices. Cada procesador esclavo recibe la lista y calcula lasmatrices Toeplitz correspondientes, enviando al procesador maestro dichas matrices. Finalmente el procesadormaestro hace el producto Toep(Qn) x Toep(Qn−1) x ... x Toep(Q1). Puede verse el codigo Maple paralelo en elanexo E. Posteriormente se ha tratado de mejorar el codigo paralelo en las secciones E.2 y E.3.

Finalmente se ha descartado el modelo maestro-esclavo y se ha tratado a todos los procesadores por igual.Esta nueva filosofıa nos lleva a un anillo de nodos (ver seccion E.4) que calculan y pasan a su vecino los datos.

ReceiveSend

Send

Receive Send

Receive

SendReceive

Nodo 3Nodo 0

Nodo 2Nodo 1

Figura 2: Anillo de procesadores.

En las diversas pruebas realizadas, el anillo de procesadores es inferior su tiempo de ejecucion con respectoal modelo maestro-esclavo. En consecuencia, ha sido la distribucion de procesadores elegida para implementarel algoritmo paralelo de Berkowitz.

16

4.2. Resultados experimentales

Tiempos para matriz Sylvester. En la tabla siguiente se muestran los tiempos totales obtenidos al calcularel determinante de las distintas matrices Sylvi(p,q).

Polinomios Determinant Determinant(Berko) ParaleloP1 (11,10,2,5,denso) 1h, 2m, 46s 4h, 23m, 23s 3h, 39m, 44sP2 (8,7,3,5,no denso) 59s 55m, 21s 1h, 27m, 17sP3 (9,8,3,5,denso) 6h, 31m, 6s Cancelado a las 24h 29h, 48m, 40sP4 (9,8,3,5,no denso) 2m, 15s 2h, 16m, 5s Error MapleP5 (8,7,4,5,no denso) 36m, 17s Error Maple Error MapleP6 (8,7,4,5,denso) 37m, 34s 5h, 32m, 57s Error MapleP7 (14,6,2,5,no denso) 8m, 39s 54m, 32s Error MapleP8 (14,6,2,5,denso) 30s 2h, 12m, 50s 44m, 46sP9 (8,5,3,5,denso) 14m, 20s 41m, 34s 2h, 11m, 58sP10 (8,5,3,7,no denso) 19s 7m, 57s 22mP11 (9,4,3,5,denso) 6m, 5s 14m, 14s 22m, 57sP12 (10,6,3,5,no denso) 3h, 22m, 44s 5h, 15m, 57s 5h, 25m, 6sP13 (8,5,4,5,no denso) 5m, 38s 1h, 8m, 32s 4h, 48m, 40sP14 (8,5,4,6,denso) 4h, 7m, 8s 7h, 0m, 7s Error Maple

En la primera columna tenemos las caracterısticas del par de polinomios empleados:

(grado pi, grado qi, numero de parametros, grado total maximo de los coeficientes, densidad)

En la segunda columna tenemos el tiempo secuencial resultante de aplicar la funcion Determinant de Maple11.

En la tercera columna tenemos el tiempo secuencial total obtenido con la funcion Determi-nant(method=BerkowitzAlgorithm).

En la cuarta columna tenemos el tiempo paralelo resultante de la ejecucion en un anillo con cuatro proce-sadores (procesador Intel Quad Core con cuatro nucleos). Los valores de la maquina virtual Java de Maple 11son los valores por defecto (no se han modificado).

Puede apreciarse como el codigo paralelo no mejora el calculo secuencial.

Tiempos para matriz Bezout Hıbrida. Veamos los tiempos obtenidos para la matriz Bezout Hıbrida delos polinomios.

Polinomios Determinant Determinant(Berko) ParaleloP1 (11,10,2,5,denso) 1h, 10m, 4s 1h, 46m, 31s 36m, 32sP2 (8,7,3,5,no denso) 1h, 12m, 18s 23m, 46s 17m, 7sP3 (9,8,3,5,denso) 18h, 16m, 50s 12h, 34m, 15s 6h, 56m, 51sP4 (9,8,3,5,no denso) 4h, 46m, 4s 1h, 16m, 58s Error MapleP5 (8,7,4,5,no denso) 13m, 54s Error Maple Error MapleP6 (8,7,4,5,denso) 3m, 12s 6h, 36m, 27s Error MapleP7 (14,6,2,5,no denso) 13m, 11s 31m, 39s 22m, 50sP8 (14,6,2,5,denso) 36m, 52s 4h, 50m, 3s 2h, 23m, 37sP9 (8,5,3,5,denso) 1h, 32m, 55s 58m, 26s 45m, 19sP10 (8,5,3,7,no denso) 12m, 34s 4m, 19s 4m, 34sP11 (9,4,3,5,denso) 1h, 10m, 54s 32m, 56s 30m, 56sP12 (10,6,3,5,no denso) 3h, 43m, 52s 1h, 27m, 34s 48m, 12sP13 (8,5,4,5,no denso) 2m, 36s 1h, 3m 1h, 12m, 51sP14 (8,5,4,6,denso) 43m, 46s 17h, 12m, 58s 14h, 15m, 30s

En este caso si se aprecia una mejora con el codigo paralelo, si se exceptuan las ejecuciones con error Maple.Como se han producido errores de memoria Java se ha procedido a un ajuste en los valores de la maquina

virtual Java. Estos valores se especifican en la creacion del Grid y han sido:

17

-Xms128M -Xmx800m -XX:PermSize=64M -XX:MaxPermSize=256M

Los nuevos tiempos paralelos para las matrices Bezout Hıbridas han sido:

Polinomios Determinant Determinant(Berko) ParaleloP1 (11,10,2,5,denso) 1h, 10m, 4s 1h, 46m, 31s 37m, 34sP2 (8,7,3,5,no denso) 1h, 12m, 18s 23m, 46s 17m, 40sP3 (9,8,3,5,denso) 18h, 16m, 50s 12h, 34m, 15s 7h, 9m, 50sP4 (9,8,3,5,no denso) 4h, 46m, 4s 1h, 16m, 58s 52m, 40sP5 (8,7,4,5,no denso) 13m, 54s Error Maple 34m, 47sP6 (8,7,4,5,denso) 3m, 12s 6h, 36m, 27s 34h, 6m, 4sP7 (14,6,2,5,no denso) 13m, 11s 31m, 39s 22’, 45mP8 (14,6,2,5,denso) 36m, 52s 4h, 50m, 3s 2h, 24m, 48sP9 (8,5,3,5,denso) 1h, 32m, 55s 58m, 26s 46m, 29sP10 (8,5,3,7,no denso) 12m, 34s 4m, 19s 4m, 24sP11 (9,4,3,5,denso) 1h, 10m, 54s 32m, 56s 31m, 42sP12 (10,6,3,5,no denso) 3h, 43m, 52s 1h, 27m, 34s 53m, 12sP13 (8,5,4,5,no denso) 2m, 36s 1h, 3m 1h, 13m, 42sP14 (8,5,4,6,denso) 43m, 46s 17h, 12m, 58s 14h, 32m, 26s

Con el cambio en los valores de memoria de la maquina virtual Java se evitan los errores de Maple. Sinembargo se aprecia un ligero aumento en el tiempo paralelo. Por tanto es fundamental un ajuste de los valoresJava para una ejecucion optima de codigo paralelo.

Polinomios empleados. A continuacion se detallan los polinomios empleados en este estudio.

P1 (11, 10, 2, 5) :p1 = x11 + (2ab3 − 3a5 − a2b3)x10 + (8a4b − 5a2b3)x9 + (a3b2 + 3a2b3)x8 + (4a4b + 3a3b2)x7 +

(8 − 6a + 7b + 8a2 − 5a3b2)x6 + (−9 − 7a − 8ab − 8a2b + 9b5)x5 + (7a2 + 8a3 + 4a2b + 8a4 +8a3b + 9b5)x4 + (3a + 3a2 − 5b3 − 9a3b − 5ab3 − 4ab4)x3 + (5ab2 + 2ab3 + 9a5 − 7b5)x2 +(7a2 − 6b2 + 5a5 + 6a4b + 4ab4)x + 8 + 7b + 4b3

q1 = x10 + (7a − 6b + 3a3b2)x9 + (6ab3 + 9b4 − 7ab4)x8 + (6a3 − 5a3b − 5a4b − 8a2b3)x7 + (7a3b + b5)x6 +(−a − 5b2 + 7a4 − 3a5)x5 + (8b − 7a3 + 2a2b − 4b3 − 6ab3 + 3a4b)x4 + (−3b − 9a2 − 5ab3 − 4b4)x3 +(a5 + 9a4b − 3a2b3)x2 + (−8b − 3a2 − 7ab)x − 1 + 7a + b3

P2 (8, 7, 3, 5) :p2 = x8 + (b2c3 + 6bc4)x6 + (−2b4 + 2b3c + 3b3c2)x4 + (4b2c2 − 5a2c3 − 6abc3)x2 + (8a2b2c + 2b2c3)x + 5a2bc

q2 = x7 + (9ab + 2a2c2 − 9c5)x6 + (−8b3c + 5b5 + 6b4c)x4 + (−6ab2 − 3a4)x + 6a3c2

P3 (9, 8, 3, 5) :p3 = x9 + (6a − 8ab + c2 + a5)x8 + (a4b − 2a3b2)x7 + (7b3c + 8a3bc + 7b4c)x6 + (7a + 7a2 − 3c3)x5 +

(4c4 + 7a4b)x4 + (6bc + 9a2b2 − 7bc3)x3 + (5ac − 9a3bc)x2 + (5a2b2c + 8ab4)x − ab2c

q3 = x8 + (5abc2 − 5abc3)x7 + (b2c + 8b2c2 + 7bc4)x6 + (3a4 − 8a3b − 2a2b3)x5 + (2c4 + a2bc2)x4 +(8ab2 + 8b3 + 4ac4)x3 + (−5a2c + 4a5 + 4a4b + 6b4c)x − 5ab2 + 7abc

P4 (9, 9, 3, 5) :p4 = x9 + (7b3c + 8a3bc + 7b4c)x6 + (4c4 + 7a4b)x4 + (5ac − 9a3bc)x2 − ab2c

q4 = x8 + (5abc2 − 5abc3)x7 + (2c4 + a2bc2)x4 + (−5a2c + 4a5 + 4a4b + 6b4c)x − 5ab2

P5 (8, 7, 4, 5) :p5 = x8 + (ab + abc + 8acd2 + 2a4b)x6 + (ac2 − 2abc2 + 9bc2d + 8b2c3 − 3bd4)x5 + (d2 + 8b2c + 4bc2 −

6c2d + 5abd2 − 9a4c)x4 + 6ad + 8a3

q5 = x7 + (4cd3 + 2a2bc2)x5 + (7a4b + 9a3c2 − b5)x2 + (4d4 + 6a4b + 5a3bd + 5ac3d)x + 8a + ab +4a2d − 7ad2 + abc2 + 5a4c

18

P6 (8, 7, 4, 5) :p6 = x8 + (6b2c2d − 2b2cd2)x7 + (3a3bc + 8b2c3)x6 + (3adbc2 − 2b3c2)x5 + (6a2bc2 + 7ab2cd + 3b2cd2)x4 +

(4a2c2 + 9bcd2)x2 + (9a2bd2 − 2bc3d)x − 8abcd2

q6 = x7 +(4abcd+3a3bd)x6 +(adbc2 +8b3d2)x4 +(2adbc2 +ac4)x3 +(5abcd2−4bc2d2)x2 +(a2bc+7abc3)x+5abc

P7 (14, 6, 2, 5) :p7 = x14 + (b3 + 4ab3 + 6b4)x11 + (−2 + 6a + 4a2 − 4ab2)x10 + (a2b3)x10 + (5a + 6ab + 6a2b2 − 4a5 − 4ab4)x9 +

(9a3 − 8a2b + 5b3)x7 + (a3b − 6a2b2 − 3ab3 + 6a2b3)x2

q7 = x6 + (3a2b2 − 4b4 + 6a5)x4 + (a − 4a2 + 3a2b)x3 + (−8a − 9b + 8a3 − 8a3b − 5a2b2 − 4ab4)x2 +(3b4 + 2a3b2 + 3a2b3)x + 5 + 9a − 7a3 + 7a3b + 9ab3

P8 (14, 6, 2, 5) :p8 = x14 + (−8a + 5b2 − 9a2b + 9b3 + 7a2b2 + 8a5)x13 + (−8a − 2ab − 3a3 + 2a2b2 + b4 − 7a3b2)x12 +

(−5 − 4a3 − 4a2b − 6b3 − 5a2b2 − 3b5)x11 + (7 − 9ab2 − b3 − 2a4 − 2b4 + 5a3b2)x10 +(b − 8a2b + 3ab2 − 9b3 − 9a5 − a2b3)x9 + (−8 − 8ab + 9b3 − 2a3b − 3a5 − 9a2b3)x8 +(9b − 3a2 − 5b3 + 4a2b2)x7 + (8a − 2ab3 − 3a4b − 6a3b2 + 4ab4)x6 +(2b − 3a4 + 7a4b − 8a3b2 − 3ab4)x5 + (5ab − 9b2 + 7a4 + 6b4 − a4b + 4a3b2)x4 +(5b + 7a3 − 2b3 + 5a4 + a2b2 + 6b5)x3 + (−8 − 5a2 + 6b2 + 7b3 − a3b + 3ab4)x2 +(2 + 8a3 − 3a2b + 8a3b − 4b4 + 6b5)x + 3b + 4a2 + 5a2b2 + 6a4b

q8 = x6 + (5a2 − 8ab + 6b2 − 5ab3 + 4a2b3 − 4b5)x5 + (8a − 5b2 − 2ab2 + 2a3b + 9b4 + 8b5)x4 +(−6b2 − 2a3b + 6b4 + 3a5 + 8a4b + 6ab4)x3 + (4 − 9a2b + ab2 − 5a2b2 + 8a3b2 − 6a2b3)x2 +(b − 9ab + 4a4 + 7a2b2 + 9a4b + 5a3b2)x + 6ab − 7b2 − 4a2b3 − 2b5

P9 (8, 5, 3, 5) :p9 = x8+(−8bc+9a3b−2a2b2)x7+(−4+2a+b2c−9ab3−9ab3c)x6+(−5ac+8ab3−3b3c+7b2c2−7c4+2ab4)x5+

(5a− 9a2b + 3b2c)x4 + (3b2 + 7c2 + 3bc3 + 6a5)x3 + (c + 9c3 − a2c2 + 5ac3)x2 + (3b2 + 6bc + a2b + 4b3)x +4bc3 + 3b3c2 + 5c5

q9 = x5+(6b2c+6bc2+5a3b)x4+(−8a2b−8a4c−5a2b2c+3b5)x3+(−8ac+2ac3)x2+(−3b+5ab2−3c4−9a3b2)x+4ac2 − 9abc2 − 5a2c3

P10 (8, 5, 3, 7) :p10 = x8 + (−a3 − 9ac2)x5 + (b4 + 7b2c2 + 6a3bc + 4b5)x4 + (−7ab + 2bc− 5a2b − 6a3c + abc2 − 7c5)x3 +

7a2b2 + 5a2c2 − a3c2

q10 = x5 + (−2a − 4a2bc + 4a2b5)x3 + (8a2c − 4a3c − 7abc2 + 7b4 − 9a3b2 + 3a2c3)x2 + 9 − 2a2c2 − 3b3c

P11 (9, 4, 3, 5) :p11 = x9 + (6b2c2 − bc3)x8 + (7ab3c − 2ac4 + 6b5)x7 + (6a3bc)x6 + (5ab3 + 2a4b − 5a2bc2 − 8ab3c)x5 +

(a3c2 + 5a2b3 + 6a2bc2)x4 + (77ab2c + abc2)x3 + (8bc3 + 3a3c2 − 7b5)x2 + (b3 − 9a3b)x +9c2 − 3abc + 6bc4

q11 = x4 + 2ab2cx3 + (4a3c + 4a5 + 6b4c)x2 + (7bc + 2ab4 + 9ab2c2)x + 6c4 + 7a5 − 6abc3

P12 (10, 6, 3, 5) :p12 = x10 + (7ab + 4a2c + 5ac2 + 5a3b2)x8 + (5b2 − 8c2 − 3bc3 − 6b4c)x2 + 2bc + c2 − 2a2c − abc

q12 = x6+(−5b+4a2+6bc)x5+(4ab2c+3b4+5a4b)x2+(5ab2+7a4−5abc3−5bc4)x+4b+8c−8a3b+8b2c2+7bc4

P13 (8, 5, 4, 5) :p13 = x8 + (ab + abc + 8acd2 + 2a4b)x6 + (ac2 − 2abc2 + 9bc2d + 8b2c3 − 3bd4)x5 +

(d2 + 8b2c + 4bc2 − 6c2d + 5abd2 − 9a4c)x4 + 6ad + 8a3 + 3ab2d + 3b2

q13 = x5 + (ac2d2 + 3c4d)x3 + (a2cd2 − 6ac2d2 − 2ad4)x2 + (2abc2 − 3bc3 − c2d2 − 4a4b + 8acd3)x +7b3 − 9ac3

P14 (8, 5, 4, 6) :p14 = x8 + (3b2d2 − 6a3bcd − 5a2b2d)x7 + (3a2bdc + 6a3b2 + 7c2d3)x6 + (4b2d + 2bcd − 8a2c2)x5 +

(7b3c2 − 2b2c2d + 4bc3d)x4 + (ab2d2 − 9c3d2)x3 + (3b2d − 5c2d2 − 4a3c2)x2 +(c2d2 + 5a2b3 + 9abcd2)x + 8bcd2 − c2d2 − 4ac2d2

q14 = x5 + (5c3 + 6d3 − 9ab2cd)x4 + (7a2b3 + a2bd2 − 6ab3c)x3 + (9a3bc − 9abc2d − 3b2d3)x2 +(−2a3d − 9a2bcd + ab2c2)x − bcd2 − 2c2d2 + 6abd3

19

4.3. Analisis de los tiempos obtenidos

A la hora de elegir los polinomios, hemos tenido en cuenta los siguientes parametros:

Los grados de los polinomios considerados.

El numero de parametros en los coeficientes.

La densidad de los polinomios.

La diferencia entre el grado de A(x), n, y de B(x), m. Observese que para los pares P1, . . . , P6, tenemosn − m = 1 y para el resto, n − m > 1.

Estos parametros influyen a la hora de examinar los tiempos de la siguiente manera:

Los grados de los polinomios considerados: Esto condiciona la talla de las matrices. Observese las matricesdefinidas a partir de la matriz de Bezout Hıbrida son de orden n, mientras que las definidas a partir de lamatriz de Sylvester son de talla n + m − 2i, i = 0, . . . , m − 1.

El numero de parametros en los coeficientes: mayor numero de parametros, tiempos peores en general. Seha observado, por ejemplo, que cuando el numero de parametros es 4, los tiempos obtenidos usando lamatriz de Bezout son mejores.

La densidad de los polinomios. Esto puede explicar la diferencia de tiempos entre el par P7 y P8, porejemplo.

La diferencia n−m. Observese que las ultimas n−m filas de la matriz de Bezout Hıbrida vienen definidaspor los coeficientes del polinomio B(x). Por ello, a medida que esta diferencia aumenta, el comportamientode la matriz de Bezout Hıbrida deberıa mejorar e incluso ser mejor que el de la matriz de Sylvester. Esto,junto la talla de las matrices podrıa explicar el tiempo del par P12

Aun ası, hay resultados sorprendentes de difıcil explicacion. Por ejemplo, los tiempos del par P10.Otro ejemplo son los tiempos obtenidos aplicando Berkowitz, bien en paralelo, bien en secuencial. Se ha

comprobado como la implementacion paralela no mejora los tiempos de calculo utilizando la matriz de Sylvester.Sin embargo, sı hay un descenso al utilizar la matriz Bezout Hıbrida.

Ademas los tiempos obtenidos por la orden Determinant no deja de ser sorprendentes. Obviamente, situvieramos acceso a los codigos de Maple, los entenderıamos mejor. Este es un de los inconvenientes de trabajarcon un programa comercial.

Otra consecuencia del estudio experimental es tener constancia de la importancia de los valores asignadosa la maquina virtual Java de Maple. Si se aumentan los valores por defecto, es posible evitar los errores deejecucion por falta de memoria.

20

5. Trabajo futuro

Vista la influencia de Java en el rendimiento de la herramienta Maple tenemos pensado analizar con mayordetalle la gestion de memoria de la maquina virtual Java para conseguir un rendimiento adecuado con otrasbaterıas de polinomios. Es un campo de gran interes pues el lenguaje Java tambieen lo usan otras herramientasporque es multiplataforma (Windows, Linux, Mac, etc.) y por su versatilidad como middleware en entornosparalelos.

Ademas tenemos planteado escribir un artıculo sobre el calculo de Subresultantes aplicando la imple-mentacion del algoritmo de Berkowitz presentada en la memoria. Sin embargo, para dicho calculo, no se vaa aplicar tan solo la metodologıa descrita en la seccion 3, sino tambien el algoritmo presentado en [15], gener-alizacion de [3].

En [15], el autor presenta un algoritmo para el calculo de Subresultantes que consiste en lo siguiente: haciendooperaciones elementales y eliminando ciertas filas y columnas de la matriz de Bezout Hıbrida, se definen una seriede matrices, cuyos menores principales dominantes son coeficientes de las Subresultantes. Al aplicar el algoritmode Berkowitz a dichas matrices, se van calculando todos los menores principales dominantes en el proceso. Porello, el autor propone el algoritmo de Berkowitz cuando los polinomios tengan coeficientes polinomiales.

Nuestra idea es implementar dicho algoritmo, calcular tiempos y compararlos con las formulas determinan-tales de la seccion 3.

21

A. Anexo A: Mathematica vs Maple

A.1. Motivacion.

En la asignatura Herramientas Informaticas para Computacion Cientıfica se han visto distintas aplicacionescomo Matlab, Mathematica y Maple. Con estas herramientas se pueden abordar distintos trabajos de investi-gacion, como por ejemplo el de esta tesis de master, donde es preciso trabajar con polinomios. Como un primerpaso fundamental es la eleccion de la herramienta con la cual vamos a trabajar con calculo simbolico, es fun-damental decidir si vamos a emplear Mathematica o Maple. Para ayudar en la decision se ha decidido realizaruna serie de calculos de determinantes para tener un criterio en la eleccion de la herramienta.

A.2. Calculo de determinates.

Maple y Mathematica representan las matrices de forma parecida, aunque con caracteres delimitadoresdistintos. Ası, en Mathematica podemos tener esta matriz:

m1 = {{a6, a5, a4, a3, a2, a1, 0, 0, 0, 0, 0, 0, 0},

{0, a6, 0, a5, a4, 0, a3, a2, a1, 0, 0, 0, 0},

{a6, 0, a5, a4, 0, a3, a2, a1, 0, 0, 0, 0, 0},

{0, 0, a6, 0, 0, a5, a4, 0, 0, a3, a2, a1, 0},

{0, 0, 0, 0, a6, 0, 0, a5, a4, 0, a3, a2, a1},

{0, 0, b6, 0, 0, b5, b4, 0, 0, b3, b2, b1, 0},

{0, 0, 0, b6, 0, 0, b5, b4, 0, 0, b3, b2, b1},

{0, b6, 0, b5, b4, 0, b3, b2, b1, 0, 0, 0, 0},

{0, 0, 0, 0, b6, 0, 0, b5, b4, 0, b3, b2, b1},

{0, 0, 0, 0, c6, 0, 0, c5, c4, 0, c3, c2, c1},

{0, c6, 0, c5, c4, 0, c3, c2, c1, 0, 0, 0, 0},

{c6, 0, c5, c4, 0, c3, c2, c1, 0, 0, 0, 0, 0},

{0, 0, 0, 0, c6, 0, 0, c5, c4, 0, c3, c2, c1}}

y en Maple se representa de este otro modo:

m1 := Matrix([[a6, a5, a4, a3, a2, a1, 0, 0, 0, 0, 0, 0, 0],

[0, a6, 0, a5, a4, 0, a3, a2, a1, 0, 0, 0, 0],

[a6, 0, a5, a4, 0, a3, a2, a1, 0, 0, 0, 0, 0],

[0, 0, a6, 0, 0, a5, a4, 0, 0, a3, a2, a1, 0],

[0, 0, 0, 0, a6, 0, 0, a5, a4, 0, a3, a2, a1],

[0, 0, b6, 0, 0, b5, b4, 0, 0, b3, b2, b1, 0],

[0, 0, 0, b6, 0, 0, b5, b4, 0, 0, b3, b2, b1],

[0, b6, 0, b5, b4, 0, b3, b2, b1, 0, 0, 0, 0],

[0, 0, 0, 0, b6, 0, 0, b5, b4, 0, b3, b2, b1],

[0, 0, 0, 0, c6, 0, 0, c5, c4, 0, c3, c2, c1],

[0, c6, 0, c5, c4, 0, c3, c2, c1, 0, 0, 0, 0],

[c6, 0, c5, c4, 0, c3, c2, c1, 0, 0, 0, 0, 0],

[0, 0, 0, 0, c6, 0, 0, c5, c4, 0, c3, c2, c1]]):

Como podemos apreciar, mientras Mathematica usa como delimitador las llaves, en Maple se empleanlos corchetes. Otro elemento diferenciador es como se construye la matriz. Ası, en Mathematica se escribedirectamente la matriz. Sin embargo Maple necesita la funcion Matrix() para crear una matriz.

Otra posibilidad es la conversion de matrices entre herramientas. Por ejemplo Maple permite exportar lamatriz de este modo:

ExportMatrix(matriz, m1, target=delimited,delimiter=,format=rectangular);

Para importarla desde Mathematica empleamos:

23

m1 = Import[”matriz”,”Table”]

Es el momento de empezar a comparar ambas herramientas. Si calculamos el determinante de la matrizanterior, los tiempos obtenidos han sido:

Mathematica 6 Maple 11Tiempo 52,608 seg 0,092 seg

La toma de tiempos se ha realizado en un ordenador tipo PC con un procesador Intel Core2 Duo a 2,33 GHzcon 2 GB de memoria RAM con las versiones Mathematica 6.0 y Maple 11 para Linux de 32 bits.

Como se aprecia en la tabla, Maple es mas rapido en el calculo del determinante anterior.

A.3. La resultante de dos polinomios.

Tanto en Maple como en Mathematica es posible calcular la resultante de dos polinomios P y Q respecto auna indeterminada x de varias formas posibles. Veamos tres de ellas:

Con la funcion correspondiente.

Calculando el determinante de la matriz Sylvester.

Calculando el determinante de la matriz de Bezout.

Funcion resultante En Mathematica la funcion disponible es Resultant[] y en Maple es resultant().Veamos un ejemplo de uso en Mathematica:

p2 = x^8 + (b^2 c^3 + 6 b c^4) x^6 + (-2 b^4 + 2 b^3 c + 3 b^3 c^2) x^4 +

(4 b^2 c^2 - 5 a^2 c^3 - 6 a b c^3) x^2 + (8 a^2 b^2 c + 2 b^2 c^3) x + 5 a^2 b c

q2 = x^7 + (9 a b + 2 a^2 c^2 - 9 c^5) x^6 + (-8 b^3 c + 5 b^5 + 6 b^4 c) x^4 +

(-6 a b^2 - 3 a^4) x + 6 a^3 c^2

Resultant[p2,q2,x]

TimeUsed[]

El mismo codigo en Maple serıa:

t1 := time();

p2:=x^8+(b^2*c^3+6*b*c^4)*x^6+(-2*b^4+2*b^3*c+3*b^3*c^2)*x^4+

(4*b^2*c^2-5*a^2*c^3-6*a*b*c^3)*x^2+(8*a^2*b^2*c+2*b^2*c^3)*x+5*a^2*b*c;

q2:=x^7+(9*a*b+2*a^2*c^2-9*c^5)*x^6+(-8*b^3*c+5*b^5+6*b^4*c)*x^4+

(-6*a*b^2-3*a^4)*x+6*a^3*c^2;

resultant(p2,q2,x);

time()-t1;

Los tiempos obtenidos han sido:

Mathematica 6 Maple 11Tiempo 590,929 seg 23,669 seg

24

Resultante con matriz Sylvester En Mathematica, dentro de las opciones de la funcion Resultant[] pode-mos indicar que el metodo empleado sea la matriz de Sylvester.

p2 = x^8 + (b^2 c^3 + 6 b c^4) x^6 + (-2 b^4 + 2 b^3 c + 3 b^3 c^2) x^4 +


q2 = x^7 + (9 a b + 2 a^2 c^2 - 9 c^5) x^6 + (-8 b^3 c + 5 b^5 + 6 b^4 c) x^4 +

(-6 a b^2 - 3 a^4) x + 6 a^3 c^2

Resultant[p2,q2,x,Method->"SylvesterMatrix"]

TimeUsed[]

En el caso de Maple primero se crea la matriz de Sylvester con la funcion SylvesterMatrix() y despues secalcula su determinante:

t1 := time();

p2:=x^8+(b^2*c^3+6*b*c^4)*x^6+(-2*b^4+2*b^3*c+3*b^3*c^2)*x^4+


q2:=x^7+(9*a*b+2*a^2*c^2-9*c^5)*x^6+(-8*b^3*c+5*b^5+6*b^4*c)*x^4+

(-6*a*b^2-3*a^4)*x+6*a^3*c^2;

sm := LinearAlgebra:-SylvesterMatrix(p2,q2,x):

LinearAlgebra:-Determinant(sm);

time()-t1;

Los tiempos obtenidos han sido:

Mathematica MapleTiempo 590,966 seg 0,007 seg

Resultante con matriz Bezout Como en caso anterior, cambiamos la matriz Sylvester por la matriz deBezout. Para Mathematica el codigo es:

p2 = x^8 + (b^2 c^3 + 6 b c^4) x^6 + (-2 b^4 + 2 b^3 c + 3 b^3 c^2) x^4 +


q2 = x^7 + (9 a b + 2 a^2 c^2 - 9 c^5) x^6 + (-8 b^3 c + 5 b^5 + 6 b^4 c) x^4 +

(-6 a b^2 - 3 a^4) x + 6 a^3 c^2

Resultant[p2,q2,x,Method->"BezoutMatrix"]

TimeUsed[]

En el caso de Maple primero se crea la matriz de Bezout con la funcion BezoutMatrix() y despues secalcula su determinante:

t1 := time();

p2:=x^8+(b^2*c^3+6*b*c^4)*x^6+(-2*b^4+2*b^3*c+3*b^3*c^2)*x^4+


q2:=x^7+(9*a*b+2*a^2*c^2-9*c^5)*x^6+(-8*b^3*c+5*b^5+6*b^4*c)*x^4+

(-6*a*b^2-3*a^4)*x+6*a^3*c^2;

bm := LinearAlgebra:-BezoutMatrix(p2,q2,x):

LinearAlgebra:-Determinant(bm);

time()-t1;

La medicion de tiempos ha dado como resultado:

Mathematica MapleTiempo 98,964 seg 0,008 seg

25

A.4. Analisis y conclusiones.

Si reunimos todos los tiempos anteriores en una tabla

Mathematica 6 Maple 11Determinante 52,608 seg 0,092 segResultante 590,929 seg 23,669 segDeterminante Sylvester 590,966 seg 0,007 segDeterminante Bezout 98,964 seg 0,008 seg

puede apreciarse como Maple obtiene mejores tiempos. En consecuencia se ha decidio emplear Maple comoherramienta para el desarrollo de esta tesis de master.

26

B. Anexo B: OpenMaple API

B.1. Funciones especıficas de OpenMaple

Permiten iniciar una sesion maple, reiniciarla y pararla. Ademas de controlar la entrada y salida estandar.

Funcion Cometido EjemploStartMaple Inicia una sesion OpenMaple MKernelVector kv;

kv=StartMaple(argc,argv,&cb,NULL,NULL,err);char err[2048];

RestartMaple Limpia la memoria interna de Maple MKernelVector kv;RestartMaple(kv,err)

StopMaple Finaliza una sesion OpenMaple StopMaple(kv);

Cuadro 4: Funciones especıficas de OpenMaple.

El argumento &cb de StartMaple es una estructura de este tipo:

typedef struct {

void (M_DECL *textCallBack)(void *data, int tag, char *output);

void (M_DECL *errorCallBack)(void *data, M_INT offset, char *msg);

void (M_DECL *statusCallBack)(void *data, long kilobytesUsed,

long kilobytesAlloc, double cpuTime);

char * (M_DECL *readLineCallBack)(void *data, M_BOOL debug);

M_BOOL (M_DECL *redirectCallBack)(void *data, char *name, char *mode);

char * (M_DECL *streamCallBack)(void *data, char *name,

M_INT nargs, char **args);

M_BOOL (M_DECL *queryInterrupt)(void *data);

char * (M_DECL *callBackCallBack)(void *data, char *output);

} MCallBackVector, *MCallBack;

Si se desea el valor por defecto se emplea NULL. Sin embargo Maple aconseja definir al menos la funciontextCallBack que gestiona la salida de OpenMaple. Una posible definicion de textCallBack es:

/* callback used for directing result output */

static void M_DECL textCallBack(void *data, int tag, char *output)

{

printf("%s\n",output);

}

El resto de funciones pueden consultarse en la ayuda de Maple.

B.2. Funciones de conversion de Maple a estructuras de datos en C

Funcion Cometido EjemploMapleToM INT Convierte a un entero de la maquina i = MapleToM INT(kv,n);MapleToM BOOL Convierte a un booleno cbool = MapleToM BOOL(kv,mbool);MapleToString Convierte a string cstring = MapleToString(kv,mstring);

Cuadro 5: Funciones de conversion de Maple a C.

Otras funciones disponibles son MapleToComplexFloat32, MapleToComplexFloat64, MapleToFloat32,MapleToFloat64, MapleToInteger32, MapleToInteger64 y MapleToPointer.

27

Funcion Cometido EjemploM INT i;char expr[256];ALGEB var;

ToMapleName Crea una printf(“Enter variable: ”); fgets(expr,sizeof(expr),stdin);variable for(i = 0; i ¡strlen(expr); i++)Maple ;

expr[i-1] = ’\0’;var = ToMapleName(kv,expr,TRUE);

ToMapleInteger Convierte a un i = 12;entero Maple var = ToMapleInteger(kv,i);

ToMapleFloat Convierte a un double j=12.34;decimal Maple var = ToMapleFloat(kv,j);Convierte a double re, img;

ToMapleComplex un complejo re = 1; img = 5;Maple var = ToMapleComplex(kv,re,img; /* 1.+5.*I */Convierte a i =-1; var = ToMapleBoolean(kv,i); /* FAIL */

ToMapleBoolean un booleano i = 0; var = ToMapleBoolean(kv,i); /* false */Maple i = 1; var = ToMapleBoolean(kv,i); /* true */

ALGEB lista;ToMapleExpressionSequence Convierte a lista = ToMapleExpressionSequence(kv,2,

una secuencia ToMapleInteger(kv,1),ToMapleInteger(kv,2)); /* [1,2] */ALGEB rel;

ToMapleRelation Crea una rel = ToMapleRelation(kv,“>”,relacion ToMapleInteger(kv,2),ToMapleInteger(kv,1)); /* 1 < 2 */

ALGEB f, result;ToMapleFunction Crea una i = 0; f = EvalMapleStatement(kv,“cos:”);

funcion f = ToMapleFunction(kv,f,1,ToMapleInteger(kv,i));result = MapleEval(kv,f); /* sin(0)=1 */

Cuadro 6: Funciones de conversion de C a Maple.

B.3. Funciones de conversion de C a objetos Maple

Otras funciones disponibles son ToMapleChar, ToMapleComplexFloat, ToMapleInteger64, NewMapleEx-pressionSequence, ToMapleNULL, ToMapleNULLPointer, ToMaplePointer, ToMapleString y ToMapleUneval.

B.4. Funciones para mostrar objetos Maple

Funcion Cometido EjemploMaplePrintf Muestra datos en C MaplePrintf(kv,”var=”%d”,MapleToFloat32(kv,var));MapleALGEB Printf Muestra objetos Maple MapleALGEB Printf(kv, ”var=%a”, var);

Cuadro 7: Funciones de salida de datos Maple.

Otras funciones disponibles son MapleALGEB SPrintf y MapleUserInfo.

B.5. Funciones de comprobacion de tipo de dato Maple

Otras funciones disponibles son MapleNumArgs, IsMapleComplexNumeric, IsMapleNumeric, IsMapleInte-ger8, IsMapleInteger16, IsMapleInteger32, IsMapleInteger64, IsMapleNULL, IsMaplePointer, IsMaplePointer-

28

Funcion Cometido EjemploALGEB x;

IsMapleName Verifica si x = ToMapleName(kv,”x”,TRUE);existe el if(IsMapleName(kv,x))objeto Maple printf(“ Name OK”);

MapleAssign(kv,x,ToMapleInteger(kv,10));IsMapleAssignedName Verifica si if(IsMapleAssignedName(kv,x))

esta asignado printf(“ Assign OK”);M INT i;

IsMapleInteger Verifica si i = 10;es un entero if(!IsMapleInteger(kv,ToMapleInteger(kv,i)))

MapleRaiseError(kv,”integer expected”);ALGEB list;i = 2;

IsMapleList Verifica si list = MapleListAlloc(kv,i);el objeto es for(i = 1; i < 3; i++)una lista MapleListAssign(kv,x,i,ToMapleInteger(kv,i));

if(IsMapleList(kv,x))printf(“ List OK”);

ALGEB M; /* rtable like a matrix */M = EvalMapleStatement(kv,”ImportMatrix(matriz,

IsMapleRTable Verifica si source=MatrixMarket):”);es un rtable if(IsMapleRTable(kv,M))

printf(“ Rtable OK”);

Cuadro 8: Funciones de comprobacion de tipo de dato Maple.

NULL, IsMapleProcedure, IsMapleSet, IsMapleStop, IsMapleString, IsMapleTable, IsMapleUnassignedName yIsMapleUnnamedZero.

B.6. Funciones del sistema

Funcion Cometido EjemploALGEB inicio, fin;

Establece u inicio = MapleKernelOptions(kv,cputime”,NULL);MapleKernelOptions obtiene las calculo(); /* Calculo en Maple */

variables fin = MapleKernelOptions(kv,cputime”,NULL);del kernel MapleALGEB Printf(kv,”%a”,EvalMapleProc(kv,

ToMapleName(kv,”−”,TRUE),2,fin,inicio));

Cuadro 9: Funciones del sistema.

Otras funciones disponibles son MapleHelp y MapleLibName.

B.7. Funciones de manipulacion de tablas y listas

Otras funciones disponibles son MapleTableAlloc, MapleTableAssign, MapleTableDelete, MapleTableHasEn-try y MapleTableSelect.

B.8. Funciones de asignacion y seleccion

Otras funciones disponibles son MapleAssignIndexed, MapleSelectImaginaryPart y MapleSelectRealPart.

29

Funcion Cometido EjemploM INT n;

MapleListAlloc Crea una ALGEB list;lista n = 3;

list = MapleListAlloc(kv,n);MapleListAssign Asigna un MapleListAssign(kv,list,1,

valor ToMapleInteger(kv,1));MapleListSelect Selecciona n = MapleToM INT(kv,

un elemento MapleListSelect(kv,list,1);

Cuadro 10: Funciones de manipulacion de tablas y listas.

Funcion Cometido EjemploAsigna un ALGEB var;

MapleAssign valor a una var = ToMapleName(kv,”var”,TRUE);variable Maple MapleAssign(kv,var,ToMapleInteger(kv,100));

M INT i, n, index[1];ALGEB list, r;

Selecciona de n = 3; list = MapleListAlloc(kv,n);MapleSelectIndexed una tabla, for(i=1;i¡=n;i++)

lista o rtable MapleListAssign(kv,list,i,ToMapleInteger(kv,i));index[0]=2;r = MapleSelectIndexed(kv,list,1,index); /* 2 */

Cuadro 11: Funciones de asignacion y seleccioon.

B.9. Funciones de gestion de memoria y objetos

Funcion Cometido Ejemplochar expr[128];char *buff;ALGEB alloc, var;

MapleAlloc Reserva printf(“Enter variable: ”); fgets(expr,sizeof(expr),stdin);memoria alloc = MapleAlloc(kv,strlen(expr)*sizeof(char));

buff = (char*)alloc;sprintf(buff,” %s”,expr);var = ToMapleName(kv,buff,TRUE);

MapleDispose Borra reserva MapleDispose(kv,alloc);MapleGcIsProtected Mira si esta ALGEB rt;

protegido del gc /* If garbage icollection run, clean rt valueand the program crash */

MapleGcProtect Protege al if(!MapleGcIsProtected(kv,rt))objeto del gc MapleGcProtect(kv,rt);

Cuadro 12: Funciones de gestion de memoria y objetos.

Otras funciones disponibles son MapleGcAllow, MapleGcMark, MaplePointerSetDisposeFunction, Maple-PointerSetMarkFunction, MaplePointerSetPrintFunction, MaplePointerSetType y MaplePointerType.

30

Funcion Cometido EjemploEvalMapleStatement Evaua una ALGEB r, x;

cadena de r = EvalMapleStatement(kv,”int(1/(x4+1),x);”);texto en Maple MapleAssign(kv,ToMapleName(kv,”x”,TRUE),

MapleEval Evalua un ToMapleInteger(kv,0));objeto maple r = MapleEval(kv,r);

M INT n, bounds[4];ALGEB rt, f, ncol;RTableSettings rts;EvalMapleStatement(kv,”with(LinearAlgebra):”);n = 4;

Eval’ua un RTableGetDefaults(kv,&rts);EvalMapleProc objeto rts.num dimensions = n;

funcion rts.subtype = RTABLE MATRIX;bounds[0] = 1; bounds[1] = n;bounds[2] = 1; bounds[3] = n;rt = RTableCreate(kv,&rts,NULL,bounds);f = EvalMapleStatement(kv,ColumnDimension:”);ncol = EvalMapleProc(kv,f,1,rt);M INT i;

MapleRaiseError Control de i = 10;error if(!IsMapleInteger(kv,ToMapleInteger(kv,i)))

MapleRaiseError(kv,”integer expected”);M INT i;

MapleCheckInterrupt Gestiona for(i = 0; i < 5000000; i++)control-c MapleCheckInterrupt(kv);

Cuadro 13: Funciones de evaluacion y manejo de errores.

B.10. Funciones de evaluacion y manejo de errores

Otras funciones disponibles son DoubleTohfData, EvalhfDataProc, EvalhfMapleProc, ImaginaryhfData,MapleEvalhf, MapleRaiseError1, MapleRaiseError2, MapleTohfData, MapleTrapError, MapleUnique, Realhf-Data y ToMaplehfData.

B.11. Funciones de manipulacion de tablas rectangulares (rtables)

Otras funciones disponibles son RTableAppendIndFn, RTableCopy, RTableCopyImPart, RTableCopyReal-Part, RTableDataBlock, RTableIndFn, RTableIndFnArgs, RTableIsReal, RTableLowerBound, RTableNumDi-mensions, RTableSetIndFn, RTableSetType, RTableSparseCompact, RTableSparseIndexRow, RTableSpar-seIndexSort, RTableSparseResize, RTableSparseSetNumElems, RTableSparseSize, RTableUpperBound yRTableZipReIm.

31

Funcion Cometido EjemploRTableGetDefaults Valores RTableSettings rts;

por defecto RTableGetDefaults(kv,&rts);M INT n, bounds[4];INTEGER32 *data;ALGEB rt;RTableSettings rts;n = 4;RTableGetDefaults(kv,&rts);rts.num dimensions = 2;

RTableCreate Crea una rts.subtype = RTABLE MATRIX;tabla rts.data type = RTABLE INTEGER32;rectangular bounds[0] = 1; bounds[1] = n;

bounds[2] = 1; bounds[3] = n;rt = RTableCreate(kv,&rts,NULL,bounds);data = (INTEGER32*)RTableDataBlock(kv,rt);data[MATRIX OFFSET FORTRAN RECT(1,1,n,n)] = 1;MapleALGEB Printf(kv,”%a”,rt);/* Matrix(4,4,[[1,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],datatype=integer[4]) */n = 4;RTableGetDefaults(kv,&rts);

RTableSetAttribute Incluye un rts.num dimensions = 2;atributo rts.subtype = RTABLE MATRIX;

rts.data type = RTABLE INTEGER32;RTableSetAttribute(kv,&rts,”stack”);RTableAppendAttribute(kv,&rts,”12345”);bounds[0] = 1; bounds[1] = n;

RTableAppendAttribute Asigna un bounds[2] = 1; bounds[3] = n;valor al rt = RTableCreate(kv,&rts,NULL,bounds);atributo MapleALGEB Printf(kv,”%a”,rt);

/* Matrix(4,4,[[1,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],datatype=integer[4],attributes = [stack, ‘12345‘]) */

RTableGetSettings Propiedades /* rtable created by RTableCreate function */de la rtable RTableGetSettings(kv,&rts,rt);

M INT index[2];RTableSelect Recupera index[0] = 1; index[1] = 1;

un dato rtd = RTableSelect(kv,rt,index);MaplePrintf(kv,”%d”, rtd.dag); /* 1 */M INT index1[2], index2[2];RTableData rtd1, rtd2;INTEGER32 *data;data = (INTEGER32*)RTableDataBlock(kv,rt);data[MATRIX OFFSET FORTRAN RECT(1,1,n,n)] = 1;data[MATRIX OFFSET FORTRAN RECT(2,1,n,n)] = 2;

RTableAssign Asigna en /* [[1,0,0,0],[2,0,0,0],[0,0,0,0],[0,0,0,0]] */una rtable index1[0] = 1; index1[1] = 1;

rtd1 = RTableSelect(kv,rt,index1);index2[0] = 2; index2[1] = 1;rtd2 = RTableSelect(kv,rt,index2);RTableAssign(kv,rt,index1,rtd2);RTableAssign(kv,rt,index2,rtd1);/* [[2,0,0,0],[1,0,0,0],[0,0,0,0],[0,0,0,0]] */

RTableNumElements Numero de i = RTableNumElements(kv,rt);elementos /* i = 16 */

Cuadro 14: Funciones de manipulacion de rtables.32

C. Anexo C: Algoritmo de Berkowitz en OpenMaple

C.1. Codigo C del Algoritmo de Berkowitz en secuencial

#include <stdio.h>

#include <stdlib.h>

#include <signal.h>

#include <string.h>

#include <sys/types.h>

#include <time.h>

#include "maplec.h"

/* callback used for directing result output */

static void M_DECL textCallBack(void *data, int tag, char *output)

{

printf("%s\n",output);

}

int min(int a, int b)

{

if(a <= b)

return a;

else

return b;

}

/* berkoseq */

ALGEB M_DECL berkoseq(MKernelVector kv, ALGEB M, ALGEB X)

{

M_INT ncol; /* C column dimension of M */

M_INT nfil; /* C row dimension of M */

M_INT index1[2]; /* X,Y position on Matrix */

M_INT r,i,j,k; /* Contadores */

static M_BOOL was_protected = FALSE;

ALGEB f; /* Maple function */

ALGEB result; /* Result to apply function */

ALGEB n; /* Maple column dimension of M */

ALGEB rt; /* Rtable like matrix */

RTableData rtd; /* Rtable cell data */

ALGEB temp1,temp2; /* Maple temp values */

ALGEB temp3,temp4; /* Cell value from rtd */

ALGEB variable; /* Var input from user */

ALGEB Vect; /* Vector coeficientes pol caracteristico anterior */

ALGEB C;

ALGEB S;

ALGEB Q;

ALGEB Temp; /* Vector with temporal data */

/* Determine column dimension of M */

f = EvalMapleStatement(kv,"ColumnDimension:");

n = EvalMapleProc(kv,f,1,M);

ncol = MapleToInteger32(kv,n);

/* Determine row dimension of M */

f = EvalMapleStatement(kv,"RowDimension:");

33

nfil = MapleToInteger32(kv,EvalMapleProc(kv,f,1,M));

/* Check dimension */

if(ncol != nfil)

{

printf("Error: not a square matrix");

exit(1);

}

/* OpenMaple code for A[1,1] */

rt = M;

index1[0] = 1; index1[1] = 1;

rtd = RTableSelect(kv,rt,index1);

temp1 = ToMapleName(kv,"temp1",TRUE);

MapleAssign(kv,temp1,rtd.dag);

/* Prevent garbage collection If garbage icollection run,

clean rt value and the program crash */

if(!MapleGcIsProtected(kv,rt))

MapleGcProtect(kv,rt);

/* OpenMaple code for X */

variable = ToMapleName(kv,"variable",TRUE);

MapleAssign(kv,variable,X);

/* Trivial value 1x1 */

switch(ncol)

{

case 1 : /* OpenMaple code for A[1,1]-X */

return EvalMapleStatement(kv,"temp1-variable:");

case 2 : /* OpenMaple code for normal((A[1, 1]-X)*(A[2, 2]-X)-A[2, 1]*A[1, 2]) */

index1[0] = 2; index1[1] = 2;




index1[0] = 2; index1[1] = 1;




index1[0] = 1; index1[1] = 2;




return EvalMapleStatement(kv,"normal((temp1-variable)*(temp2-variable)-(temp3*temp4)):");

default : /* Berko code */

/* OpenMaple code for Vect:=table([1=-1, 2=A[1,1]]) */

Vect = MapleTableAlloc(kv);

MapleTableAssign(kv,Vect,ToMapleInteger(kv,1),ToMapleInteger(kv,-1));

MapleTableAssign(kv,Vect,ToMapleInteger(kv,2),temp1);

/* Prevent garbage collection */

34

if(!MapleGcIsProtected(kv,Vect))

MapleGcProtect(kv,Vect);

/* OpenMaple code for C[1]:=-1 */

C = MapleTableAlloc(kv);

MapleTableAssign(kv,C,ToMapleInteger(kv,1),ToMapleInteger(kv,-1));


if(!MapleGcIsProtected(kv,C))

MapleGcProtect(kv,C);

/* Other vars */

S = MapleTableAlloc(kv);

Q = MapleTableAlloc(kv);

Temp = MapleTableAlloc(kv);





if(!MapleGcIsProtected(kv,S))

MapleGcProtect(kv,S);

if(!MapleGcIsProtected(kv,Q))

MapleGcProtect(kv,Q);

if(!MapleGcIsProtected(kv,Temp))

MapleGcProtect(kv,Temp);

/* OpenMaple code: for r from 2 to n do */

for(r=2; r <= ncol; r++)

{

/* OpenMaple code: for i to r-1 do S[i] := A[i,r] end do; */

for(i=1; i < r; i++)

{

index1[0] = i; index1[1] = r;

rtd = RTableSelect(kv,rt,index1); /* A[i,r] */

MapleTableAssign(kv,S,ToMapleInteger(kv,i),rtd.dag); /* S[i]:=A[i,r] */

}

/* OpenMaple code: C[2] := A[r,r]; */

index1[0] = r; index1[1] = r;

rtd = RTableSelect(kv,rt,index1); /* A[r,r] */

MapleTableAssign(kv,C,ToMapleInteger(kv,2),rtd.dag); /* C[2]:=A[r,r] */

/* OpenMaple code: for i to r-2 do */

for(i=1; i<r-1; i++)

{

/* OpenMaple code: C[i+2] := normal(add(A[r,j]*S[j], j=1 .. r-1)) */

for(j=1; j<r; j++)

{

/* First, product */

index1[0] = r; index1[1] = j;

rtd = RTableSelect(kv,rt,index1); /* A[r,j] */


MapleAssign(kv,temp3,MapleTableSelect(kv,S,ToMapleInteger(kv,j))); /* S[j] */

35

result = EvalMapleStatement(kv,"temp2*temp3:");

if(j==1)

MapleTableAssign(kv,Temp,ToMapleInteger(kv,j),result);

else

{

MapleAssign(kv,temp3,result);

MapleAssign(kv,temp2,MapleTableSelect(kv,Temp,ToMapleInteger(kv,j-1)));

/* Second, add */

result = EvalMapleStatement(kv,"temp2+temp3:");


}

} /* end for(j=1; j<r; j++) */

/* And finally, normal */

f = EvalMapleStatement(kv,"normal:");

result = EvalMapleProc(kv,f,1,MapleTableSelect(kv,Temp,ToMapleInteger(kv,j-1)));

/* Assign to C[i+2] */

MapleTableAssign(kv,C,ToMapleInteger(kv,i+2),result);

/* OpenMaple code: for j to r-1 do */

for(j=1; j<r; j++)

{

/* OpenMaple code: Q[j] := normal(add(A[j,k]*S[k], k=1 .. r-1)) */

for(k=1; k<r; k++)

{


index1[0] = j; index1[1] = k;

rtd = RTableSelect(kv,rt,index1); /* A[j,k] */


MapleAssign(kv,temp3,MapleTableSelect(kv,S,ToMapleInteger(kv,k))); /* S[k] */


if(k==1)

MapleTableAssign(kv,Temp,ToMapleInteger(kv,k),result);

else

{


MapleAssign(kv,temp2,MapleTableSelect(kv,Temp,ToMapleInteger(kv,k-1)));

/* Second, add */


MapleTableAssign(kv,Temp,ToMapleInteger(kv,k),result);

}

} /* end for(k=1; k<r; k++) */



result = EvalMapleProc(kv,f,1,MapleTableSelect(kv,Temp,ToMapleInteger(kv,k-1)));

/* Assign to Q[j] */

MapleTableAssign(kv,Q,ToMapleInteger(kv,j),result);

} /* end for(j=1; j<r; j++) */

/* OpenMaple code: for j to r-1 do S[j] := Q[j] end do */

for(j=1; j<r; j++)

MapleTableAssign(kv,S,ToMapleInteger(kv,j),

MapleTableSelect(kv,Q,ToMapleInteger(kv,j)));

} /* end for(i=1; i<r-1; i++) */

36

/* OpenMaple code: C[r+1] := normal(add(A[r, j]*S[j], j = 1 .. r-1)) */

for(j=1; j<r; j++)

{


index1[0] = r; index1[1] = j;

rtd = RTableSelect(kv,rt,index1); /* A[r,j] */


MapleAssign(kv,temp3,MapleTableSelect(kv,S,ToMapleInteger(kv,j))); /* S[j] */


if(j==1)


else

{



/* Second, add */



}

} /* end for(j=1; j<r; j++) */




/* Assign to C[r+1] */

MapleTableAssign(kv,C,ToMapleInteger(kv,r+1),result);

/* OpenMaple code: for i to r+1 do

Q[i] := normal(add(C[i+1-j]*Vect[j], j = 1 .. min(r, i)))

end do; */

for(i=1; i<=r+1; i++)

{

for(j=1; j<=min(r,i); j++)

{


MapleAssign(kv,temp2,MapleTableSelect(kv,C,ToMapleInteger(kv,i+1-j)));

MapleAssign(kv,temp3,MapleTableSelect(kv,Vect,ToMapleInteger(kv,j)));


if(j==1)


else

{



/* Second, add */



}

} /* end for(j=1; j<=min(r,i); j++) */




/* Assign to Q[i] */

MapleTableAssign(kv,Q,ToMapleInteger(kv,i),result);

} /* end for(i=1; i<=r+1; i++) */

37

/* OpenMaple code: for i to r+1 do Vect[i]:=Q[i] end do */

for(i=1; i<=r+1; i++)

MapleTableAssign(kv,Vect,ToMapleInteger(kv,i),

MapleTableSelect(kv,Q,ToMapleInteger(kv,i)));

} /* end for(r=2; r <= ncol; r++) */

} /* switch */

/* OpenMaple code: add(Vect[i+1]*X^(n-i), i = 0 .. n) */

for(i=0; i<=ncol ; i++)

{


MapleAssign(kv,temp2,MapleTableSelect(kv,Vect,ToMapleInteger(kv,i+1))); /* Vect[i+1] */

MapleAssign(kv,temp1,ToMapleInteger(kv,ncol-i));

MapleAssign(kv,temp3,EvalMapleStatement(kv,"variable^temp1:")); /* X^(n-i) */


if(i==0)

MapleTableAssign(kv,Temp,ToMapleInteger(kv,i),result);

else

{


MapleAssign(kv,temp2,MapleTableSelect(kv,Temp,ToMapleInteger(kv,i-1)));

/* Second, add */


MapleTableAssign(kv,Temp,ToMapleInteger(kv,i),result);

}

}

return result;

}

int main(int argc, char *argv[])

{

M_INT i;

static M_BOOL was_protected = FALSE;

char expr[2048], err[2048]; /* command input and error string buffers */

MKernelVector kv; /* Maple kernel handle */

MCallBackVectorDesc cb = { textCallBack,

0, /* errorCallBack not used */

0, /* statusCallBack not used */

0, /* readLineCallBack not used */

0, /* redirectCallBack not used */

0, /* streamCallBack not used */

0, /* queryInterrupt not used */

0 /* callBackCallBack not used */

};

ALGEB var;

ALGEB A; /* Matrix */

ALGEB result;

/* initialize Maple */

if((kv=StartMaple(argc,argv,&cb,NULL,NULL,err)) == NULL)

{

printf("Fatal error, %s\n",err);

return(1);

38

}

EvalMapleStatement(kv,"with(LinearAlgebra):");

printf("Enter variable: ");

fgets(expr,sizeof(expr),stdin);

for(i=0;i<strlen(expr);i++)

;

expr[i-1]=’\0’;

var = ToMapleName(kv,expr,TRUE);

/* Read matrix from file matriz */

printf("Leyendo la matriz del fichero matriz...\n");

A = EvalMapleStatement(kv,"ImportMatrix(matriz,source=MatrixMarket):");

/* Prevent garbage collection. If garbage icollection run,

clean A values and the program crash */

if(!MapleGcIsProtected(kv,A))

MapleGcProtect(kv,A);

result = berkoseq(kv,A,var);

MapleALGEB_Printf(kv,"El polinomio caracteristico es %a\n",result);

StopMaple(kv);

return(0);

}

C.2. Compilacion y ejecucion del codigo C

Suponiendo que tenemos la matriz polinomial en un fichero denominado matriz, el fichero con el codigoOpenMaple se denomina berkoseq.c y tenemos instalado Maple 11 en el directorio /usr/local/maple11, entonceslas ordenes necesarias para compilar y ejecutar el programa son:

export MAPLE=/usr/local/maple11

export PATH=$PATH:/usr/local/maple11/bin

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/maple11/bin.IBM_INTEL_LINUX

gcc -I/usr/local/maple11/extern/include -L/usr/local/maple11/bin.IBM_INTEL_LINUX \

-lmaplec -lc -lm -o berkoseq berkoseq.c

./berkoseq

39

D. Anexo D: Implementacion OpenMaple de Berkowitz vs Berkow-

itzAlgorithm de Maple

D.1. Tiempo para BerkowitzAlgorithm de Maple

Si tenemos una matriz polinomial cuadrada de dimension 14 en un fichero denominado matriz podemoscalcular el polinomio caracterıstico empleando el algoritmo de Berkowitz mediante estas ordenes Maple.

with(LinearAlgebra[Generic]):

Z[‘0‘], Z[‘1‘], Z[‘+‘], Z[‘-‘], Z[‘*‘], Z[‘=‘] := 0, 1, ‘+‘, ‘-‘, ‘*‘, ‘=‘:

A := ImportMatrix(matriz, source = MatrixMarket):

st := time():

C := BerkowitzAlgorithm[Z](A):

time()-st;

En un ordenador con un procesador Intel Core2 Duo a 2,33 GHz el tiempo invertido ha sido 0,04 segundos.

D.2. Tiempo para la implementacion OpenMaple de Berkowitz

Si tenemos una matriz polinomial cuadrada de dimension 14 en un fichero denominado matriz podemoscalcular el polinomio caracterıstico con nuestro programa escrito en OpenMaple tecleando lo siguiente:

export MAPLE=/usr/local/maple11

export PATH=$PATH:/usr/local/maple11/bin

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/maple11/bin.IBM_INTEL_LINUX

gcc -I/usr/local/maple11/extern/include -L/usr/local/maple11/bin.IBM_INTEL_LINUX \

-lmaplec -lc -lm -o berkoseq berkoseq.c

./berkoseq

En el mismo ordenador el tiempo ha sido 186,01 segundos.

D.3. Matriz polinomial empleada en la toma de tiempos

La matriz de Sylvester polinomial cuadrada de dimension catorce ha sido creada a partir de estos dospolinomios:

p1 = x11 + (2ab3 − 3a5 − a2b3)x10 + (8a4b − 5a2b3)x9 + (a3b2 + 3a2b3)x8 + (4a4b + 3a3b2)x7 +(8 − 6a + 7b + 8a2 − 5a3b2)x6 + (−9 − 7a − 8ab − 8a2b + 9b5)x5 + (7a2 + 8a3 + 4a2b + 8a4 +8a3b + 9b5)x4 + (3a + 3a2 − 5b3 − 9a3b − 5ab3 − 4ab4)x3 + (5ab2 + 2ab3 + 9a5 − 7b5)x2 +(7a2 − 6b2 + 5a5 + 6a4b + 4ab4)x + 8 + 7b + 4b3

q1 = x10 + (7a − 6b + 3a3b2)x9 + (6ab3 + 9b4 − 7ab4)x8 + (6a3 − 5a3b − 5a4b − 8a2b3)x7 + (7a3b + b5)x6 +(−a − 5b2 + 7a4 − 3a5)x5 + (8b − 7a3 + 2a2b − 4b3 − 6ab3 + 3a4b)x4 + (−3b − 9a2 − 5ab3 − 4b4)x3 +(a5 + 9a4b − 3a2b3)x2 + (−8b − 3a2 − 7ab)x − 1 + 7a + b3

41

E. Anexo E: Grid Computing Toolbox for Maple

E.1. Primer programa paralelo. Reparto de matrices Toeplitz entre los proce-sadores

El primer enfoque ha sido el modelo maestro-esclavo. Es decir, el procesador maestro (procesador 0) deter-mina el trabajo del resto de procesadores. El trabajo de los n-1 procesadores esclavos es calcular las matricesToeplitz. Una vez que el procesador 0 recibe todas las matrices Toeplitz realiza el producto de dichas matrices.Por tanto el calculo de matrices Toeplitz se realiza en paralelo en los n-1 procesadores y el producto en elprocesador maestro.

Suponiendo una matriz inicial de dimension 14x14 y nuestro procesador Intel Quad Core con cuatro nucleos,entonces el reparto de calculos de matrices Toeplitz entre los tres nucleos esclavos serıa:

Procesador Indice matriz Toeplitz1 2, 5, 8, 11, 142 3, 6, 9, 123 4, 7, 10, 13

Cuadro 15: Reparto calculo matrices Toeplitz entre procesadores esclavos.

Pasemos a ver el codigo Maple. En un principio se deseaba crear el grid desde un fichero Maple y despuesejecutar el codigo paralelo. Por este motivo se usan dos ficheros:

1. call berkopar1.mpl : crea el grid y ejecuta el programa berkopar1.mpl.

2. berkopar1.mpl : codigo paralelo con el reparto y multiplicacion de matrices Toeplitz.

Creacion del grid y ejecucion del programa paralelo. Veamos el contenido del ficherocall berkopar1.mpl.

# Procedimiento para leer un fichero de texto con comandos Maple

LoadCode := proc(n)

local sep, fname,r,l;

sep := kernelopts(dirsep);

fname := cat(kernelopts(homedir), sep, "maple", sep, "grid", sep, n, ".mpl");

r := "";

l := FileTools[Text][ReadLine](fname);

while l <> NULL do

r := cat(r, l, "\n");


end do;

FileTools[Text][Close](fname);

r;

end:

# Procedimiento para Launch()

checkAbort := proc()

false;

end proc;


fname := cat(kernelopts(homedir), sep, "maple", sep, "grid", sep, "matriz");

# Se importa la matriz desde el fichero $HOME/maple/grid/matriz

M := ImportMatrix(fname,source=MatrixMarket);

# Se carga el fichero berkopar1.mpl con el codigo Maple paralelo

43

Code := LoadCode(berkopar1);

# Valores previos para iniciar el grid

host := "localhost":

port := 2000:

numNodes := 4:

broadcastAddress := "127.0.0.255":

broadcastPort := 4400:

logFile := "logs/gridlog.txt":

# Se lanza el grid

Grid[Server][StartServer](port,numNodes,broadcastAddress,broadcastPort,logFile);

Grid[Setup](host,port);

# Se comprueba el numero de nodos disponibles en el grid

(TotalNodes, AvailableNodes) := Grid[Status]();

# Se ejecuta el programa paralelo, la matriz M se pasa a todos los procesadores

r := Grid[Launch](AvailableNodes, Code, printf, checkAbort, ["M"]):

# Se para el grid y termina el programa

Grid[Server][StopServer](port,numNodes);

Para ejecutar este fichero basta teclear:

maple < call_berko1.mpl

Codigo del primer programa paralelo. Veamos el contenido del fichero berko1.mpl.

# Usamos el paquete LinearAlgebra de Maple

with(LinearAlgebra);

# n

# Producto R x A x S

prod_ras := proc(Indice::integer, ExpoA::integer)

local i, j, k, pras, prxa, polcar, T, R, S;

if ExpoA = 0 then

# R x S

pras := normal(add(M[Indice,j]*M[j,Indice], j = 1..Indice-1));

else

# n

# R x A

for k from 1 to Indice-1 do

T[k] := M[Indice, k];

end do;

for k from 1 to ExpoA do

for i from 1 to Indice-1 do

prxa[i] := normal(add(T[j]*M[j,i], j = 1..Indice-1));

end do;

for j from 1 to Indice-1 do

T[j] := prxa[j];

end do;

end do;

# n

# R x A x S

pras := normal(add(prxa[j]*M[j,Indice], j = 1..Indice-1));

44

end if;

end proc:

# Procedimiento para el calculo del vector columna C y la matriz Toeplitz

calculotoep := proc(n::integer)

local i, C;

C := Vector(n+1);

C[1] := -1;

C[2] := M[n,n];

for i from 3 to n+1 do

C[i] := prod_ras(n, i-3);

end do;

C:

end proc:

# Procedimiento principal que distribuye el calculo de matrices Toeplitz

berkopar := proc(X::name)

local thisnode, totalnodes, i, j, k, n, Toeplitz, Vect, C, Q, msg, mensaje, posicion, calculo, polcar;

thisnode := Grid:-Util:-MyNode();

totalnodes := Grid:-Util:-NumNodes();

if totalnodes < 3 then error "Grid mus have at least 3 nodes" end if;

# El nodo cero lleva controla el calculo de matrices Toeplitz

if thisnode = 0 then

n := ColumnDimension(M);

# La matriz debe ser cuadrada

if RowDimension(M) <> n then error "Not a square matrix" end if;

# Casos particulares

if n = 1 then RETURN([1, 1] - X)

elif n = 2 then RETURN(normal((M[1, 1] - X)*(M[2, 2] - X) - M[2, 1] * M[1, 2])) end if;

# Array con las matrices Toeplitz

Toeplitz := Array(1..n);

# Toeplitz(Q1)

Toeplitz[1] := Vector([-1, M[1,1]]);

# Envio de listas de trabajos a los procesadores

k := trunc(n/(totalnodes-1));

for j from 2 to totalnodes do

if (totalnodes-1)*k+j > n then k := k-1 end if;

msg := [seq((totalnodes-1)*i+j, i=0..k)];

Grid:-Util:-Send(j-1, msg);

end do;

# Recibe de los procesadores Toeplitz(Q2) a Toeplitz(Qn)

for j from 2 to n do

mensaje := Grid:-Util:-Receive();

45

posicion := op(1, mensaje);

calculo := op(2, mensaje);

printf("\n Recibido Toeplitz(%a)", posicion);

Toeplitz[posicion] := calculo;

end do;

printf("\n Comienza producto matrices Toeplitz");

# Producto de matrices Toeplitz

Vect := Toeplitz[1];

for k from 2 to n do

C := Toeplitz[k];

Q := Vector(k+1);

for i from 1 to k+1 do

Q[i] := normal(add(C[i+1-j]*Vect[j], j = 1..min(k,i)));

end do;

Vect := Vector(k+1);

for i from 1 to k+1 do

Vect[i] := Q[i];

end do;

end do;

polcar := add(Q[i+1]*X^(n-i), i = 0 .. n);

polcar:

else

# Procesador esclavo recibe lista de trabajos

msg := Grid:-Util:-Receive(0);

for i from 1 to nops(msg) do

# Determino el Toeplitz(i) a calcular

posicion := op(i, msg);

# Se calcula el Toeplitz(i)

calculo := calculotoep(posicion);

# msg = [posicion_en_Toeplitz, calculo_realizado]

mensaje := [posicion, calculo];

Grid:-Util:-Send(0, mensaje);

end do;

end if;

end proc;

# testrun

berkopar(x);

46

E.2. Segundo programa paralelo. Reparto de matrices y producto casi paralelo

El primer programa paralelo (ver seccion E.1) no realizaba el producto de matrices Toeplitz en paralelo conel calculo de matrices Toeplitz por parte de los procesadores esclavos. Este segundo programa trata subsanaresa circunstancia. Para conseguirlo, el procesador 0 no se queda ocioso a la espera de recibir la totalidad dematrices Toeplitz. Conforme disponga de un par de matrices procede a calcular su producto. Por tanto se mejorael paralelismo. Sin embargo, como el calculo de las primeras matrices Toeplitz es rapido el procesador maestrodispone de matrices para el producto y los cuatro nucleos estan trabajando. Conforme el ındice de la matrizToeplitz es mayor los procesadores esclavos tardan mas en calcular la matriz Toeplitz y el procesador 0 quedaocioso a la espera de la matriz Toeplitz siguiente para realizar el producto.

Veamos el codigo Maple. Como el caso anterior, nuestro programa se compone de dos ficheros:

1. call berkopar2.mpl : crea el grid y ejecuta el programa berko2.mpl.

2. berkopar2.mpl : codigo paralelo con el reparto y multiplicacion de matrices Toeplitz en paralelo.

El fichero call berkopar2.mpl es similar a call berkopar1.mpl (ver seccion E.1). Por ese motivo no se incluye.

Codigo del segundo programa paralelo. Veamos el contenido del fichero berko2.mpl.



# Producto de matrices Toeplitz (matriz por fila)

prod_toepmxf := proc(Toeplitz2::list, Toeplitz1::list, ncol::integer)

local i, j, C;

C := [seq(i, i=1..ncol+1)];

for i from 1 to ncol+1 do

C[i] := normal(add(Toeplitz2[i+1-j]*Toeplitz1[j], j = 1..min(ncol,i)));

end do;

C:

end proc;

# n




if ExpoA = 0 then

# R x S


else

# n

# R x A



end do;




end do;


T[j] := prxa[j];

end do;

47

end do;

# n

# R x A x S


end if;

end proc:



local i, C;

C := [seq(i, i=1..n+1)];

C[1] := -1;

C[2] := M[n,n];



end do;

C:

end proc:



local thisnode, totalnodes, i, j, k, n, inicio, fin, Toeplitz, C, Q, msg, mensaje,

posicion, calculo, indice, polcar;












Toeplitz := [seq(‘empty‘, i=1..n)];

Toeplitz[1] := [-1, M[1,1]];







end do;

48

indice := 2;

fin := n-(2*(totalnodes-1));

# Recibe de los procesadores algunas matrices Toeplitz





printf("\n Recibido Toeplitz(%a)", posicion);

Toeplitz[posicion] := calculo;

if Toeplitz[indice] <> ‘empty‘ and indice <= fin then

Toeplitz[indice] := prod_toepmxf(op(indice,Toeplitz),op(indice-1,Toeplitz),indice);

printf("\n P0: Calculo Toeplitz(%a) x Toeplitz(%a) = %a", indice,indice-1,Toeplitz[indice]);

Toeplitz[indice-1] := 0;

indice := indice+1;

end if;

end do;

# Se han recibido todas la matrices Toeplitz pero no se han multiplicado todas

for i from 1 to nops(Toeplitz) while op(i,Toeplitz)=0 do end do;

if i < fin then

for j from i+1 to fin do

Toeplitz[j] := prod_toepmxf(op(j,Toeplitz),op(j-1,Toeplitz),j);

printf("\n P0: Calculo Toeplitz(%a) x Toeplitz(%a) = %a", j,j-1,Toeplitz[j]);

Toeplitz[j-1] := 0;

end do;

end if;

else

# Procesador esclavo recibe lista de trabajos



# Determino el Toeplitz(i) a calcular


# Se calcula el Toeplitz(i)


# msg = [posicion_en_Toeplitz, calculo_realizado]



end do;

end if;

end proc;

# testrun

berkopar(x);

49

E.3. Tercer programa paralelo. Reparto de matrices Toeplitz y mejora del pro-ducto de matrices

Si deseamos mayor paralelismo en el producto de matrices del segundo programa paralelo (ver seccion E.2)entonces debemos permitir a los procesadores esclavos el producto de matrices. Esta circunstancia hace mascomplejo el codigo paralelo. El motivo es la dimension de la matriz. Por ejemplo, si la dimension de la matrizde entrada es 5x5 no es posible asignar un par de matrices Toeplitz a los tres procesadores esclavos para suproducto. Otra circunstancia es que siempre empezamos el producto de matrices P2=Toep(Q2)xToep(Q1) dondeToep(Q1) es un vector. Despues seguimos con P3 = Toep(3)xP2 y ası sucesivamente. En todo momento es unproducto matriz por vector. El resultado es otra matriz Toeplitz. Sin embargo si multiplicamos directamenteToep(10)xToep(9) es un producto de matrices y la matriz resultante no es una matriz Toeplitz. En consecuenciamejoramos el paralelismo pero el codigo Maple se hace mas complejo.

El codigo Maple se compone ahora de tres ficheros:

1. call berkopar3.mpl : crea el grid y ejecuta uno u otro programa dependiendo de la dimension de la matrizde entrada.

2. berkopar3seq.mpl : codigo paralelo similar al primer programa paralelo.

3. berkopar3.mpl : codigo paralelo para repartir productos de matrices Toeplitz entre los procesadores es-clavos.

Creacion del grid y ejecucion del programa paralelo. Veamos el contenido del ficherocall berkopar3.mpl.

# Procedimiento para leer un fichero de texto con ordenes Maple

LoadCode := proc(n)

local sep, fname,r,l;


fname := cat(kernelopts(homedir), sep, "maple", sep, "grid", sep, n, ".mpl");

r := "";


while l <> NULL do

r := cat(r, l, "\n");


end do;

FileTools[Text][Close](fname);

r;

end:

# Procedimiento para Launch()

checkAbort := proc()

false;

end proc;



# Se importa la matriz desde el fichero $HOME/maple/grid/matriz con la matriz


# Se comprueba la dimension para codigo secuencial o paralelo

n := LinearAlgebra:-ColumnDimension(M);

if n < 7 then

Code := LoadCode(berkopar3seq);

else

50

Code := LoadCode(berkopar3);

end if;

# Valores previos para iniciar el grid

host := "localhost":

port := 2000:

numNodes := 4:

broadcastAddress := "127.0.0.255":

broadcastPort := 4400:

logFile := "logs/gridlog.txt":

# Se lanza el grid

Grid[Server][StartServer](port,numNodes,broadcastAddress,broadcastPort,logFile);

Grid[Setup](host,port);

# Se comprueba el numero de nodos disponibles en el grid

(TotalNodes, AvailableNodes) := Grid[Status]();

# Se ejecuta el programa Maple correspondiente

r := Grid[Launch](AvailableNodes, Code, printf, checkAbort, ["M"]):

# Se para el grid

Grid[Server][StopServer](port,numNodes);

Para ejecutar este fichero basta teclear:

maple < call_berkopar3.mpl

Codigo si la matriz de entrada es de dimension inferior a siete. Veamos el contenido del ficheroberkopar3seq.mpl.




prod_toepmxf := proc(Toep2::Vector, Toep1::Vector, ncol::integer)

local i, j, C;

C := Vector(ncol+1);


C[i] := normal(add(Toep2[i+1-j]*Toep1[j], j = 1..min(ncol,i)));

end do;

C:

end proc;

# n




if ExpoA = 0 then

# R x S


else

# n

# R x A

51



end do;




end do;


T[j] := prxa[j];

end do;

end do;

# n

# R x A x S


end if;

end proc:



local i, C;

C := Vector(n+1);

C[1] := -1;

C[2] := M[n,n];



end do;

C:

end proc:



local thisnode, totalnodes, i, j, k, n, inicio, fin, Toep, MiToep, Vect, C, Q, msg,

mensaje, posicion, calculo, indice, polcar;













52

Toep := [seq(i, i=1..n)];

# Toep(Q1)

Toep[1] := Vector([-1, M[1,1]]);







end do;

indice := 2;

# Recibe de los procesadores Toep(Q2) a Toep(Qn)





Toep[posicion] := calculo;

printf("\n Recibido Toep(%a)", posicion);

if (Toep[indice] <> indice) then

printf("\n Toep(%a) x Toep(%a)", indice,indice-1);

Toep[indice] := prod_toepmxf(op(indice,Toep),op(indice-1,Toep),indice);

# Se elimina el Toep anterior para ahorrar memoria

Toep[indice-1] := 0;

indice := indice+1;

end if;

end do;

# Resto de Toeps no multiplicadas

for i from indice to n do

printf("\n Toep(%a) x Toep(%a)", i,i-1);

Toep[i] := prod_toepmxf(op(i,Toep),op(i-1,Toep),i);

# Se elimina el Toep anterior para ahorrar memoria

Toep[i-1] := 0;

end do;

Q := op(n,Toep);

polcar := add(Q[i+1]*X^(n-i), i = 0 .. n);

polcar:

else

# Recibe la lista de trabajos



# Determino el Toep(i) a calcular


# Se calcula el Toep(i)


# msg = [posicion_en_Toep, calculo_realizado]



end do;

end if;

end proc;

53

# testrun

berkopar(x);

Codigo si la matriz de entrada es de dimension superior a seis. Veamos el contenido del ficheroberkopar3.mpl.



prod_fxmtoep := proc(MToep::list, VToep::list)

local nc, lr, dimnc, dimvtoep, i, j, C;

nc := op(1,MToep);

lr := op(2,MToep);

dimnc := nops(nc);

dimvtoep := nops(VToep);

C := [seq(i, i=1..dimnc+1)];

for i from 1 to dimnc do

C[i] := normal(simplify(add(nc[i+1-j]*VToep[j], j = 1..min(dimvtoep,i))));

end do;

C[dimnc+1] := normal(simplify(add(lr[i]*VToep[i], i=1..dimvtoep)));

C:

end proc;

# Multiplicacion de dos matrices Toep(n) * Toep(n-1)

muldosToep := proc(A::list,B::list)

local n, i, j, nc, lr;

n := nops(B);

nc := [seq(i, i=1..n)];

lr := [seq(i, i=1..n-1)];

# Columna C de dimension n

for i to n do

nc[i] := normal(add(simplify(A[i-j+1]*B[j]), j = 1 .. i));

end do;

# Ultima fila de dimension n-1

lr[1] := normal(add(A[n+2-j]*B[j], j = 1 .. n));

for i from 2 to n-1 do

lr[i] := normal(simplify(nc[n-i+2]+B[n-i+2]));

end do;

[nc,lr]:

end proc;


prod_toepmxf := proc(Toep2::list, Toep1::list, ncol::integer)

local i, j, C;

C := [seq(i, i=1..ncol+1)];



end do;

C:

54

end proc;




if ExpoA = 0 then

# R x S


else

# n

# R x A



end do;




end do;


T[j] := prxa[j];

end do;

end do;

# n

# R x A x S


end if;

end proc:



local i, C;

C := [seq(i, i=1..n+1)];

C[1] := -1;

C[2] := M[n,n];



end do;

C:

end proc:



local thisnode,totalnodes,i,j,k,n,wl,recv,LToep,msg,pos,toep,toep2,mensaje,indice,fin,val

or,VToep,polcar;






55





wl := [seq(i, i=1..totalnodes-1)];

recv := 0;

for i from 1 to totalnodes-1 do

wl[i] := [seq(n-(j-1)-(2*(i-1)), j=1..2)];

recv := recv+1;

end do;

i := 3;

for j from n-(2*(totalnodes-1)) to 2 by -1 do

wl[i] := [j,op(wl[i])];

i := i-1;

if i=0 then i:=3 end if;

recv := recv+1;

end do;

LToep := [seq(i, i=1..recv+1)];

LToep[1] := [-1, M[1,1]];

for i from 1 to totalnodes-1 do

Grid:-Util:-Send(i,wl[i]);

printf("\n Enviado a P%a lista %a",i,wl[i]);

end do;

printf("\n Recepciones = %a", recv);

indice := 2;

fin := n-(2*(totalnodes-1));

# Recibe de los procesadores

for j from 2 to recv+1 do


pos := op(1,mensaje);

toep := op(2,mensaje);

printf("\n Recibido Toep(%a)", pos);

LToep[pos] := toep;

if LToep[indice] <> indice then

if type(op(1,LToep[indice]),’list’) then

LToep[indice] := prod_fxmtoep(op(indice,LToep),op(indice-1,LToep));

else

LToep[indice] := prod_toepmxf(op(indice,LToep),op(indice-1,LToep),indice);

end if;

printf("\n Calculo Toep(%a)xToep(%a)", indice, indice-1);

LToep[indice-1] := 0;

indice := indice+1;

end if;

end do;

# No multiplicados

for i from 1 to nops(LToep) while op(i,LToep)=0 do end do;

for j from i+1 to recv+1 do

if type(op(1,LToep[j]),’list’) then

LToep[j] := prod_fxmtoep(op(j,LToep),op(j-1,LToep));

else

LToep[j] := prod_toepmxf(op(j,LToep),op(j-1,LToep),j);

56

end if;

printf("\n*Calculo Toep(%a) x Toep(%a)", j, j-1);

LToep[j-1] := 0;

end do;

VToep := op(recv+1,LToep);

polcar := add(VToep[i+1]*x^(n-i), i = 0 .. n);

printf("\n **"); # para mostrar printf anterior

polcar:

else


k := nops(msg);

for i from 1 to k-1 do

pos := op(i,msg);

valor := calculotoep(pos);

if i < k-1 then

mensaje := [pos,valor];

Grid:-Util:-Send(0,mensaje);

else

pos := op(k,msg);

toep2 := calculotoep(pos);

pos := pos-((totalnodes-1)-thisnode);

mensaje := [pos,muldosToep(valor,toep2)];

Grid:-Util:-Send(0,mensaje);

end if;

end do;

end if;

end proc;

# testrun

berkopar(x);

57

E.4. Cuarto programa paralelo. Anillo de nodos

Introduccion. En la seccion E.3 vimos como el codigo Maple se ha vuelto muy complejo. En base al ejemplode un anillo de nodos incluido en Maple Grid Computing Toolbox se ha desarrollado el codigo paralelo paraimplementar el algoritmo paralelo de Berkowitz mediante un anillo de procesadores.

Otro aspecto notable es el cambio en la creacion del grid y el modo de ejecucion del codigo paralelo. Elmotivo han sido los fallos Java:

Error, (in unknown) java.lang.OutOfMemoryError: Java heap spacePese a cambiar el valorJAVAHEAP=400en el fichero ejecutable maple de Maple 11 no se modifica dicho valor cuando se utiliza Maple Grid Computing

Toolbox. Si deseamos el cambio es necesario especificarlo en el script gridserver.sh utilizado para crear el grid.Por tanto, si deseamos evitar el error anterior con ciertas matrices de entrada es necesario cambiar tanto elmodo de creacion del grid como el modo de ejecutar el programa paralelo.

Creacion del grid y ejecucion del programa paralelo. En la seccion 2.2 se detallan las distintas formasde crear el grid. Una de ellas es con el script gridserver.sh. Si deseamos cambiar los valor JAVAHEAP de 400 a800 entonces en el fichero gridserver.sh modificamos la lınea

$java_pgm -Drootdir=$grid_path -jar $grid_jar $pgm_args

por

$java_pgm -Drootdir=$grid_path -jar $grid_jar $pgm_args

Una vez modificado podemos crear el grid con la orden

gridserver.sh -n 4 -b 4400 -a 127.0.0.255 -p 2000 -f logs/gridlog.txt startmultiple

Para ejecutar el fichero Maple anillo.mpl debemos especificar:

launcher.sh -h localhost -n 4 anillo.mpl

Anillo de nodos. Veamos el contenido del fichero anillo.mpl.



# Para evitar el fallo

# Error, (in unknown) java.lang.OutOfMemoryError: Java heap space

kernelopts(limitjvmheap=true);

kernelopts(jvmheaplimit=131072);

# Numero de nodo dentro del grid (de 0 a 3)

me := Grid:-Util:-MyNode();

# Numero de nodos (4 en un Intel QuadCore)

n := Grid:-Util:-NumNodes();

# Necesitamos al menos dos nodos en el grid

if n<2 then error "Grid must have at least 2 nodes" end if;

# Se lee la matriz de un fichero denominado $HOME/maple/grid/matriz




# Se determina la dimension de la matriz

58

cdm := ColumnDimension(M);

# Total de rotaciones en el anillo

fin := floor(cdm/n);

# Procedimiento para calcular el vector columna C que define la matriz Toep


local i, C;

if n = 1 then

C := [-1, M[1,1]];

else

C := [seq(i, i=1..n+1)];

C[1] := -1;

C[2] := M[n,n];



end do;

end if;

C:

end proc:

# n

# Procedimiento que calcula el producto R x A x S



if ExpoA = 0 then

# R x S


else

# n

# R x A



end do;




end do;


T[j] := prxa[j];

end do;

end do;

# n

# R x A x S


end if;

end proc:

# Producto de matrices Toep. Como una de ellas es una matriz columna

# es el producto de una matriz por un vector

prod_toepmxf := proc(Toep2::list, Toep1::list, ncol::integer)

local i, j, C;

59

C := [seq(i, i=1..ncol+1)];



end do;

C:

end proc;

# Calculo del destinatario dependiendo de la posicion en el anillo

if me=n-1 then

destination := 0;

else

destination := me+1;

end if;

# Inicio del programa. El nodo 0 determina sus operaciones iniciales

# y envia al siguiente nodo para el calculo

if me=0 then

t1 := Grid[Time](); # Inicio del cronometro

# Supongamos que tenemos una matriz 14x14 y cuatro procesadores. Como el resto

# de la division 14/4 es 2 entonces el procesador 0 debe calcular dos Toep

# adicionales para resolver en tres pasos las 14 matrices Toep en un anillo de

# cuatro procesadores. El reparto de calculos entre procesadores es:

#

# Proc3 Proc2 Proc1 Proc0

#

# Toep(Q6) Toep(Q5) Toep(Q4) Toep(Q3),Toep(Q2),Toep(Q1)

# Toep(Q10) Toep(Q9) Toep(Q8) Toep(Q7)

# Toep(Q14) Toep(Q13) Toep(Q12) Toep(Q11)

calculo := cdm mod n;

msg := calculotoep(1);

for i from 1 to calculo do

toep := calculotoep(i+1);

printf("\n Calculating Toep(%a)xToep(%a)", i+1,i);

msg := prod_toepmxf(toep,msg,i+1);

end do;

Grid:-Util:-Send(destination, msg);

end if;

# Siguiendo con el ejemplo, en este momento Proc0 ha calculado

#

# P2 = Toep(Q2)xP1 siendo P1 = Toep(Q1)

# P3 = Toep(3)xP2

#

# finalmente el procesador 0 envia la matriz P3 al procesor 1

#

# Cada nodo calcula en paralelo su Toep correspondiente, lo multiplica

# por el Pi recibido, obteniendo un nuevo Pi+1 que se envia al siguiente

# nodo en el grid. El reparto de multiplicaciones es:

#

# Proc0 Proc1 Proc2 Proc3

#

# P3 --> P4=Toep(Q4)xP3 --> P5=Toep(Q5)xP4 --> P6=Toep(Q6)xP5 -->

60

# P7=Toep(Q7)xP6 --> P8=Toep(Q8)xP7 --> P9=Toep(Q9)xP8 --> P10=Toep(Q10)xP9 -->

# P11=Toep(Q11)xP10 --> P12=Toep(Q12)xP11 --> P13=Toep(Q13)xP12 --> P14=Toep(Q14)xP13 -->

for i from 1 to fin do

r := Grid:-Util:-Receive();

indice := nops(r);

if indice <= cdm then

printf("\n Calculating Toep(%a)", indice);

toep := calculotoep(indice);

msg := prod_toepmxf(toep,r,indice);

printf("\n Calculating Toep(%a)xToep(%a)", indice+1,indice);

Grid:-Util:-Send(destination, msg);

unassign(’toep’,’msg’); # Libera memoria

end if;

end do;

# El nodo 0 ha recibido la matriz solucion P14

if me=0 then

t2 := Grid[Time]();

printf("\n Total time : %a \n", t2-t1); # Tiempo invertido en el grid

polcar := add(r[i+1]*x^(cdm-i), i = 0..cdm):

else

polcar := NULL;

end if;

polcar:

61

Referencias

[1] J. Abdeljaoued, Algorithmes rapides pour le calcul du polynome caracteristique. Tesis Doctoral, Univer-site de Franche-Comte, Francia (1997).

[2] J. Abdeljaoued, Berkowitz Algorithm, Maple and computing the characteristic polynomial in an arbitrarycommutative ring. Computer Algebra MapleTech 4(3), Birkhauser (1997).

[3] J. Abdeljaoued, G.M. Diaz-Toca and L. Gonzalez-Vega, Minors of Bezout matrices, Subresultants and pa-rameterization of the degree of the polynomial greatest common divisor. International Journal of ComputerMathematics 81, 10, 1223–1238 (2004).

[4] J. Abdeljaoued, G.M. Diaz-Toca and L. Gonzalez-Vega, Bezout matrices, Subresultant Polynomials andParameters. Preprint (2008).

[5] J. Abdeljaoued, H. Lombardi, Methodes matricielles. Introduction a la complexite algebrique. Springer(2004).

[6] S. Basu, R. Pollack and M.–F. Roy, Algorithms in Real Algebraic Geometry. Algorithms and Computationsin Mathematics 10, Springer–Verlag (2003).

[7] S. J. Berkowitz: On computing the determinant in small parallel time using a small number of processors.Information Processing Letters 18, 147–150 (1984).

[8] G. E. Collins, Subresultants and reduced polynomial remainder sequences. Journal of Association for Com-puting Machinery 14, 128–142 (1967).

[9] G. Diaz-Toca, L. Gonzalez-Vega and H. Lombardi, Generalizing Cramer’s rule: Solving uniformly linearsystems of equations. SIAM J. Matrix Anal. Appl. 27 (3), 621–637 (2005).

[10] G. Diaz–Toca and L. Gonzalez–Vega: Various New Expressions for Subresultants and Their Appplications.Applicable Algebra in Engineering, Communication and Computing 15, 233–266 (2004).

[11] Dana Petcu, Working with multiple Maple kernels connected by Distributed Maple or PVMaple. RISC-LinzTechnical Report Series No. 01-18 (2001).

[12] D. Dubu, Interconnecting Computer Algebra Systems within the Grid. Master Thesis (2004).

[13] I.Z. Emiris and V.Y. Pan, Improved algorithms for computing determinants and resultants.Journal of Com-plexity 21, 43–71 (2005).

[14] L. Gonzalez-Vega, H. Lombardi, T. Recio and M.-F. Roy, Specialisation de la suite de Sturm et sous-resultants (I). Informatique Theorique et Applications 24(6), 561–588 (1990).

[15] M. Kerber, Division–free Computation of Subresultants Using Bezout Matrices. To appear in IJCM (2008).

[16] Y.B. Li, A new approach for constructing subresultants. Applied Mathematics and Computation 183 471–476 (2006).

[17] H. Lombardi, M.–F. Roy and M. Safey, New structure theorem for subresultants. Journal of SymbolicComputation 29, 663–689 (2000).

[18] R. Loos: Generalized polynomial remainder sequences. Computer Algebra, Computing Suplementum 4,115–138, Springer–Verlag (1982).

[19] A. Marco and J-J. Martinez, Parallel computation of determinants of matrices with polynomial entries.Journal of Symbolic Computation 37, 749–760 (2004).

[20] D. Petcu, Working with multiple Maple kernels connected by Distributed Maple or PVMaple.

[21] J. von zur Gathen and T. Lucking, Subresultants revisited. Theoretical Computer Science 297, 199–239(2003).

62

[22] D. Wang: Subresultants with the Bezout Matrix. Computer Mathematics. Proceedings of the Fourth AsianSymposium on Computer Mathematics (ASCM 2000), 19–28, World Scientific (2000).

63

Paralelizaci´on en sistemas multicore con Maple · acceder desde programas escritos en lenguaje C,...

Documents

Transcript of Paralelizaci´on en sistemas multicore con Maple · acceder desde programas escritos en lenguaje C,...