FR2610429A1 - Procede d'analyse a elements finis utilisant une matrice maximisee en largeur de bande - Google Patents

Abstract

LA PRESENTE INVENTION CONCERNE UN PROCEDE D'ANALYSE A ELEMENTS FINIS DANS LEQUEL UNE AUGMENTATION DE LA VITESSE DE RESOLUTION DANS DES PROCESSUS DE RETROSUBSTITUTION ET D'ELIMINATION DIRECTE EST OBTENUE EN SUPPRIMANT LES DEPENDANCES ENTRE LES DONNEES QUI OBLIGENT A DES CALCULS SEQUENTIELS. CECI EST REALISE PAR UN PROCESSUS DIT DE MAXIMISATION DE LARGEUR DE BANDE QUI TRANSFORME UNE MATRICE DE SYSTEME EN UNE MATRICE TRIANGULAIRE DECOMPOSEE AYANT DES DONNEES NON NULLES SUR SA DIAGONALE ET DANS SON COIN SUPERIEUR DROIT OU INFERIEUR GAUCHE, LES DONNEES NON NULLES ETANT SEPAREES PAR UNE BANDE PRESELECTIONNEE K DE DONNEES TOUTES NULLES. UNE TELLE MATRICE MAXIMISEE EN LARGEUR DE BANDE EST CREEE EN RENUMEROTANT LES NOEUDS DU RESEAU A ELEMENTS FINIS (OU PLUS PRECISEMENT LES NUMEROS DE RANGEES ET DE COLONNES DES INCONNUES ASSOCIEES AUX NOEUDS) DE SORTE QUE LES NUMEROS ASSIGNES AUX NOEUDS ADJACENTS DIFFERENT D'UNE QUANTITE SUFFISANTE POUR CREER LA BANDE SUSMENTIONNEE DE DONNEES TOUTES NULLES.

Description

PROCEDE D'ANALYSE A ELEMENTS FINIS UTILISANT UNE MATRICE

MAXIMISEE EN LARGEUR DE BANDE

De nombreux systèmes physiques peuvent être décrits mathématiquement en terme de systèmes d'équations linéaires qui

sont eux-mêmes résolus par des procédés de manipulation matri-

cielle. L'analyse à éléments finis est un procédé qui concerne

la description de divers systèmes physiques en termes de

systèmes d'équations et le développement de méthodologies pour la solution de tels systèmes. Le terme "système physique" est

interprété ici comme se référant à une structure, des dispo-

sitifs, un appareil, ou des corps de matière (solide liquide ou gaz) ou simplement une région de l'espace dans laquelle un

phénomène particulier physique, chimique ou autre prend place.

L'analyse à éléments finis était à l'origine un procédé d'analyse structurelle mais est aujourd'hui couramment utilisée dans la conception de moteurs, générateurs, systèmes d'imagerie à résonance magnétique, systèmes d'allumage de moteurs d'avions, rupteurs et transformateurs de circuit, pour en

nommer quelques uns; ces techniques sont utilisées pour ana-

lyser des contraintes, températures, structures moléculaires, champs électro-magnétiques, courants, forces physiques, etc dans tout type de systèmes physiques. C'est devenu une partie classique d'un cycle de conception pour de nombreux produits qui ne sont pas facile à analyser par d'autres procédés. La présente invention s'applique en particulier à l'analyse et à

la conception de tels produits.

Des systèmes d'équations linéaires qui doivent être résolus par des techniques d'analyse à éléments finis sont très souvent importants et pour cette raison difficile à résoudre par le calcul. Par exemple, un système d'équations pour une analyse par éléments finis bidimensionnelle importante mais non extraordinaire peut comprendre 25 000 inconnues. Quand de telles équations sont basées sur un réseau d'éléments finis présentant des contributions pour une majorité de noeuds, aucun choix n'existe sauf d'utiliser une force de calcul importante pour arriver à une solution. Toutefois, dans certains cas, ces

équations sont à la fois importantes et creuses et en consé-

quence permettent de prétraiter ou de transformer les équations d'une façon qui exige moins de calcul pour les résoudre. Le

terme "creux" est utilisé ici pour se référer à la caractéris-

tique selon laquelle un très petit pourcentage des éléments d'une matrice a des valeurs non nulles. Quand un caractère creux extrême existe dans un très grand système, il existe

plusieurs techniques qui peuvent être utilisées pour transfor-

mer le système d'équations en un système qui est plus faci-

lement traité du point de vue du calcul. Toutefois, en dépit de ces transformations, les techniques de calcul classiques peuvent être ou bien malaisées ou bien très inefficaces selon

la dimension et d'autres caractéristiques des équations matri-

cielles résultantes.

Comme on peut le comprendre à partir de l'exposé ci-

dessus, le champ de l'analyse à éléments finis s'est développé dans une grande mesure en raison de la disponibilité de machines de calcul plus importantes et plus puissantes pouvant

être utilisées pour résoudre de tels systèmes. Il existe main-

tenant divers systèmes d'ordinateur spécialisés conçus pour réaliser des calculs d'applications particulières qui sont particulièrement coûteux à réaliser sur des ordinateurs à usage

général. Un tel système est basé sur le concept d'une architec-

ture systolique et fournit une méthodologie générale pour

mapper des calculs de haut niveau dans des structures maté-

rielles. Dans un système systolique, les données s'écoulent à partir de la mémoire de l'ordinateur de façon rythmée, en passant par de nombreux éléments de traitement en chaîne ou en pipe-line avant de retourner vers la mémoire, permettant ainsi

des calculs multiples pour chaque accès en mémoire et entraî-

nant une grande augmentation de la vitesse d'exécution de problèmes de calculs intensifs sans augmentation associée des exigences d'entrée/sortie. Quelques procédés pour mettre en

oeuvre des architectures systoliques pour traiter des opé-

rations matricielles sont traités dans un article de H.T. Kung et C.E. Leiserson ayant pour titre "Systolic Arrays (for VLSI)", Sparse Matrix Proc. 1978, Society for Industrial and Applied Mathematics, 1979, page 256-282. Une autre analyse de ce problème et une solution suggérée sont traitées dans l'article de IEEE Transactions on Computers, Vol. C-32, N 3, mars 1983 ayant pour titre "An Efficient Parallel Algorithm for

the Solution of Large Sparse Linear Matrix Equations".

Il est également décrit dans l'art antérieur des techniques et des agencements pour réaliser diverses opérations matricielles qui sont communes au procédé d'analyse à éléments finis, entre autres, et qui utilisent un processeur parallèle

pour accroître la vitesse de calcul. Ceci est réalisé essen-

tiellement en agençant la solution de sorte qu'elle est réali-

sée en tant qu'étapes multiples, itératives et sensiblement identiques en parallèle par une série de processeurs en parallèle. Par exemple, on connaît un procédé pour mémoriser

une matrice creuse importante dans une architecture à proces-

seurs multiples d'une façon qui rend les éléments de la matrice

facilement disponibles pour diverses opérations matricielles.

Une telle opération, utilisée pour résoudre des systèmes d'équations linéaires, consiste en une rétro-substitution mais

d'autres apparaîtraient d'elles-mêmes à l'homme de l'art.

2610 4 2 9

La solution aux équations linéaires produites par suite des techniques selon la présente invention peut être mise

en oeuvre par un grand nombre d'architectures à multipro-

cesseurs parallèles comme on y a fait allusion ci-dessus.

Toutefois, à titre d'exemple, la solution des équations fournie par suite des procédés décrits ici sera spécifiquement décrite en relation avec l'architecture à réseau systolique telle que généralement décrite dans l'article de E. Arnould et al. ayant pour titre "Systolic Array Computer" présenté lors de "l'IEEE International Conference on Acoustics, Speech, and Signal Processing", 26-29 mars 1985, pp. 232-235. Le procédé selon l'invention peut être mis en oeuvre dans l'appareil décrit dans le brevet des Etats-Unis d'Amérique 4 493 048 ayant pour titre "Systolic Array Apparatuses for Matrix Computations" délivré au

nom de H.T. Kung et al., dont la description sera considérée

ici comme connue.

Un inconvénient principal des procédés et appareils de l'art antérieur pour résoudre des systèmes d'équations linéaires par les procédés de rétro-substitution ou d'élimination directe résulte des dépendances des données qui obligent à des calculs séquentiels plutôt que simultanés ou

parallèles pour résoudre les inconnues.

Un premier objet de l'invention est en conséquence de prévoir un appareil et des procédés mis en oeuvre par une machine perfectionnée pour analyser des systèmes physiques

représentés par des systèmes d'équations linéaires.

Un autre objet est d'actionner un multiprocesseur de

façon à permettre une analyse plus efficace de produits manu-

facturés et de systèmes physiques représentés par des systèmes

d'équations linéaires.

Un autre objet de l'invention est de prévoir un procédé et un appareil pour augmenter la vitesse et l'efficacité de solution des systèmes d'équations linéaires sur des processeurs parallèles, en particulier sur des processeurs

pipe-line parallèles.

Un autre objet est de prévoir un procédé pour traiter des matrices creuses de grande dimension sur un système d'ordinateur à réseau systolique qui supporte un degré de

concurrence élevé, tout en utilisant seulement des communi-

cations et des commandes régulières et simples pour permettre

une mise en oeuvre efficace.

Un autre objet de l'invention est de prévoir une méthodologie pour améliorer la vitesse de calcul de résolution de systèmes d'équations linéaires en utilisant des techniques de rétro-substitution et d'élimination directe qui agissent sur une matrice décomposée triangulaire qui a été "maximisée en largeur de bande" c'est-à-dire transformée pour fournir une bande d'éléments nuls de largeur minimale présélectionnée pour séparer les éléments non nuls de la diagonale principale d'éléments non nuls dans les coins supérieur droit ou inférieur

gauche de la matrice.

Un autre objet de l'invention est de prévoir une technique pour renuméroter les noeuds d'une maille d'éléments finis en mettant en oeuvre le procédé d'analyse à éléments finis de sorte qu'une matrice décomposée caractérisée par la

bande susmentionnée est produite.

Un autre objet de l'invention réside dans la prévision d'un processeur parallèle et d'une méthodologie associée pour faire fonctionner le processeur parallèle de façon à résoudre un système d'équations linéaires caractérisés par une matrice de système du type susmentionné de façon

nouvelle et plus efficace.

Ces objets et d'autres de l'invention sont réalisés au moyen de procédés et de techniques pour mettre en oeuvre le procédé d'analyse à éléments finis pour analyser un système physique qui comprend l'étape consistant à produire une matrice triangulaire décomposée pour le système sous une forme ayant une bande d'éléments ou données de valeurs nulles, de largeur

de bande minimale prédéterminée, adjacente à sa diagonale prin-

cipale et à résoudre le système d'équations associé à une telle

2 6 1 0 4 2 9

matrice sur un processeur parallèle qui agit sur les éléments de données de la matrice mémorisée dans les mémoires d'une pluralité de processeurs. En raison de la présence de la bande de données de valeurs nulles dans la matrice décomposée comme cela a indiqué ci-dessus, des calculs de composantes de vecteurs inconnus du système peuvent être effectués de façon nouvelle, très concurrentielle et efficace, le calcul des valeurs pour une pluralité de composantes d'un vecteur inconnu prenant place par rétro-substitution, indépendamment de la solution antérieure des autres valeurs des composantes du vecteur inconnu. Ayant supprimé le besoin de retarder le processus de rétro-substitution pour une composante vectorielle

donnée jusqu'à ce que des calculs antérieurs d'autres compo-

santes du vecteur aient été achevés, la vitesse de résolution

d'ensemble du système est notablement accrue.

Ces objets, caractéristiques et avantages ainsi que d'autres de la présente invention seront exposés plus en détail

par la description suivante de modes de réalisation parti-

culiers faite en relation avec les figures jointes parmi lesquelles: la figure 1 représente une architecture systolique fondamentale qui est utilisée selon un procédé de l'invention pour traiter des matrices creuses de grande dimension; la figure 2 est un organigramme illustrant le procédé pour mémoriser une matrice creuse de grande dimension dans la mémoire d'un multiprocesseur parallèle; la figure 3 illustre un procédé de mappage d'une matrice particulière, constituant un exemple de matrice creuse de grande dimension à largeur de bande maximisée dans le réseau multiprocesseur parallèle selon une mise en oeuvre de l'invention; la figure 4 illustre un procédé pour mettre en oeuvre la procédure de rétrosubstitution pour résoudre un système triangularisé d'équations linéaires, ce système étant

partiellement caractérisé par la matrice à décomposition spéci-

2610 4 2 9

figue maximisée en largeur de bande mémorisée selon la figure 3 dans un système multiprocesseur parallèle; les figures 4A-4M représentent le débit des données dans le multiprocesseur de la figure 4 pendant des cycles successifs de la machine pour mettre en oeuvre le processus de

rétro-substitution de l'exemple de matrice mémorisée à décom-

position maximisée en largeur de bande de la figure 4; la figure 5A représente la structure à caractère creux d'une matrice à décomposition minimisée en largeur de bande selon l'art antérieur; et la figure 5B représente la structure à caractère creux d'une matrice à décomposition maximisée en largeur de

bande selon l'invention.

Comme on l'a noté ci-dessus, le procédé selon l'invention concerne des perfectionnements au procédé d'analyse

à éléments finis utilisé pour analyser des systèmes physiques.

Dans de tels systèmes, une variable de champ (qui peut être la pression, la température, le déplacement, la contrainte ou une autre quantité) possède une infinité de valeurs composantes car c'est une fonction de chaque point géométrique dans le corps ou zone de solution en cours d'investigation. Dans la première étape du procédé, le problème est rendu discret en un problème impliquant un nombre fini (quoique important) d'inconnues en divisant la zone de solution en éléments et en exprimant la variable inconnue en terme de fonctions d'approximation dans chaque élément. Ainsi, le problème complexe se réduit à considérer une succession de problèmes très simplifiés. Les fonctions d'approximation sont définies en termes de valeurs inconnues des variables aux points spécifiés appelés noeuds. Ainsi, des équations matricielles sont écrites

pour exprimer des propriétés des éléments individuels.

En général, chaque noeud est caractérisé par une ou

plusieurs équations linéaires qui associent une variable résul-

tante ou liée R à une variable de champ Y aux noeuds donnés et par une constante de "raideur" K. Plus particulièrement, pour une zone de solution en cours d'investigation comprenant un grand nombre de tels noeuds, la variable vectorielle résultante ou liée [R] ayant des composantes [Ri, R2... RN] (qui sont connues aux noeuds limites et ont une valeur nulle), s'exprime en termes d'une variable de champ vectoriel [Y] ayant des composantes [Y1, Y2...YN] multipliées par la matrice de raideur [K] constituée des constantes ou coefficients des variables de champ aux noeuds, par l'équation matricielle

[K] [Y] = [R]

Par exemple, dans l'analyse d'un système de ressorts linéaires, les composantes de variables résultantes représentées par [RI, R2...] peuvent être des valeurs de forces appliquées au système à des noeuds choisis; des composantes de variables de champ [Y1, Y2...] peuvent être les valeurs de déplacement aux noeuds, et les constantes peuvent être des valeurs de raideurs de ressorts des éléments de ressorts étudiés qui associent la force au déplacement au niveau des noeuds. Les constantes

forment une matrice de coefficients ou de raideurs [K].

L'étape suivante dans le procédé d'analyse à élé-

ments finis consiste à résoudre l'équation matricielle ci-

dessus. Les techniques de machine pour manipuler l'équation

matricielle pour arriver automatiquement à une solution rapi-

dement et efficacement pour des systèmes ou équations repré-

sentant des systèmes physiques complexes analysés par le procédé d'analyse fini constituent l'objet principal de cette demande de brevet. On pourra trouver une explication plus détaillée des principes du procédé d'analyse à éléments finis dans les textes ayant pour titre "The Finite Element Method for Engineers" de K. H. Huebner et " Finite Element Procedure

Procedure in Engineering Analysis" de Klaus-Jurgen Bathe.

De nombreuses techniques ou procédés existent pour résoudre le système d'équations susmentionné. De telles techniques sont souvent extrêmement complexes dans leur nature

et impliquent l'utilisation de décompositions et de transfor-

mations des équations du système qui sont obtenues à partir des

équations initiales du système et qui en sont ou bien des équi-

valents ou bien des approximations.

Ainsi, par exemple, les équations initiales du système

[K] [Y] = [R]

peuvent être transformées sous la forme

[A] [X] = [Q]

o [A] est une matrice décomposée connue d'ordre N par N et Q est un vecteur à N dimensions transformé connu et X est un

vecteur à N dimensions transformé inconnu.

Le but de ces transformations et de la décomposition du système initial d'équations est de permettre l'application de diverses techniques matricielles connues sur le système

transformé pour arriver à une solution finale.

La rétro-substitution et l'élimination directe sont deux telles techniques numériques qui sont utilisées pour résoudre des équations d'un système obtenu à partir de

l'utilisation du procédé d'analyse à éléments finis.

La technique de rétro-substitution est utilisée quand la matrice décomposée ou transformée [A] est sous forme triangulaire (supérieure ou inférieure), tous les éléments diagonaux de la matrice étant non nuls. Elle peut être utilisée une fois pour arriver à un résultat final définitif ou elle

peut être utilisée itérativement en tant que partie de tech-

niques d'approximation plus complexes pour arriver à une solu-

tion finale. Des techniques pour transformer un système linéaire de forme générale en un système de forme triangulaire sont bien connues. Un tel système triangulaire a la forme: AllXl+A12X2+...A1,N-!XN-1+A1NXN = Q1 A22X2+...A2,N-lXN-1+A2NXN = Q2 AN-1,N-1XN-i+AN-1,NXN = QN-1

ANNXN = QN

Pour résoudre ce système d'équations par rétro-

substitution il faut d'abord noter que la dernière équation peut être résolue immédiatement pour XN puisque ANN (un élément diagonal de la matrice de raideur) et QN (une composante de la variable vectorielle liée) sont connus. Connaissant XN, la pénultième équation peut être résolue puisqu'une seule inconnue

existe, à savoir XN-1.

Plus particulièrement,

XN-1 = [QN-1 - AN-1,NXN] / AN-1,N-1,

XN et XN-1 étant connus, l'antépénultième d'équation peut être résolue puisqu'elle contient une seule vraie inconnue à savoir XN-2. Ainsi, de façon générale pour i = N, N-l...1 N Xi = [Qi - Z XjAi, j)] /Aj,i j = i+l Il faut noter que quand i = N, la sommation se lit qui j j= N+i s'interprète comme une somme sur un nombre de terme nul et

donne par convention la valeur 0.

Plusieurs caractéristiques du processus de rétro-

substitution présentent la possibilité de traitements très

simultanés qui seront mis oeuvre ici par le système à multi-

traitement. La première est que XN, une fois calculé, doit être multiplié par chaque élément de la ième colonne de la matrice pendant le processus de résolution des composantes de X. Ainsi, en calculant XN-1, XN est multiplié par AN-1,N (le coefficient de la N-lième rangée et de la Nième colonne. De même, pour

calculer XN-2...X1, XN est multiplié respectivement par AN-

2,N...A1,N (les coefficients restants de la Nième colonne). De façon similaire, le terme XN-1, une fois calculé, est multiplié

par chaque coefficient de la N-lième colonne.

il Un appareil et un procédé pour mettre en oeuvre efficacement le processus indiqué ci-dessus dans un processeur

parallèle ont déjà été décrits. Toutefois, dans l'art anté-

rieur, le procédé et l'appareil décrits sont limités de façon inhérente du fait que les dépendances des données dans les opérations requises limitent le parallélisme disponible pour résoudre les équations à éléments finis et ralentissent donc la

sortance de l'architecture à multiprocesseur pour cette appli-

cation. Plus spécifiquement, le procédé d'obtention de XN-1,XN-

2...--, etc, est de façon inhérente séquentiel puisque, à partir de l'équation ci-dessus, Xi ne peut pas calculer tant que les Xi+I, Xi+2..., XN précédents n'ont pas été obtenus, leurs multiplications en colonnes associées réalisées et les

résultats partiels accumulés.

Quand le système physique analysé est caractérisé par une équation de système de grande dimension et creuse (ayant un grand pourcentage de coefficients nuls) l'utilisation de décompositions et de transformations regroupées est devenu

tout à fait courante.

Le regroupement se réfère au confinement des termes non nuls dans des parties spécifiques ou bandes de la matrice de décomposition et peut mieux se comprendre en relation avec les figures 5A et 5B. En figures 5A et 5B, les structures de coefficients ou empreintes de deux matrices sont représentées avec des lignes diagonales pleines représentant des éléments non nuls, des zones hachurées représentant des zones contenant des éléments de données nuls et non nuls, et des zones blanches ou noires représentant seulement des éléments nuls. On voit ainsi qu'à la figure 5A les éléments non nuls de la matrice sont confinés (en même temps qu'un certain nombre d'éléments

nuls) vers une bande de largeur b de chaque côté de la diago-

nale. Seuls des éléments nuls sont disposés dans les coins supérieur droit et inférieur gauche de la matrice. Ce type de regroupement est appelé "minimisation de largeur de bande" puisque les termes non nuls sont confinés à de petites bandes

2 6 10 4 2 9

adjacentes aux diagonales. Comme cela est clair pour l'homme de l'art, une telle matrice minimisée en largeur de bande réduit le nombre de calculs requis pour réaliser le processus de rétro-substitution puisqu'un nombre connu de calculs (ceux impliquant tous les éléments dans les zones noires) entraînera inévitablement un résultat nul. Toutefois, la minimisation de largeur de bande n'élimine pas le problème de la dépendance des données notées ci-dessus puisque le calcul des composantes inconnues de la variable de champ doit néanmoins se dérouler

séquentiellement.

La formation d'un système d'équations minimisé en largeur de bande et de sa matrice de décomposition minimisée

correspondante est réalisée pendant le procédé d'analyse à élé-

ments finis par suite de techniques connues pour numéroter ou renuméroter les noeuds dans la maille du système à éléments finis de sorte que les différences entre les numéros assignés

aux noeuds adjacents dans la maille à éléments finis sont mini-

misées. Les techniques (incluant des programmes appropriés) pour la numérotation des noeuds pour atteindre une minimisation de largeur de bande comme on l'utilise couramment sont exposées en détail dans les publications suivantes: R. Rosen, "Matrix Bandwidth Minimization", Proceedings National Conference A.C.M., 1968, pp. 585-595. E.H. Cuthill et J.M. McKee,

"Reducing the Bandwidth of Sparse Symmetric Matrices", Procee-

dings National Conference A.C.M., 1969, pp. 151-172.

Pour supprimer les dépendances de données susnotées dans la résolution des équations à éléments finis et pour permettre une solution plus efficace de ces équations sur des processeurs parallèles, il est décrit ici une technique de "maximisation de largeur de bande", c'est-à-dire l'utilisation d'une matrice décomposée ayant une structure de coefficients telle que représentée en figure 5B. Une telle matrice (et son système correspondant d'équations) est caractérisée par la présence de données uniquement non nulles sur sa diagonale

principale; une bande d'éléments de données tous nuls adja-

cents à la diagonale; et toutes les autres données (incluant toutes les autres données non nulles) rassemblées dans les coins supérieur droit et inférieur gauche de la matrice. La production d'une matrice ayant une empreinte maximisée en largeur de bande est facilement obtenue en numérotant ou en renumérotant les noeuds de la maille d'éléments finis de sorte que des noeuds adjacents diffèrent d'au moins une certaine valeur minimale prédéterminée K. Ainsi, seuls des éléments de données nuls sont présents dans la bande de largeur K de part

et d'autre de la diagonale de la matrice.

Quand un système d'équations caractérisé par une matrice décomposée de la forme représentée en figure 5B est résolu en utilisant le procédé de rétro-substitution sur un

système à processeurs en parallèle fonctionnant selon la métho-

dologie décrite plus en détail ci-après, une grande efficacité est obtenue, en comparaison des procédés et de l'appareil de l'art antérieur. Cette efficacité résulte de la modification

des dépendances des données quand on effectue la rétro-substi-

tution ou l'élimination directe sur une architecture parallèle

selon les procédés de fonctionnement décrits.

La rétro-substitution pour la solution des K premières composantes de la variable de champ inconnue devient maintenant

XN = QN/ANN

XN-1 = QN-1/AN-1,N-1,

et de même, XN-2 = QN-2/AN-2,N-2, etc. Ceci est lié au fait que AN-1,N, AN-2,N, et AN-2,N-1 sont tous des zéros en raison de la forme généralisée maximisée

en largeur de bande de la matrice de solution, telle que repré-

sentée en figure 5B. De façon générale, la rétro-substitution devient, quand une matrice à largeur de bande maximisée est utilisée: pour i = N, N-i,...., l, faire Xi = (Qi - Xj Aij)/A ii j=i+k+1

Ainsi, tout Xi donné peut être calculé indépen-

damment des K Xj précédents tels que j - i 2 K. K est choisi de sorte qu'il correspond à la quantité de temps maximale que cela prend de collecter toutes les données dont dépend un Xi donné, comme cela sera représenté plus en détail en relation avec le fonctionnement du multiprocesseur de la figure 4, en mettant en oeuvre le processus de rétro-substitution en utilisant

l'exemple ci-dessous de matrice maximisée en largeur de bande.

Un bref rappel de l'algèbre matricielle en temps qu'outil de calcul pour résoudre des systèmes d'équations linéaires peut être utile comme arrièreplan à la compréhension

de l'invention.

Une matrice, dans sa forme la plus usuelle est un réseau rectangulaire de quantités scalaires comprenant "m"

rangées agencées orthogonalement par rapport à "n" colonnes.

L'ordre d'une matrice est défini par son nombre de rangées

multiplié par son nombre de colonnes et ainsi la matrice repré-

sentée ci-dessous est appelée matrice d'ordre "m x n". Le mode habituel compact de représentation d'une matrice est illustré ci-dessous. All A12.. . Aln A21 A22... A2n [A] =Amn = Aml Am2... Amn Une rangée d'une matrice est une ligne horizontale ou un réseau unidimensionnel de quantités, alors qu'une colonne est une ligne verticale ou un réseau unidimensionnel de quantités. Les quantités Ail, A12, etc. sont dites être les éléments de la matrice [A]. Une matrice dans laquelle n = 1 est une matrice colonne ou un vecteur colonne, et une matrice o m = 1 est une matrice rangée ou un vecteur rangée. Encore plus généralement, un vecteur est défini comme un ensemble ordonné d'ordre n de valeurs. Une ligne d'une matrice est ou bien une

rangée ou bien une colonne.

Une matrice carrée est une matrice dans laquelle le

nombre de rangées n est égal au nombre de colonnes m. La diago-

nale d'une matrice carrée comprend les éléments Ail, A22, A33, ANN. Une matrice triangulaire est une matrice qui contient tous ses éléments nonnuls dans sa diagonale supérieure ou inférieure. Une matrice triangulaire supérieure a tous ses éléments non nuls dans et au-dessus de sa diagonale principale; une matrice triangulaire inférieure a tous ses

éléments nuls dans et en-dessous de sa diagonale principale.

Dans l'opération [A] x [X] = [Q] o [A] est d'ordre m x n et [X] est d'ordre n x p, [Q] est d'ordre m x p. Ainsi, quand [X] est un vecteur colonne, [Q] sera de même un vecteur colonne, ayant le même nc.mbre, m, de rangées que le nombre de rangées de [A]. Il faut ncter que la multiplication matricielle

est définie seulement quand la matrice [X] dans l'exemple ci-

dessus a le même nombre de rangées que le nombre de colonnes de la matrice [A]. Dans la multiplication ci-dessus le produit [Q] est obtenu en multipliant les colonnes de [X] par les rangées de [A] de sorte que, en général, n Qij= Y Aik x Xkj k =1 Selon ce qui précède, la multiplication d'une matrice [A] par un vecteur multiplication [X] procède normalement de la façon ci-dessous: Ail A12 A13 Xl Q1 [Q] = A21 A22 A23 x X2 = Q2

A31 A32 A33 X3 Q3

AllXi + A12X2 + A13X3

= A21X1 + A22X2 + A23X3

A31X1 + A32X2 + A33X3

En considérant le temps de traitement requis, il serait plus efficace d'effectuer toutes les multiplications impliquant chaque élément vectoriel spécifique en même temps d'o il résulte moins d'accès intermédiaires en mémoire. Ainsi, il serait plus efficace d'un point de vue entrée/sortie de réaliser la multiplication ci-dessus en réalisant toutes les

opérations AllXl, A21X1 et A31XI simultanément ou en parallèle.

Toutefois, ce processus entraîne une dispersion des divers éléments constituants de [Q] qui devraient ensuite être

rassemblés ou recombinés.

Une seconde caractéristique du processus de multi-

plication ci-dessus est que chaque élément Q1, Q2 et Q3 du produit vectoriel ne résulte pas d'une accumulation ou d'une sommation basée sur la rangée d'origine (ou indice de rangée) de l'élément matriciel. En se référant spécifiquement à

l'exemple ci-dessus, on notera que Q1 est la somme des résul-

tats partiels obtenus en multipliant des éléments différents de

la rangée 1 de la matrice par des éléments associés ou corres-

pondant du vecteur [X].

Ainsi, une observation essentielle qui est utilisée dans la méthodologie de la présente invention est que, pour augmenter la simultanéité de fonctionnement et minimiser les exigences d'entrée/sortie de la multiplication ci-dessus, il faut procéder temporellement de la façon suivante: lère opération ---> Xi * [Ail A21 A31] 2ème opération ---> X2 * [A12 A22 A32] 3ème opération ---> X3 * [A13 A23 A33] ce qui précède étant suivi par une sommation de produits partiels basés sur l'indice de rangée des éléments matriciels qui étaient impliqués dans ces produits partiels de la façon suivante: Qi All X1 + A2 X2 + A13 X3

Q2 = A21 X1 + A22 X2 + A23 X3

Q3 A31 X1 + A32 X2 + A33 X3

L'observation plus générale qui découle de ce qui précède est de nombreuses opérations matricielles peuvent être réalisées d'une façon très parallèle en réalisant la première opération commune sur une série d'éléments différents provenant d'une ligne d'une matrice suivie d'une seconde opération basée

sur l'origine de ces éléments dans la matrice; la multi-

plication matricielle est seulement l'une de ces opérations. De façon encore plus spécifique, on observe à partir de ce qui

précède que la première opération dans le processus de multi-

plication consiste à multiplier le premier élément du vecteur de multiplication par un second vecteur contenant les éléments All, A21 et A31, les éléments de la colonne matricielle. Ainsi, la multiplication prend place en multipliant chaque élément du vecteur multiplication par un vecteur transposé comprenant les éléments d'une colonne de la matrice. L'élément vectoriel X2 est de même multiplié par le transposé de la seconde colonne de

la matrice contenant les éléments A12, A22 et A32. Les résul-

tats de chacune de ces multiplications comprennent des résul-

tats partiels qui, s'ils sont accumulés sur la base de la rangée d'origine de l'élément matriciel mémorisé, fourniront le vecteur résultant. Ceci se voit à partir de l'exposé précédent dans lequel Q1 est produit par l'accumulation des résultats partiels basés sur l'origine de la rangée, à savoir AllXl, A12X2, A13X3. Ceci reste vrai pour les éléments restants du

vecteur résultant.

Ayant noté que la multiplication matricielle ci-

dessus, entres autres, peut avoir lieu simultanément comme cela a 'été noté ci-dessus, il reste à ajouter l'observation suivante qui est applicable aux matrices creuses de grande dimension que l'on rencontre dans l'analyse à éléments finis. Pour de telles matrices creuses de grande dimension, des éléments nuls ne contribuent pas au résultat final de sorte qu'ils peuvent être négligés dans un processus de mappage de la matrice dans une mémoire. Alors que les éléments nuls sont toujours négligés dans le mappage de la matrice à manipuler dans un processeur parallèle, il peut également être possible de négliger d'autres éléments qui, bien qu'ils ne soient pas strictement nuls, sont

négligeables du point de vue des calculs. Par exemple, des élé-

ments tombant dans une gamme proche de zéro pourraient être ignorés une fois que l'on a déterminé que leur insertion ne change pas les résultats finaux des opérations. Ainsi, le terme

zéro doit être interprété dans ce sens.

Ayant noté que les opérations matricielles peuvent être réalisées d'une façon très parallèle, il reste à prévoir

un procédé de mappage d'une matrice dans un processeur paral-

lèle de façon à permettre une optimisation du traitement parallèle. Un tel processus de mappage est décrit dans l'art

antérieur et sera décrit en relation avec un système de traite-

ment parallèle comprenant un réseau de cellules de traitement dont chacune est associée avec sa propre mémoire et couplée pour communiquer avec l'autre, de façon générale comme cela est

représenté en figure 1.

La figure 1 représente un système systolique fonda-

mental dans lequel un ordinateur hôte 10 reçoit typiquement des données et fournit des résultats au réseau systolique de cellules de traitement en pipe-line 15 par l'intermédiaire d'un

module d'interface 12. L'hôte peut bien str comprendre un ordi-

nateur, une mémoire, un dispositif à fonctionnement en temps réel, etc et, alors que, comme cela est représenté en figure 1, l'entrée et la sortie sont couplées à l'hôte, l'entrée peut provenir d'un premier hôte physique alors que la sortie est orientée vers un autre. L'objet essentiel de l'architecture systolique est de réaliser des calculs multiples pour chaque

accès d'entrée/sortie du réseau avec l'hôte. Ces calculs mul-

tiples sont réalisés en parallèle en prévoyant que les cellules continuent le traitement pour une entrée reçue selon un mode "pipe-line". Ainsi, tandis que la cellule 1 reçoit et traite

une nouvelle donnée, les autres cellules continuent à fonc-

tionner sur des données ou des produits partiels de données reçus lors d'un accès d'entrée précédent. L'objectif de l'utilisation d'un réseau systolique est de décomposer le problème à résoudre en sous-étapes qui peuvent être traitées en parallèle par les diverses cellules du système de façon rythmique. Un processus de charge et de mémorisation d'une matrice creuse de grande dimension dans une série de cellules

de traitement en parallèle a été décrit dans l'art antérieur.

Ce processus est répété ici par un souci de description

complète de l'invention et sera mieux compris en relation avec l'organigramme du processus de la figure 2 en relation avec l'architecture de réseau systolique de la figure 1. Il faut noter que le chargement et la mémorisation d'une matrice sont réalisés sous la commande du module d'interface 12 de la figure 1 qui communique pendant ce processus avec l'hôte 10 qui a déjà mémorisé la matrice sous forme d'un réseau ou autre structure

de données appropriée.

* L'objet de la procédure de mappage de la figure 2 est, comme on y a fait allusion ci-dessus, de mémoriser la matrice dans la mémoire du multiprocesseur de façon à permettre

un haut degré de simultanéité dans la réalisation des opéra-

tions matricielles de façon générale et plus particulièrement dans l'opération de rétro-substitution. Ceci est réalisé en

26 1 0 4 29

mémorisant les éléments d'une ligne de la matrice d'une façon

qui permet une opération simultanée par un opérande commun.

Plus particulièrement, et comme ceci sera décrit ci-après, dans une opération de rétro-substitution, les éléments d'une colonne de matrice sont mémorisés pour une multiplication successive ou

simultanée par un opérande commun qui est une composante asso-

ciée et calculée précédemment du vecteur de la variable de

champ que l'on cherche à résoudre.

Une seconde caractéristique du processus de mémori-

sation est qu'il est produit un indice destiné à être associé à chaque élément matriciel mémorisé qui identifie l'origine (par rapport à une autre ligne de la matrice) de chaque élément. En particulier, chaque élément en provenance d'une colonne de la matrice est associé à un indice qui identifie la rangée de la matrice d'o il provient. En outre, les éléments mémorisés de chaque colonne sont formés en des vecteurs transposés dont chacun peut être associé à un opérande commun pour assurer une régularité et une simultanéité de fonctionnement. Ainsi, plusieurs vecteurs transposés sont formés, chacun étant situé à un emplacement donné de mémoire dans une pluralité de cellules,

de sorte qu'un opérande commun peut agir simultanément sur eux.

L'origine et le caractère de l'opérande commun dépendra dans une grande mesure de l'opération matricielle particulière effectuée. En outre, des éléments nuls et négligeables du point

de vue des calculs sont négligés dans le processus de mémori-

sation de sorte que les opérations matricielles impliquant les éléments mémorisés peuvent être effectuées rapidement. Enfin, les éléments diagonaux contenus dans tout vecteur transposé donné sont rassemblés pour mémorisation dans une cellule de processeur unique de sorte qu'ils sont facilement disponibles pour un traitement régulier et répétitif pendant certaines opérations matricielles, l'opération de rétro-substitution

étant une telle opération matricielle.

Les objectifs ci-dessus sont atteints en transposant

les éléments des colonnes matricielles en des vecteurs trans-

posés répartis entre les cellules de traitement. Puisque les éléments nuls de la matrice sont négligés pendant un processus de mémorisation, le nombre de cellules de traitement nécessaire pour réaliser le processus de mémorisation ci-dessus est tout à fait raisonnable. De nombreux systèmes physiques sont repré-

sentés par des matrices ayant de l'ordre de 25000 x 25000 élé-

ments dans lesquelles le nombre d'éléments non nuls dans toute colonne donnée est de l'ordre de vingt ou moins. Ainsi, par suite de l'opération de mémorisation ci-dessus, chaque élément non nul d'une colonne d'une matrice devient un élément d'un vecteur ou réseau transposé mémorisé par le même emplacement mémoire de chaque cellule de traitement. En outre, tous les éléments diagonaux sont mémorisés, pour un accès facile, dans une cellule de processeur unique. Les éléments non nuls des colonnes étant accessibles simultanément par ce processus, les opérandes communs peuvent être amenés à agir sur le vecteur transposé d'une façon extrêmement parallèle, comme cela sera

exposé plus en détail ci-après.

Un procédé de mappage d'une matrice dans un réseau multiprocesseur pour atteindre les objectifs ci-dessus va maintenant être décrit en relation avec la figure 2. D'abord, en partant à l'étape 200, la constante "C" est choisie égale au

nombre total d'éléments de traitement dans le réseau de proces-

seurs en parallèle utilisé. Le nombre de cellules de traitement peut varier, mais une efficacité maximale est obtenue quand ce nombre est approximativement identique au nombre maximal d'éléments non nuls dans une colonne unique de la matrice creuse de grande dimension à mapper. Les variables de processus "rangée" ("Row") (indice de rangée pour l'élément matriciel traité) et "Col" (indice de colonne pour l'élément matriciel traité) sont toutes deux choisies égales à 1, de sorte que le processus commence pour l'élément matriciel de la première rangée, première colonne et part de là. La variable "Proc" (traitement) (la cellule de traitement suivante devant être chargée par un élément non-diagonal non-nul) est choisie égale à 2; bien sûr, "Proc" peut varier entre 1 et C mais la cellule

1 a été étiquettée pour contenir seulement les éléments diago-

naux de la matrice comme cela sera décrit plus en détail ci-

après. La variable MémoireL (l'emplacement mémoire dans le "Proc" chargé) est choisie égale à "Mem", le premier empla-

cement mémoire à remplir dans chacune des cellules de proces-

seur. L'hôte fixe également ou a déjà fixé à partir d'un processus antérieur les variables M (le nombre de rangées dans la matrice) et N (le nombre de colonnes dans la matrice. Enfin, à l'étape 200, une table appelée Table[col], qui est placée dans l'hôte, est initialisée à des valeurs toutes à 1 lors d'une étape préliminaire pour corréler une colonne de la

matrice (un vecteur transposé spécifique à créer) avec un opé-

rande spécifique devant être utilisé ensuite dans une opération

matricielle, comme cela sera exposé plus en détail ci-après.

Ensuite, à l'étape 201, l'élément matriciel A[Row][Col] (qui, au début du processusest l'élément matriciel A[1] [1] est testé pour déterminer s'il est égal à zéro. Si c'est un zéro, l'élément ne doit pas être chargé dans une cellule quelconque du processeur et sera négligé. Le processus

passe alors aux étapes 202 et 203 auxquelles Row est incré-

menté, et la valeur de Row est testée pour déterminer si Row est supérieur à M. Sinon, le processus revient à l'étape 201 pour tester un nouvel élément matriciel. Si R est supérieur à M, une colonne complète de la matrice a été vérifiée et, à l'étape 205 des zéros sont mémorisés dans la ou les cellules restantes du processeur à partir de Proc vers C et à ou aux emplacements de mémoire Mem. Le but de cette étape, comme on y

a fait allusion précédemment, est de permettre à chaque empla-

cement de mémoire spécifique dans chaque processeur de contenir des éléments mémorisés associés à une seule colonne de la matrice. De cette façon, les éléments associés à une colonne donnée définiront un vecteur transposé (chacun mémorisé à un emplacement de mémoire connu sur une pluralité de cellules) et une pluralité de tels vecteurs transposés (1 pour chaque

colonne à la matrice) sera créée.

Après que toutes les cellules de traitement pour un Mem particulier ont été remplies, le processus continue à balayer la matrice en partant du haut de la colonne suivante. Ainsi, à l'étape 207, Row est ramené à 1 et Col est incrémenté pour préparer le test des éléments suivants dans une nouvelle colonne de la matrice. Toutefois, si, à l'étape 209, Col est supérieur à N (le nombre de colonnes de la matrice), alors le

processus de chargement de la matrice se termine (étape 210).

Si Col n'est pas supérieur à N, alors le processus revient à l'étape 201 au niveau de laquelle le nouvel élément matriciel

est testé.

Si le test à l'étape 201 indique que l'élément matriciel testé n'est pas nul, l'élément est testé à l'étape 230 pour déterminer si c'est un élément diagonal, c'est-à-dire

si Row = Col. Si c'est un élément diagonal, l'élément est mémo-

risé à Proc = 1 (cellule de processeur 1) et un indice est produit pour accompagner la diagonale mémorisée indiquant la rangée de la matrice d'o il provient. Après mémorisation d'une diagonale dans la cellule 1, le processus passe aux étapes 202

et suivantes d'une façon similaire à celle décrite ci-dessus.

Si l'élément matriciel testé à l'étape 230 n'est pas

un élément diagonal, il est mémorisé dans la cellule de proces-

seur Proc à l'emplacement de mémoire Mem, étape 220. Alors, un indice égal à Row est produit et mémorisé en association avec l'élément nouvellement mémorisé au même emplacement mémoire Mem, étape 221. Ayant mémorisé un élément non-nul, le processus incrémente Proc, étape 225, de sorte que l'élément non-nul suivant sera mémorisé dans la cellule de processeur suivante,

puis reste Proc à l'étape 227 pour déterminer s'il est supé-

rieur à C (le nombre total de cellules du processeur). S'il en est ainsi, ceci indique que tous les emplacements Mem dans les cellules de processeur ont été remplis et, à l'étape 228, Proc est remis à 2 et la mémoire est incrémentée à une nouvelle

valeur. En outre, l'entrée sur Table[col] dans la table précé-

demment créée est incrémentée de sorte que, pendant le traite-

ment ultérieur, l'hôte peut corréler chaque colonne de la matrice avec une longueur donnée d'un vecteur transposé. Cette entrée permet au module hôte ou au module d'interface d'associer les éléments non nuls correspondant à cette Col en cours à un vecteur transposé recouvrant plus d'un emplacement de mémoire (le nombre exact étant contenu dans la table). : Si des emplacements mémoire multiples sont requis, l'hôte prendra

ultérieurement une action appropriée dans le processus de fonc-

tionnement de la matrice pour assurer que les opérandes correspondants sont fournis dans les processeurs le nombre de fois convenable. Alors, le processus passe à l'étape 202, comme

cela a été précédemment exposé.

Si l'étape 227 détermine que Proc n'est pas supé-

rieur à C, le processus passe directement à l'étape 202 puisque des emplacements mémoire additionnels Mem existent dans les cellules de traitement. Le processus continue à partir de

l'étape 202 comme cela a été précédemment exposé.

On va maintenant décrire en relation avec la figure 3 un exemple de fonctionnement du processus de la figure 2 pour

mapper une matrice particulière [A] décomposée de façon trian-

gulaire supérieure qui a été maximisée en largeur de bande selon les enseignements de l'invention (représentée à la partie supérieure de la figure 3) dans des emplacements de mémoire, à

titre d'exemple, des cellules de processeur 1-4 d'un multipro-

cesseur similaire à celui de la figure 1. Comme on peut le voir, la matrice [A] décomposée est d'ordre 12 x 12 mais est de forme triangulaire supérieure. En particulier, la matrice ne comprend que des éléments nuls en dessous de la diagonale et tous ses éléments diagonaux sont non-nuls. En outre, toutes les données non-nulles dans la matrice [A] sont regroupées dans un groupe situé dans le coin supérieur droit et une bande de largeur K d'éléments nuls est située entre la diagonale et le groupement dans le coin supérieur droit. Ainsi, les données non-nulles dans la matrice sont forcées d'être ou bien sur la diagonale ou bien dans le coin supérieur droit, avec uniquement des éléments nuls dans l'intervalle. La technique de création d'une telle matrice à largeur maximisée et celle de choix de la valeur K sera décrite plus en détail ciaprès. Il faut garder en mémoire, quand on regarde ce mappage d'échantillons que la

matrice [A] est considéré comme un exemple d'une matrice trian-

gulaire supérieure maximisée en largeur de bande, de type creux et de grande dimension, que l'on rencontre couramment quand on caractérise des systèmes physiques par des procédés à éléments finis. Une telle matrice peut typiquement être de l'ordre de 25000 x 25000 avec seulement environ 0,1 % d'éléments non-nuls

dans toutes les rangées ou colonnes.

En commençant à l'étape 200 de la figure 2, les variables et les constantes sont initialisées pour se conformer aux valeurs particulières de la matrice [A] mappées de la façon suivante: C = 4; Row = 1; Col = 1; N = 12; M = 12; MémoireL = Mem; Proc = 2 est placé à 2 pour laisser de côté la

première cellule de processeur qui mémorise les éléments diago-

naux comme cela sera exposé plus en détail ci-après. Le premier élément A[Row][Col] est testé à l'étape 201. Puisque sa valeur

n'est pas nulle, il est ensuite testé à l'étape 230 pour déter-

miner si c'est un élément diagonal de la matrice, à savoir si Row = Col. Puisque l'élément [A] est un élément diagonal, il est mémorisé à un emplacement de la cellule de processeur 1 en même temps qu'un indice 1 désignant son origine à partir de la rangée 1 (étape 231). Le processus passe à l'étape 202 au niveau de laquelle Row est incrémenté. Alors, à l'étape 203, puisque Row est inférieur à M, le processus revient à l'étape 201 pour examiner l'élément suivant A[2][1]. Puisqu'il n'y a pas d'autres éléments non-nuls dans la première colonne de la matrice, le processus se déplacera de façon répétée par les étapes 201, 202, 203 jusqu'à ce que Row soit incrémenté pour être supérieur à 12, instant auquel le processus passera de l'étape 203 à l'étape 205. A l'étape 205, des zéros seront

mémorisés dans les cellules de processeur 2-4 à Mem pour indi-

quer que la colonne 1 de la matrice ne comprend pas d'éléments non nuls autre que l'élément diagonal A qui a été mémorisé dans la cellule 1. Le chargement de ces zéros dans les cellules 2-4 à Mem est représenté en figure 3. Ayant achevé le balayage et le chargement de la colonne 1 de la matrice [A], Col est incrémenté-à 2 et Row est remis à 1 pour commencer à analyser la colonne 2, si Col est

inférieur à N (étape 209).

Les colonnes 2-7 sont alignées de façon identique à la colonne 1 et, puisque chacune ne contient également qu'un seul élément non-nul B-G, respectivement, et puisque, en outre, ces éléments non-nuls sont des éléments diagonaux, ils sont mémorisés dans la cellule 1 aux emplacements Mem+l à Mem+6;

des zéros sont mémorisés à Mem+l dans les cellules 2-4.

Le processus repasse alors au balayage de la huitième

colonne de la matrice à l'étape 201 o A[8] [1] est testé.

Puisque sa valeur (H) (n'est pas nul) et puisque ce n'est pas un élément diagonal (étape 230), il est mémorisé à l'étape 220 dans la cellule 2 à Mem+7. Un indice est également mémorisé au

même emplacement dans la cellule 2 en même temps que "H", indi-

quant qu'il provient de la première rangée de la matrice.

Ensuite, Proc est incrémenté et testé aux étapes 225, 227 et, puisque Proc est inférieur à 12, Row est incrémenté et testé et un nouvel élément A[8] [2] est traité. Puisque cet élément est nul, le processus revient à nouveau à l'étape 201 pour tester A[8][3]. Puisque des zéros sont situés dans la colonne 8 dans les rangées 2-7, le processus effectue des boucles successives alors qu'aucune valeur n'est mémorisée. Puisque l'élément

A[8][8] est non nul et puisqu'en outre c'est un élément diago-

nal, il est mémorisé dans la cellule 1 à Mem+7 comme cela est représenté en figure 3, d'une façon identique à celle de l'lément de donnée précédent A dans la cellule 1. Puisque tous les autres éléments de la colonne 8 de la matrice sont des zéros, les cellules 3 et 4 sont remplies de zéros à Mem+7 par les étapes 205 pour préparer l'analyse de la colonne 9 de la matrice. Le processus continu, comme précédemment, mémorisant de façon itérative chaque élément diagonal non nul de la matrice dans des emplacements mémoire différents de la cellule 1, mémorisant des éléments non-nuls de chaque colonne dans des

cellules 2-4 tout en sautant les éléments nuls et en remplis-

sant les cellules de mémoire restantes dans tout emplacement de mémoire donné par des zéros avant de commencer à balayer une

nouvelle colonne. Le résultat en est que tous les éléments non-

nuls de la première colonne de la matrice [A] sont reformés selon un vecteur transposé dans l'emplacement mémoire identifié par Mem sur les cellules du processeur. En outre, tout élément diagonal inclus dans un tel vecteur transposé est situé dans une cellule de processeur spécifique, à savoir la cellule 1. En effet, l'emplacement mémoire Mem sert d'identificateur pour un vecteur nouvellement créé ou transposé comprenant tous les éléments non nuls de la colonne 1 de la matrice, l'élément diagonal étant mémorisé dans la cellule 1. Dans le cas de la colonne 1 de la matrice [A] de la figure 3, un seul élément y est contenu, tous les autres éléments étant de valeur nulle, les emplacements Mem restants sont remplis de zéros, comme cela

est représenté en figure 3. L'importance de ce vecteur trans-

formé nouvellement créé Mem et le traitement particulier des éléments diagonaux pour permettre le traitement simultané des

éléments matriciels seront brièvement rappelés ci-après.

Comme cela a été indiqué ci-dessus, le nombre de cellules de traitement est choisi pour être égal au nombre maximal d'éléments non nuls de toute colonne donnée de la matrice [A], mais peut être de 2 ou de tout nombre supérieur à 2 sans sortir du domaine de l'invention. Toutefois, quel que soit le nombre de cellules de processeur utilisées, il peut arriver que le nombre d'éléments non nuls dans une colonne particulière de la matrice dépasse le nombre de cellules du processeur. Pour traiter ce cas, les étapes 227 et 228, comme on y a fait allusion ci-dessus, sont incluses dans le processus de mappage tel que représenté dans l'organigramme de la figure 2. Si, par exemple, la première colonne de la matrice [A] avait

plus d'éléments non-nuls que de cellules de processeur dispo-

nibles, c'est-à-dire si Proc>C (étape 227), avant que tous les éléments d'une colonne aient été mémorisés, alors Proc est remis à 2 et Mem est incrémenté vers Mem+l (étape 228) et le vecteur transposé identifié par "Mem" est poursuivi ou étendu dans l'emplacement de mémoire Mem+l. Le résultat en est que la colonne 1 de la matrice est transformée en un vecteur transposé situé aux emplacements de mémoire Mem et Mem+l. En ce cas, l'entrée Table[Col] à la colonne 1 serait de "2" pour indiquer

ce fait.

Toutes les cellules restantes du processeur corres-

pondant à un emplacement de mémoire donné qui ne sont pas remplies par des éléments non-nuls après le test de tous les éléments à partir d'une colonne donnée de la matrice sont remplies de zéros selon l'étape 205. Par exemple, à la colonne de la matrice [A], il y a seulement trois éléments non-nuls, M, N et P, laissant Mem+9 de la cellule de processeur No4 à remplir par des zéros pour compléter le vecteur transformé à Mem+9. Une situation similaire est représentée en figure 3 par le remplissage de l'emplacement de mémoire Mem+10 de la cellule de processeur N04 par un zéro par suite de l'étape 205 puisque seulement trois éléments non-nuls (S, R, T) sont présents à lacolonne 4 de la matrice.

La fonction de la partie des emplacements mémoire

désignée par Q' en figure 3 sera expliquée ci-dessous en rela-

tion avec la réalisation d'une opération de rétro-substitution pour résoudre un système d'équations linéaires impliquant la

matrice [A] de la figure 3.

La réalisation de l'opération de rétro-substitution sur le système décomposé d'équations caractérisé par la matrice décomposée [A] à largeur de bande maximisée selon l'invention va maintenant être décrite en relation avec la figure 4. En bref, il faut comprendre que, avant l'initialisation de l'opération de rétro-substitution, un système d'équations qui décrit le produit ou système physique en cours d'investigation est produit sous la forme:

[K] [Y] = [R]

Comme cela sera rappelé à partir de l'exposé des pages 7 et 8 ci-dessus, [Y] est un vecteur qui représente une variable de champ inconnue qui décrit un attribut du système physique en cours d'investigation, tel qu'un déplacement pour un système à ressorts. Le vecteur [Y] est constitué de composantes Yi...YN pour lesquels le système d'équations doit être résolu. Le vecteur [R] comprend d'autre part des composantes Rl.

RN qui..DTD: sont connues et représentent une variable vectorielle résul-

tante, telle qu'une force dans un système de ressorts analysé.

[K] est la matrice de raideur et comprend une pluralité de

valeurs ou de constantes qui concernent les variables vecto-

rielles résultantes connues pour la variable vectorielle de champ inconnue à des noeuds spécifiques du réseau d'éléments

finis, comme cela a été exposé précédemment.

Selon l'invention, le système d'équations susmen-

tionné est transformé ou décomposé en une autre série d'équations:

[A] [X] = [Q]

o [A] est sous la forme d'une matrice dite maximisée en largeur de bande, comme indiqué ci-dessus, [Q] étant un vecteur connu et [X] un vecteur inconnu. Le système étant transformé sous cette nouvelle forme maximisée en largeur de bande, la solution selon l'invention peut être traitée sur un processeur parallèle du type représenté en figure 4 et, en utilisant (au

moins pour certaines parties du processus de solution) la tech-

nique de rétro-substitution ou d'élimination directe, est mis en oeuvre sur un réseau systolique. L'objet du processus de solution est en conséquence de calculer les valeurs de X1...XN par la technique de rétrosubstitution, étant donné que [A] et [Q] sont connus. Dans sa forme fondamentale, et comme cela a

été mentionné précédemment, la matrice décomposée A est typi-

quement de grande étendue et creuse. Spécifiquement, elle peut être de l'ordre de 25 000 x 25 000 éléments, le nombre d'éléments non-nuls étant seulement un très petit pourcentage

du nombre total d'éléments. Elle est également de forme trian-

gulaire supérieure, maximisée en largeur de bande, c'est-à-dire d'une forme dans laquelle tous les éléments diagonaux de la matrice sont nonnuls, tous les autres éléments non nuls (y compris certains éléments nuls) étant groupés dans le coin supérieur droit de la matrice; et une bande de largeur K est située au-dessus de la diagonale contenant seulement des données de valeurs nulles, comme cela est représenté en figure 5B. La matrice [A] de la figure 3 est précisément une telle

matrice triangulaire supérieure maximisée en largeur de bande.

Ayant fourni une matrice décomposée de forme triangulaire supérieure maximisée en largeur de bande, la première étape dans le processus de solution est de charger et de mémoriser la matrice [A] comme cela a été décrit ci-dessus en relation avec les figures 2 et 3. Le résultat du processus de chargement et de mémorisation de la matrice [A] de la figure 3 est représenté en figure 4 dans laquelle la mémoire 56 de la

cellule 1 est représentée comme mémorisant les éléments diago-

naux A...Y dans la mémoire 58 à Mem --> Mem+11, respectivement, chaque valeur diagonale mémorisée étant accompagnée d'un indice mémorisé dans la partie de mémoire 59 désignant l'origine de la rangée de son élément diagonal associé. Une mémoire de résultat

partiel 57 est également prévue dans la cellule 1 pour mémo-

riser des valeurs calculées Xl...X12 du vecteur X et son fonc-

tionnement va être expliqué plus en détail ci-après.

De façon similaire, les mémoires 56 des cellules 2 à 4 sont représentées comme mémorisant les valeurs non-nulles restantes de la matrice, chaque élément non-nul mémorisé contenu dans la partie de mémoire 58 étant accompagné d'un indice mémorisé dans la partie 59 pour identifier la rangée de la matrice d'o son élément d'accompagnement non-nul est provenu. Chacune des cellules 2-4 comprend également une mémoire 57 pour mémoriser les résultats partiels accumulés en accord avec l'indice mémorisé accompagnant l'élément matriciel

mémorisé, comme cela sera exposé plus en détail ci-après.

Il faut également noter que les éléments de chaque

colonne de la matrice [A] ont été réordonnés pendant le proces-

sus de mappage pour former ce que l'on peut appeler des "vecteurs transposés" à chacun des emplacements mémoire Mem à Mem+11, les éléments diagonaux étant rassemblés dans la cellule 1. Par exemple, la colonne 9 de la matrice A est transformée en un vecteur transposé à l'emplacement mémoire Mem+8, l'élément diagonal L de la rangée 9 de la matrice étant situé dans la cellule 1 et les éléments non nuls restants de la colonne 4 (J, K) contenus dans les cellules 3 et 2, respectivement, avec des indices (1,2), respectivement, indiquant les rangées de la

matrice d'o ils proviennent. Comme cela a été exposé ci-

dessus, une valeur de données nulle est utilisée pour compléter

les emplacements mémoire non remplis dans chaque vecteur trans-

posé restant après que tous les éléments non nuls dans sa

colonne associée ont été mémorisés.

Etant donné la matrice triangulaire [A] de la figure 3, le système d'équations à résoudre pour [X] est le suivant: (1) AX1+...... HXS+JX9 = Qi

(2) BX2±-...... KX9+MX10+ UX12 = Q2

(3) CX3+...... NXO10+RXll+... = Q3 (4) DX4+...... SXll+VX12 = Q4

(5) EX5+... 0... WX12 = Q5

(6) FX6+...... = Q6

(7) GX7+...... = Q7

(8) IX8+......= Q8

(9) LX9+...... = Q

(10) PX1+...... = Q1o (11) TXll+... = Qil

(12) YX12 = Q12

Comme on se rappelera à partir de l'exposé des pages 10 et 11 ci-dessus, pour résoudre ce système d'équations par le processus de rétrosubstitution, X12 est d'abord calculé en résolvant l'équation (12). De même, de façon similaire Xll, X10, X9, Xg, X7 et X6 peuvent être résolus immédiatement par une simple division d'une composante connue associée de [Q] par

un élément diagonal correspondant T, P, L, I, G et F, respec-

tivement. I1 faut noter que la puissance de l'invention repose dans le fait que plusieurs composantes de la variable de champ inconnue [X] peuvent être calculées simultanément et sans relation avec d'autres valeurs de la variable de champ. Ceci permet une augmentation importante de la vitesse de calcul sur un processeur parallèle comme cela sera exposé plus en détail ci-après. Ceci est en un contraste évident avec le procédé de fonctionnement décrit dans l'art antérieur dans lequel, sauf pour le calcul d'une composante unique de la variable de champ, les valeurs restantes doivent être calculées en relation avec les valeurs successives déjà connues. Ce problème de dépendance des données réduit considérablement la vitesse de solution des équations par rétrosubstitution sur un processeur parallèle,

comme pourra le noter l'homme de l'art.

Ayant calculé les valeurs de X12-X8 par une simple division, comme cela a été souligné ci-dessus, ces valeurs connues et le résultat des opérations impliquant ces valeurs connues sont alors rendus disponibles pour calculer les valeurs restantes pour les composantes inconnues de [X], à savoir X5-X1 dans les équations (5) à (1), respectivement, page 32. Il faut également noter ici que le système d'équations ci-dessus constitue seulement un exemple d'un système considérablement plus grand et plus creux que l'on rencontre dans des problèmes physiques typiques étudiés en utilisant un procédé d'analyse à

éléments finis. Le système d'équations de la page 32 est uti-

lisé ici dans un but illustratif seulement et pour expliquer

l'avantage de l'invention découlant de la réduction des dépen-

dances de données dans le processus de solution en utilisant une matrice de solution maximisée en largeur de bande du type

représenté en figure 5B.

Le processus réel pour trouver la solution du

système d'équations ci-dessus par la technique de rétro-substi-

tution dans un système multiprocesseur systolique va maintenant être exposé plus en détai: en relation avec la figure 4 et les figures 4A-4M. Le système multiprocesseur comprend 4 cellules de processeur, les cellules 1 à 4, dont chacune est affectée à

la réalisation d'opérations spécifiques contribuant à la solu-

tion. Chacune des opérations réalisées par une cellule de

processeur est illustrée par un module fonctionnel en figure 4.

Il faut garder à l'esprit que les modules fonctionnels auxquels

on se réfère en figure 4 ne désignent pas des composants maté-

riels distincts d'une machine mais plutôt des opérations de traitement de données différentes pouvant être réalisées par chaque cellule du système, comme cela sera facilement compris de l'homme de l'art. Les opérations réalisées par chaque cellule sont conçues pour être régulières et répétitives

pendant chaque cycle de la machine, mais réalisées sur des opé-

randes différents présentés aux entrées de chaque module pendant des cycles successifs de la machine. Le débit

2610 4 2 9

d'informations à travers le processeur va maintenant être expliqué en relation avec les figures 4A-4M qui indiquent la circulation des opérandes et des résultats au travers de la machine pendant ses 12 premiers cycles dans un processus de rétro-substitution de l'exemple d'ensemble d'équations ci- dessus. Le système multiprocesseur de la figure 4 comprend une mémoire 56 associée à chacune des cellules 1-4. Les mémoires 56 peuvent être considérées comme comprenant un premier réseau 58 qui mémorise les valeurs des éléments de la

matrice de solution A, comme cela a été exposé précédemment.

Chaque élément mémorisé dans le réseau de mémoire 58 est asso-

cié à un indice mémorisé dans le réseau 59 qui désigne la

rangée de la matrice d'o provient sa valeur mémorisée asso-

ciée. Ainsi, en considérant la cellule 1, l'élément A de la matrice provient de la rangée 1 de la matrice, l'élément L de la matrice de la rangée 9 de la matrice, etc. Une mémorisation dans la cellule 1 signifie également que l'élément mémorisé est un élément diagonal sur la base de la technique de mémorisation

de la figure 2.

Chaque cellule contient également un autre module de mémoire 57 qui mémorise des résultats en provenance de calculs indexés réalisés par le système, comme cela sera expliqué plus

en détail pendant la description du processus de calcul ci-

après. Le module de mémorisation 57 de la cellule 1 fonctionne comme la mémoire pour les valeurs composantes finales X1 à X12 de la variable de champ inconnue. Les modules de mémorisation 57 des cellules 2, 3 et 4 mémorisent des résultats selon les indices et, au lieu de laisser entrer et sortir des valeurs de façon séquentielle, répondent à des entrées indexées pour

fournir des sorties indexées, c'est-à-dire des valeurs mémo-

risées associées à ces indices. Dans ce but, les modules

mémoire 57 des cellules 2-4 fournissent deux sorties séparées.

Une première sortie sur la ligne 57' est envoyée à une série

d'accumulateurs 53 au niveau desquels différentes valeurs asso-

ciées au même indice sont ajoutées les unes aux autres avant

d'être restaurées en mémoire dans un ordre basé sur cet indice.

La sortie sur la ligne 57' représente la valeur mémorisée dans la mémoire associée à l'indice de la diagonale en cours de traitement dans l'accumulateur 53. Pour cette raison, la sortie sur la ligne 57' est commandée par une entrée d'indice sur la ligne 59' dans la mémoire 59. Une seconde sortie sur une ligne 57" en provenance de chacune des mémoires 57 des cellules 2-4 est fournie à un additionneur 65. Les sorties indexées sur la

ligne 57" représentent des accumulations Q'N de calculs réali-

sés antérieurement qui doivent être sommés pour mettre en oeuvre le processus de rétro-substitution. Les sorties indexées sur la ligne 57" sont choisies par des sorties simultanées pour chacune des mémoires 57 des cellules 2-4 en provenance d'un compteur 61 qui décompte en partant à l'indice 11 pendant le

cycle 1 jusqu'à l'indice 1 du cycle 11, comme cela est repré-

senté dans les dessins joints aux figures 4A-4M.

Les valeurs mémorisées dans les mémoires 58 des cellules 1-4 sont fournies à des modules fonctionnels dans le

reste du système associé à chaque cellule. Les valeurs mémo-

risées dans les mémoires 58 sont envoyées à partir des adresses supérieures extrêmes ou de numéros élevés (comme on le voit en figure 4) vers les modules fonctionnels 52 et 51 par des lignes 58' de façon étagée. Ceci signifie que la valeur mémorisée à

Mem+ll dans la mémoire 58 de la cellule 1 est fournie au divi-

seur 52 lors d'un cycle de machine donné qui est, respecti-

vement, antérieur de 1, 2 et 3 cycles au moment o les valeurs correspondantes pour les mêmes emplacements mémoire dans les

mémoires 58 des cellules 2, 3 et 4 sont fournies à leurs multi-

plieurs associés 51 par des lignes 58', comme cela sera exposé plus en détail ci-après. La cellule 1 contient un module de soustraction 60 et un module de division 52. Le module de soustraction 60 est alimenté par un premier opérande sur une ligne 65' par un module additionneur 65. La seconde entrée vers le soustracteur 60 est la valeur de l'une des composantes

2 6 1 0 4 2 9

vectorielles connues dans le système physique Q1-Q12- Le module de division 52 divise le résultat du processus de soustraction mis en oeuvre par le module 60 par une valeur appropriée en provenance du module de mémorisation 58 de la cellule 1, qui est un élément diagonal choisi de la matrice mémorisée. Chaque résultat sortant du diviseur 52 de la cellule 1, représente en conséquence une solution d'une composante de la variable de champ inconnue X dans le système d'équation résolu. Ces résultats, selon le procédé de rétro-substitution décrit ci-dessus, sont mémorisés pour une lecture ultérieure dans le module de mémoire 57 de la cellule 1. En outre, chaque résultat de la division dans le module 52, doit être multiplié par le vecteur transposé représentant la colonne correspondante de la matrice. Par exemple, X12 doit être multiplié par chaque élément non-nul en provenance de la colonne 12 de la matrice; les résultats de ces opérations doivent être rassemblés et ensuite sommés selon l'indice de rangée. Dans ce but, chaque

résultat en provenance du module 52 est envoyé en mode pipe-

line au module de multiplication 51 des cellules 2-4 au niveau desquelles un résultat indexé est calculé. Les sorties de chacun des modules 51 sont accumulées selon l'indice de rangée dans des modules accumulateurs 53 dont les résultats sont à nouveau mémorisés selon l'indice de rangée dans des modules de mémorisation 57 pour être fournis ultérieurement à l'additionneur 65 sur la ligne 57" sous la commande du compteur 61. Le système de la figure 4 sera mieux compris en

expliquant son fonctionnement lors de la résolution des équa-

tions simultanées de la page 32 par le processus de rétro-

substitution. Cette explication sera associée aux figures 4A-4M qui représentent le flux de données dans le système dans ce but et pour effectuer la méthodologie de rétro-substitution, telle

qu'on l'a souligné en page 32 ci-dessus.

De façon générale, le système de processeur de la figure 4 fonctionne en pipe-line. Au début du processus de

2610 4 2 9

calcul, seules les valeurs contenues dans les modules de mémoire 56 et les valeurs connues QI à Q12 sont disponibles,

les valeurs QI à Q12 étant mémorisées dans un processeur asso-

cié, tel que l'hôte 1, figure 1, pour introduction de façon appropriée dans le module 60. Lors de cycles successifs de la machine, certaines sélectionnées de ces valeurs sont rendues disponibles pour les cellules multiprocesseur pour produire des opérandes résultants. Le processus commence par le déplacement des données de façon générale de la gauche vers la droite en

figure 4 à travers les processeurs pour les remplir progres-

sivement jusqu'à ce que chaque cellule effectue un calcul significatif pendant chaque cycle. Les quelques premiers cycles de la machine servent à remplir le pipe-line de la façon décrite ci-dessous tout en calculant simultanément les quelques premières valeurs pour les composants de la variable de champ

inconnue [X].

Comme le représente la figure 4A, pendant le premier cycle de la machine, Q12 est fourni à partir d'un module de mémoire approprié ou d'un ordinateur hôte (non représenté) vers une entrée du soustracteur 60; l'autre entrée vers le module , Q'12, en provenance de l'additionneur 65 est à zéro à cet instant dans le processus puisqu'aucun des opérandes à traiter n'a atteint les autres modules fonctionnels du système à ce stade dans le processus de rétro-substitution, comme cela sera exposé ci-après; les entrées de tous les autres modules fonctionnels sont égales à zéro. La sortie Q12 du soustracteur pendant le cycle 1 est rendue disponible pour le diviseur 52

pendant le cycle 2.

Pendant le cycle 2 de la machine,- figure 4B, l'équation (12) de la page 32 est résolue pour X12. Ceci est effectué en divisant la sortie Q12 du soustracteur 60 par le premier élément diagonal Y de la matrice [A] pour fournir la première solution X12, à la sortie du module 52. Simultanément, pendant le cycle 2, la sortie Qll' de l'additionneur 65 (qui est encore à zéro) est soustraite de Qll dans le module 60 pour fournir une nouvelle entrée au diviseur 52 pour le cycle suivant. Alors, X12 est mémorisé par la ligne 64 dans la mémoire 57 et fourni sur la ligne 63 au module multiplicateur

51 de la cellule 2.

Pendant le cycle 3, figure 4C, l'équation (11) est résolue pour Xjl en divisant Qll par T et en le mémorisant par la ligne 64 dans la mémoire 57 de la cellule 1. Egalement, pendant le cycle 3, la valeur X12 précédemment calculée à la

ligne 63 pour la cellule 2 est multipliée par W (l'un des élé-

ments de la matrice en provenance de la douzième colonne de la matrice) dans le module 51 de la cellule 2 pour produire le produit partiel W.X12 à une entrée de l'accumulateur 53 de la cellule 2. Egalement, pendant le cycle 3, la sortie de l'additionneur 65 (qui est encore à zéro) est soustraite de Q1o

dans le soustracteur 60 de la cellule 1.

De façon analogue, pendant le cycle 4, l'équation (10), page 32, est résolue pour X10 en divisant Q1o par

l'élément diagonal P mémorisé ensuite dans le diviseur 52.

Simultanément, au niveau de la cellule 2, la valeur Xjl précé-

demment calculée est multiplié par S pour produire un produit partiel S. Xll, alors que le produit partiel préalablement produit W.X12 est accumulé avec d'autres produits partiels Q's, s'il y a lieu, (sur la ligne 57') dans l'accumulateur 53 de la cellule 2 pour mémorisation lors du cycle suivant en tant que Q'5 dans le module de mémoire 57 de la cellule 2 par la ligne 79. Egalement, X12 qui est passé par le retard d'un cycle 55 a maintenant atteint la cellule 3 et est multiplié par l'élément

matriciel V pour produire un produit partiel pour une accumu-

lation et une mémorisation ultérieures lors du cycle suivant

dans la mémoire 57 de la cellule 3.

Pendant le cycle 5, la sortie Q'8 en provenance de l'additionneur 65 (qui égale encore zéro) est soustraite de Q8 dans le module soustracteur 60; Q9 (la sortie du soustracteur lors du cycle précédent) est divisée par L dans le diviseur 52 pour produire X9; la valeur X10 précédemment calculée est mémorisée dans la mémoire 57 de la cellule 1 et est également multipliée dans la cellule 2 par N; la sortie précédente S.Xl du multiplieur 51 de la cellule 2 est accumulée avec des valeurs indexées et mémorisées de façon similaire (s'il y en a) de la mémoire 57 (Q'4 à la ligne 57') au niveau de l'accumulateur 53 pour mémorisation dans le module 57 lors du

cycle suivant; la valeur Xll préalablement calculée est multi-

pliée par R au niveau du module multiplieur 51 de la cellule 3

et V.X12 est accumulé avec des valeurs indexées de façon simi-

laire et préalablement mémorisées Q4' (s'il y en a) dans

l'accumulateur 53 de la cellule 3 pour une mémorisation ulté-

rieure dans la mémoire 57 lors du cycle suivant. Enfin, X12 s'est déplaçé vers la cellule 4 pour être multiplié par U dans le multiplieur 51, le produit M.X12 est fourni à l'accumulateur 53 de la cellule 4 pour accumulation et mémorisation pendant le

cycle suivant.

Les cycles ultérieurs 6-13 de la machine pour mettre en oeuvre le processus de rétro-substitution agissent d'une façon analogue à celle décrite ci-dessus. Les figures 4F-4M indiquent en détail la circulation des données à travers la machine pour résoudre les equations suivantes. En bref, lors de chaque cycle suivant une autre des équations de la page 32 est résolue pour Xl2...Xl, la valeur de cette composante étant mémorisée dans la mémoire 57 de la cellule 1 en même temps que les valeurs précédemment calculées. Chaque valeur nouvellement calculée est également successivement transmise vers la droite, comme on le voit en figure 4, vers chacun des étages ultérieurs du processeur pour calcul des produits partiels QN' qui sont accumulés par indice selon la rangée d'origine des éléments de la matrice associés à l'opération. Par exemple, pendant le cycle 6, figure 4F, trois produits partiels différents sont simultanément accumulés pour une sommation ultérieure dans l'additionneur 65 de la façon suivante: (1) un premier produit partiel N.X10 est accumulé

* (en même temps que tout autre produit partiel précédent parta-

geant le même indice de rangée que l'élément N) dans la mémoire 57 de la cellule 2. Ainsi N.X1o est mémorisé en tant que Q3' puisque l'origine de l'élément N provient de la rangée 3 de la matrice. (2) un second produit partiel R.Xîl est également accumulé dans la mémoire 57 de la cellule 3 sous l'indice 3 en

tant que Q3'.

(3) De même un troisième produit partiel V.X12 est accumulé en même temps que les autres produits partiels présentant une origine de rangée commune en tant que Q2' parce que l'élément matriciel 4 provient de la rangée 2 de la matrice. Les deux premiers des produits partiels susmentionnés restent mémorisés dans les mémoires 57 jusqu'au cycle de la machine pendant lequel le compteur 61 accède à tous

les produits partiels associés à la rangée 3 de la matrice.

Ceci prend place au cycle 9, figure 4I, pendant lequel tous les produits partiels précédemment accumulés de la rangée 3, Q3', sont fournis au module de sommation 65. Pendant le cycle 10, la sommation (N.X10 + R.Xll) est disponible pour être soustraite dans le module 60 de la valeur connue Q3 de sorte que, pendant le cycle 11, la différence Q3 - Q3' peut être divisée par l'élément diagonal D pour produire X3. On voit ainsi que l'équation (3), page 32 peut être résolue pendant le cycle 11

seulement parce que pendant les cycles antérieurs une accumu-

lation a été effectuée sur la base de l'indice de rangée de tous les produits des valeurs précédemment calculées des

vecteurs inconnus et des éléments de matrice dans la rangée 3.

Plusieurs caractéristiques importantes du procédé

selon l'invention apparaîtront à l'homme de l'art. Une caracté-

ristique essentielle est qu'elle permet la solution des 12 équations simultanées de la page 32 en 13 étapes. Ceci est seulement possible en raison de la suppression des dépendances de données résultant de la forme maximisée en largeur en bande de la matrice utilisée dans la solution en combinaison avec la configuration spécifique de mémorisation et d'accès aux valeurs de la matrice dans un multiprocesseur tel que décrit ici. Plus particulièrement, pendant les cycles de la machine o les quelques premières valeurs du vecteur inconnu X sont calculées dans une cellule par une simple étape de division (rendue pos- sible par la forme maximisée de la matrice), des opérations matricielles sont effectuées simultanément dans d'autres cellules de processeur pour produire des produits partiels (QN') utilisant les composantes précédemment calculées. Avec un tel traitement, quand les produits partiels qui dépendent des composantes d'un vecteur inconnu précédemment calculé sont nécessaires, elles sont disponibles pour des calculs. Ainsi, le processus de rétro- substitution prend place sans avoir à attendre jusqu'à ce que de tels produits partiels soient

disponibles.

Comme on y a fait allusion précédemment, la largeur de bande minimale K de la bande de zéros interposée entre la diagonale principale de la matrice et le groupement dans le coin supérieur droit des éléments non- nuls est choisie sur la base de la conception du système multiprocesseur. De façon générale, K doit être égal au nombre maximal de cycles ou d'étapes que cela prend pour que l'agencement de processeur spécifique utilisé rassemble et rende disponibles les données nécessaires pour calculer une quelconque composante Xi d'un vecteur inconnu. En ce qui concerne le système de la figure 4, on voit que, une fois que X12 par exemple a été calculé pendant le cycle 2, six étapes ou cycles supplémentaires (3-8) sont nécessaires pour achever le calcul, l'accumulation et la sommation des produits partiels associés à X12 et pour fournir ces produits au module soustracteur 60 à temps pour qu'ils soient utilisés à la solution de l'équation (5). Ainsi, une largeur de bande K égale à 6 est suffisante pour fournir un espace et une durée propre à tirer complètement avantage du

réseau systolique particulier de la figure 4. La bande particu-

lière nécessaire pour d'autres systèmes matériels dépendrait bien sr des particularités d'un tel système. Si une largeur de bande plus faible est prévue, une certaine efficacité peut être sacrifiée en fonction de la composition des données dans la

série d'équations sur lesquelles on agit.

Comme cela sera généralement compris dans la tech- nique, la caractéristique de groupement et le caractère creux reflètent la connectivité de la maille à éléments finis. Ainsi, il est connu dans la technique que transformer une matrice de

système [K] en une matrice [A] triangulaire supérieure, maxi-

misée en largeur de bande, ayant une structure creuse, comme cela est représenté en figure 5B est bien connu. Alors que la minimisation de la largeur de bande telle que représentée en figure 5A implique unerenumérotation des noeuds de la maille à éléments finis de sorte qu'une différence minimale apparaisse entre deux noeuds quelconques adjacents, la maximisation de la matrice de la figure 5B est réalisée en renumérotant les noeuds

de système pour maximiser la différence entre noeuds adjacents.

Ceci a été réalisé par l'invention sans avoir recours à une durée indue par une technique d'approximations successives dans laquelle une première passe utilise des numéros aléatoires assignés aux noeuds. Puisque la maille typique est grande, cette technique supprime habituellement 90 % des données non nulles pour la bande adjacente à la diagonale. Un raffinement supérieur est réalisé en commutant des numéros entre noeuds d'éléments nuls dans et à l'extérieur de la bande. Le fait de

répéter plusieurs fois cette technique d'approximations succes-

sives entraîne presque toujours une largeur de bande suffisante pour fournir avec efficacité la matrice décomposée selon la

technique décrite ici. En outre, puisque seul un nombre relati-

vement faible de processeurs est envisagé, la dimension de la bande est généralement petite par rapport au nombre total

d'éléments finis d'une maille de système.

Comme le notera l'homme de l'art, en utilisant la matrice maximisée en largeur de bande [A] telle que prévue ici, une pluralité de composantes inconnues XN du vecteur de variable de champ inconnu peut être calculée sans recourir à d'autres valeurs de la variable de champ. Si ces calculs sont effectués séquentiellement dans un processeur unique (cellule 1 dans l'exemple ci-dessus) d'autres parmi des processeurs (cellules 2-4) du système peuvent, pendant la même durée, calculer et accumuler les produits partiels QN' qui seront nécessaires à un instant ultérieur pour parachever le calcul de

toutes les composantes de XN (en particulier celles qui dépen-

dent des valeurs précédemment calculées). Ainsi, en augmentant le nombre de calculs pour XN qui peuvent être effectués sans dépendances de données, du temps est assuré par un système à multiprocesseur pour parachever des calculs intermédiaires sans temps de repos intermédiaire. Tous les processeurs du système, au lieu d'être forcés d'attendre pendant une certaine partie du processus de rétro-substitution, sont autorisés à effectuer des

calculs significatifs pendant chaque cycle de la machine.

Alors que le processus de calcul matriciel concerné spécifiquement ici est la rétro-substitution, le processus d'élimination directe pourrait prendre place de façon similaire et pourrait facilement être mis en oeuvre selon les principes de l'invention par l'homme de l'art. Le processus d'élimination directe nécessite cependant que la matrice initiale du système soit transformée en forme triangulaire inférieure plutôt qu'en

forme triangulaire supérieure utilisée en rétro-substitution.

La technique d'élimination directe est également décrite en

détail dans les documents susmentionnés.

Il faut également souligner à nouveau que, alors que la technique selon la présente invention a été expliquée par souci de brièveté en ce qui concerne une petite matrice et un petit système d'équations, les avantages de l'invention sont plus évidents quand ils sont appliqués à une opération sur des matrices creuses de grande dimension qui sont associées à des systèmes d'équations creuses de grande dimension représentatifs de systèmes physiques typiques devant être analysés selon le

procédé d'analyse à éléments finis.

Alors que la présente invention a été décrite en relation avec des modes de réalisation particuliers incluant l'utilisation d'une opération de rétro-substitution réalisée par une architecture d'ordinateur à réseau systolique, l'invention n'est pas limitée à un quelconque fonctionnement

particulier matriciel à multiprocesseur parallèle. En consé-

quence, des modifications apparaîtront à l'homme de l'art ayant pris connaissance de l'invention sans sortir du domaine de celle-ci.

Claims (11)

REVENDICATIONS

1. Procédé d'analyse d'un phénomène physique carac-

térisé en ce qu'il comprend les étapes suivantes: produire un premier système d'équations linéaires pour décrire le phénomène physique, ce premier système d'équation étant caractérisé par une matrice de système connue K d'ordre N x N, un vecteur d'ordre N inconnu Y et un vecteur d'ordre N connu R, cette matrice de système associant le vecteur inconnu au vecteur connu par l'expression matricielle

KY = R,

cette matrice de système étant de grande dimension et creuse; transformer le premier système en un second système d'équations caractérisé par une matrice triangulaire décomposée

ne comprenant que des éléments non-nuls sur sa diagonale prin-

cipale, une bande d'éléments nuls s'étendant sur une largeur de bande présélectionnée à partir de la diagonale principale, une absence d'éléments non-nuls dans ladite bande, alors que tous les éléments nonnuls restants sont situés dans la partie de la matrice décomposée externe à ladite bande; charger des élénents de la matrice décomposée dans

des mémoires associées à une pluralité de processeurs, les élé-

ments diagonaux de la matrice décomposée étant contenus dans la mémoire associée à l'un des processeurs; et résoudre le seconc système d'équations en calculant les valeurs d'une pluralité de composantes du vecteur inconnu d'ordre N indépendamment de toute valeur préalablement calculée d'une composante quelconque de ce vecteur inconnu d'ordre N.

2. Procédé selon la revendication 1, caractérisé en ce que l'étape de résolution comprend l'étape consistant à mettre en oeuvre le processus de rétro-substitution ou d'élimination directe du second système pour déterminer les composantes du vecteur inconnu d'ordre N par une succession de manipulations matricielles effectuées simultanément dans une

pluralité de processeurs en parallèle.

26 1 0 4 2 9

3. Procédé selon la revendication 2, caractérisé en ce que l'étape de résolution de chaque composante du vecteur

inconnu d'ordre N comprend une succession d'opérations effec-

tuées par une pluralité de processeurs.

4. Procédé selon la revendication 3, caractérisé en ce qu'il comprend l'étape consistant à mémoriser des éléments de données de la matrice décomposée dans les mémoires associées

à une pluralité de processeurs.

5. Procédé selon la revendication 1, caractérisé en ce que l'étape de transformation comprend les étapes suivantes: exprimer le phénomène en termes d'une pluralité d'éléments finis dont chacun est défini par une pluralité de noeuds selon des techniques d'analyse à éléments finis classiques pour produire le premier système d'équations; et renuméroter une pluralité desdits noeuds de façon que des noeuds adjacents diffèrent d'une valeur minimale présélectionnée, pour assurer ainsi que le second système

d'équations linéaires est caractérisé par une matrice décom-

posée du type indiqué en revendication 1.

6. Procédé selon la revendication 5, caractérisé en ce que la largeur de bande de ladite bande dans la matrice décomposée est suffisamment grande pour permettre le calcul d'une pluralité de composantes du vecteur inconnu d'ordre N indépendamment des valeurs d'une quelconque autre composante du vecteur inconnu d'ordre N.

7. Procédé selon la revendication 1, caractérisé en ce que les étapes de résolution comprennent: calculer une première pluralité de composantes inconnues du vecteur X de la variable de champ lors de cycles successifs dudit processeur en divisant des valeurs connues du vecteur Q par des éléments diagonaux correspondant de la matrice décomposée dans ledit premier processeur; et déplacer successivement les valeurs calculées du vecteur X dans chacune des autres cellules de traitement ne contenant pas les éléments diagonaux pour multiplication par d'autres éléments associés de la matrice mémorisée pour produire des produits partiels requis pour le calcul d'une seconde pluralité des composantes inconnues du vecteur X.

8. Procédé d'analyse à éléments finis pour analyser un système physique comprenant les étapes consistant à diviser la région de solution en un nombre important défini d'éléments, à exprimer les variables de champ inconnues en termes de fonctions d'approximation dans chaque élément, les fonctions étant définies en termes de valeurs des variables de champ inconnues à des noeuds spécifiés se trouvant aux frontières des éléments, à assembler les fonctions des éléments dans un premier système d'équations décrivant le comportement de la région de solution, ce premier système d'équations prenant la forme matricielle générale

KY = R

o Y représente les valeurs inconnues de la variable de champ aux noeuds du système, R représente des valeurs frontières connues, et K étant une matrice creuse de grande dimension, caractérisé en ce qu'il comprend en outre l'étape consistant à résoudre le système d'équations en transformant le système d'équations en un second système d'équations caractérisé par

une matrice triangulaire creuse de grande dimension et décom-

posée ne comprenant que des données non nulles sur sa diagonale principale, un groupement de données non nulles étant situé dans le coin supérieur droit de la matrice décomposée, et une bande de données toutes nulles étant située entre la diagonale

et ledite groupement.

9. Procédé selon la revendication 7, caractérisé en ce qu'il comprend en outre l'étape consistant à résoudre le second système d'équations en utilisant une rétro-substitution

ou une élimination directe pour déterminer les valeurs incon-

nues de la variable de champ.

10. Procédé selon la revendication 8, caractérisé en ce que la matrice est mémorisée de façon ordonnée dans une pluralité de processeurs parallèles; et une pluralité de

valeurs inconnues de la variable de champ sont calculées indé-

pendamment de toute autre valeur précalculée de la variable de champ.

11. Procédé selon la revendication 7, caractérisé en ce que la matrice décomposée est produite en renumérotant les noeuds de façon à assurer que des noeuds adjacents diffèrent d'une valeur minimale présélectionnée suffisante pour créer

ladite bande.