WO2012017177A2 - Dispositif informatique de calcul polyvalent - Google Patents

Dispositif informatique de calcul polyvalent Download PDF

Info

Publication number
WO2012017177A2
WO2012017177A2 PCT/FR2011/051859 FR2011051859W WO2012017177A2 WO 2012017177 A2 WO2012017177 A2 WO 2012017177A2 FR 2011051859 W FR2011051859 W FR 2011051859W WO 2012017177 A2 WO2012017177 A2 WO 2012017177A2
Authority
WO
WIPO (PCT)
Prior art keywords
matrix
representation
block
matrix representation
initial
Prior art date
Application number
PCT/FR2011/051859
Other languages
English (en)
Other versions
WO2012017177A3 (fr
Inventor
Laura Grigori
Frédéric NATAF
Original Assignee
Inria Institut National De Recherche En Informatique Et En Automatique
Centre National De La Recherche Scientifique (C.N.R.S)
Université Pierre Et Marie Curie (Paris 6)
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Inria Institut National De Recherche En Informatique Et En Automatique, Centre National De La Recherche Scientifique (C.N.R.S), Université Pierre Et Marie Curie (Paris 6) filed Critical Inria Institut National De Recherche En Informatique Et En Automatique
Priority to US13/813,837 priority Critical patent/US20130226980A1/en
Publication of WO2012017177A2 publication Critical patent/WO2012017177A2/fr
Publication of WO2012017177A3 publication Critical patent/WO2012017177A3/fr

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations

Definitions

  • the invention relates to the modeling and simulation of complex physical systems
  • the differential equations are solved numerically, that is, by discretizing the general equations, according to the particular parameters of the simulation.
  • These discrete systems are solved by the use of matrices of very large size, in which the discretized equations form the bases of the systems. But these basic matrices are difficult to reverse.
  • GMRES GMRES
  • preconditionnetirs Krylov subspaces
  • Preconditioners are interesting, but pose problems with some vectors for which they do not faithfully reproduce the basic matrix. To answer this problem, preconditioners "satisfying a filtering property" have been developed, which have the distinction of being faithful to the basic matrix for a particular chosen vector.
  • the invention improves the situation.
  • the invention proposes a versatile computer computing device of the type comprising:
  • a calculator-solver arranged to receive a work matrix representation and an initial matrix representation corresponding to a system of equations, as well as residue data, and to provide a solution of the system of equations from the residue data;
  • an adapter arranged to receive an initial matrix representation corresponding to a system of equations to be processed, as well as a filtering matrix representation forming a filtering matrix representation for this system of equations, and arranged to calculate a matrix representation of work corresponding to a system of equations solubie by the calculator-solver,
  • the work matrix representation being constrained to check with the initial matrix representation a stability condition comprising a comparison expression of two matrix products both having said filtering matrix representation or its transpose, and comprising respectively the initial matrix representation, and the matrix representation working,
  • the adapter is arranged to iteratively calculate in blocks an intermediate matrix from the initial matrix representation and said filter matrix representation, while the calculator-solver is arranged to work on this intermediate matrix, in blocks, so as to provide a solution of the system of equations of the initial matrix representation, without complete inversion of it, while said iterative calculation of the adapter obeys a calculation rule where a current block (i, j) of the intermediate matrix is defined by the difference between the corresponding block (ij) of the initial matrix representation and a sum of blocks each defined by a product involving two already calculated blocks of the intermediate matrix, and an auxiliary block of an approach matrix which is constrained to check with an already calculated diagonal block of the intermediate matrix a condition of equivalence, including t a comparison expression of two matrix products both having said filter matrix representation or its transposed and a previously calculated block of the intermediate matrix, and respectively comprising the inverse of said already calculated diagonal block of the intermediate matrix, and said auxiliary block of the approach matrix.
  • Such a device is particularly advantageous because it makes it
  • the invention also relates to a method comprising:
  • Operation b) comprises iteratively calculating in blocks an intermediate matrix from the initial matrix representation and said filter matrix representation by iteratively repeating, for each current index block (i, j) of the intermediate matrix:
  • Operation c) comprises working on the intermediate matrix, in blocks, so as to provide a solution of the system of equations of the initial matrix representation, without completely reversing it,
  • FIG. 1 represents a schematic view of a modeling and simulation system according to the invention
  • FIG. 2 represents a simplified flow diagram of a modeling and simulation operation by means of the system of FIG. 1,
  • FIG. 3 represents a simplified flow diagram of an operation of FIG. 2,
  • FIG. 4 represents an example of a function of calculating a preconditioner according to an operation of FIG. 3,
  • FIG. 5 represents an exemplary matrix obtained after a first reordering operation according to an optional operation of FIG. 3,
  • FIG. 6 represents an exemplary matrix obtained after a second reordering operation according to an optional operation of FIG. 3, and
  • FIG. 7 represents an example of a modification of the function of FIG. 4 to set up a parallelization of the calculations from a reordered matrix of the type of that represented in FIGS. 5 and 6.
  • Appendix A forms the formulation of certain mathematical formulas implemented in the context of the invention.
  • This Annex is set aside for the purpose of clarification, and for ease of reference, It is an integral part of the description, and may therefore serve not only to make the present invention better understood, but also to contribute to its definition, as the case may be .
  • the modeling and simulation of physical systems have become major issues. For example, in the operation of a hydrocarbon well, there is a first phase during which the oil comes out naturally. Then, as the pressure drops, it becomes necessary to act to recover the oil.
  • Preconditioners are matrices which allow to quickly approach the inverse matrix of A. By using the preconditioner M, the iterative method solves the linear system
  • Ml A x M ⁇ l b.
  • Ml A x M ⁇ l b.
  • Preconditioners satisfying a filtering property have the additional advantage of behaving identically to matrix A for a chosen vector, as is explained in formula (20) of Annex A, in which M is the preconditioner, and t the chosen vector.
  • FIG. 1 represents a polyvalent computing computing device 2 according to the invention.
  • the device 2 comprises a set of sensors 4, a digitizer 6, a discretizer 8, an adapter 10, a calculator-quiter 12, and a driver 14 which controls them.
  • the set of sensors 4 is used to obtain the data that constrains the physical system to model, and the digitizer 6 is used to transform these analog data to inject them into theoretical equations.
  • the discretizer 8 is called by the driver 14 to discretize the particularized theoretical equations with the real data, and for en. draw a system of linear equations. This system generally has a very large size, and its lines form the matrix A. Again, this element can be realized in many different ways.
  • the adapter 10 and the computer 12 are called by the driver 14 to calculate the preconditioner and to draw a solution corresponding to a particular situation that is to be modeled.
  • the "second member" data of the equation involving the matrix A are also called residue data, with reference to Newton's methods.
  • the driver 14 can call the adapter 10 and the calculator 12 to evolve. this particular solution in successive time steps, and thus to give simulation line of the evolution of the modeled physical system.
  • the adapter 10 is called once by the driver 14 to calculate a preconditioner which is used for the duration of the simulation, and the result calculated by the calculator 12 for a given time step is used. as input at the next time step.
  • the adapter 10 may be selectively called by the driver 14, depending on the evolution of the simulation, especially if it tends to modify the system at the origin of the matrix A. Again, simulation techniques are varied.
  • Figure 2 shows a simplified flow diagram showing the operations summarized above:
  • the set of sensors 4 is called to measure all the parameters necessary for the simulation
  • the digitizer 6 and the discretizer 8 are called to model the system digitally, with the measurements taken from the operation 200,
  • the adapter 10 and the computer 12 are called to perform the simulation as such.
  • Figure 3 shows a simplified flow diagram of operation 240.
  • the driver 14 transmits the matrix pulled from the operation 220 to the adapter 10.
  • the adapter 10 reorders the elements of the matrix A to allow further processing in parallel. This operation can be performed in several ways, for example by a nested dissection or by a partition into several independent domains which can also overlap and have a recursive subpartition.
  • the adapter 10 thus makes it possible to obtain a reordered matrix B which comprises null blocks.
  • the adapter 10 processes the matrix A, reordered or not, for in. derive a representation of a preconditioner M satisfying a filtering property.
  • the pilot 14 calls the computer 12 with the representation of the preconditioner M to perform the simulation.
  • the device formed by the adapter 10, the computer 12 and the driver 14, therefore makes it possible:
  • FIG. 4 represents a calculation flow diagram of the preconditioner according to the invention. This calculation is based on the decomposition of the preconditioner M in the form of the formula (30) of Appendix A.
  • This decomposition in LDU matrix is known in principle, but the calculation of the elements is different.
  • Formulas (40), (50) and (60) in Appendix A give the composition of these respective elements.
  • the decomposition of the preconditioner in the form of LDU is very advantageous because it makes it possible to solve the real system without having to invert the matrix M. More precisely, the technique contains numerous algorithms which allow a simplified resolution of a matrix equation when the decomposition LDU is used.
  • preconditioner M As such has never been done, and only its LDU components are calculated and stored.
  • solver-calder 12 selectively calls them to solve the system. It would nonetheless be possible to calculate the preconditioner M, by applying formula (30).
  • each element of the matrix C corresponds to a single element of each of the matrices L, D, U.
  • matrix C can be established according to formula (80) of Annex A, in which the term F kj satisfies formula (90) of Annex A.
  • the formula (90) expresses the fact that the matrix F which groups the terms is a matrix which approaches the diagonal block D kk for the condition relating to the index j. Thus, the fact that the formula (90) satisfies can be seen as a condition of equivalence with the inverse of the diagonal block D kk for the condition relating to the index j.
  • FIG. 4 represents an example of a function making it possible to calculate the matrices L, D and U. For this, each term of the matrix is calculated, and assigned to the matrix L, D and U as appropriate.
  • the terms are not assigned to each matrix L, D, U, but to the only matrix C which is directly used by the solver-calculator 12.
  • the matrix C is equivalent to the matrices L, D, and U, and these two variants represent only different ways of expressing the preconditioner M.
  • the first element of the matrix D can be calculated directly, since it corresponds to the first term of the diagonal of the matrix A (or B if the operation 320 is executed), likewise all the terms non-zero respectively of the first column of the matrix L and the first row of the matrix U are initialized with the corresponding term of 3a matrix A, an index i is semialized to 1. This is done in an operation 400.
  • first element or first term block. Indeed, the matrix A (or B if the operation 320 is executed) is precut in rectangular blocks whose sides are parameters that can be chosen freely. Only the diagonal blocks of A must be square. To date, the Applicants use side values such that the product of these values is equal to the size of the buffer, that is, a given block of the matrix A can be stored in the buffer memory.
  • the Applicants also use side values such that their product is smaller than the size of the buffer. If the operation 320 is executed, the size of the blocks of the matrix B is determined by this operation.
  • the global loop consists, at each iteration, in calculating the diagonal term first, then calculating by means of a local loop the other terms, by increasing index of row and column.
  • Operation 404 includes calculating the matrix F such that formula (90) is satisfied according to formula (100). If the second calculation mode is retained, the calculation of G is also performed, according to formula (110).
  • the index j of the local loop is initialized to 2 in an operation 406. This is followed by a local loop end operation 408 which is tested if j is equal to i. When i is 2, it allows to go directly to the next global loop iteration, as Cu and C 21 are known,
  • the local loop is then executed, with the calculations of Lj j which corresponds to Cy in an operation 410, and U ji which corresponds to C i in an operation 412.
  • the index of the local loop] is incremented in an operation 414, and the local loop resumes with the test of the operation 408.
  • a test checks in an operation 416 if i is equal to the number of blocks N of the matrix M.
  • the global loop is complete, and the operation ends in an operation 418 called the -solver calculator (12). Otherwise, the global loop resumes with the incrementation of the index i of the global loop in the operation 402.
  • the function of FIG. 4 makes it possible to obtain a decomposition of M in the form LDU, which serves as a basis for the known resolution methods in algebra.
  • the preconditioner M is not explicitly calculated.
  • the Applicants have developed a device implementing a preconditioner and a calculation method of a representation of this preconditiormeur which are particularly suitable for the parailelement calculations. To better understand this, it is appropriate to explain in more detail an embodiment of the operation 320. This example will be based on the case of a two-level nested dissection.
  • the matrix A is modified to give it the shape of a matrix B in the form of an "arrow".
  • the matrix A is "reordered" a first time to give it the shape of the matrix shown in FIG. 5, then the sub-matrices of the matrix of FIG. 5 are themselves reordered in the same way.
  • the matrix of FIG. 5 comprises three diagonal blocks B1, B2 and B3, two blocks B4 and B5 respectively along the bottom and right edges of the matrix B, and is zero elsewhere. This reordering of the matrix A is possible because of its low density. The same reordering is performed on the vector t.
  • the block B3 of FIG. 5 is the bioc B77 of FIG. 6, and that the blocks B4 and B5 of FIG. 5 respectively correspond to the blocks B71 to B76 and to the blocks B17 to B67 of FIG. 6.
  • the formula (80) of the Annex A it appears from the first line of this formula that blocks C1 to C17 and C21 to C71 are known directly from B. It also appears that in the application of the second line of this formula many blocks are null or known, which makes it possible to calculate some blocks Q j independently of one another, and therefore that their calculation can be performed in parallel.
  • a task dependency graph in which for a given level, all the nodes represent blocks that can be computed in parallel for the application of the formula (80) of the Annex. A. Once all the blocks linked to the nodes of a given level are calculated, it becomes possible to calculate the blocks linked to the nodes of the next level in the graph, again in parallel.
  • the task dependency graph is calculated by traversing the structure of the matrix B. The nodes of this graph represent tasks, that is, block calculations Cjj of the formula (80) of Appendix A.
  • a dynamic scheduler assigns during parallel execution tasks that are ready to run on available processors.
  • a static scheduling establishes in a first phase the order of the tasks to be performed on the different parallel computing units in order to minimize the parallel computing time.
  • a static data distribution can be used on sub-tree to sub-cube or bi-dimensional processors.
  • index i which corresponds to the level of the task dependency graph currently being traversed, is incremented in an operation 710.
  • the driver 4 calls the adapter 10 in an operation 720 to retrieve all the nodes of the level i of the task dependency graph, by means of a function Dep Gr ().
  • the result of this function is stored in a list List, which is a local variable that contains at each iteration the list of pairs (ki) that identify independent blocks at the same level of the task dependency graph.
  • the adapter calculates the CM blocks for all the pairs of the ISL list in an operation 730. This calculation is performed in parallel, as all the blocks are independent from the formula (80) of Appendix A, and are distributed across all available processors and processor cores.
  • the function of Figure 4 operates by first calculating the diagonal term, then the terms of the line and the corresponding column, here it is the task dependency graph which determines which blocks are calculated. It should be noted that in the present example, the nature of the matrix B creates a certain symmetry of the indices of each level in the task dependency graph. However, the application of a method other than nested dissection may limit this symmetry, and the blocks can be calculated in a seemingly arbitrary order.
  • the adapter 10 checks whether the index i is less than N, the number of levels of the graph. If this is the case, then the function resumes in 710 with the incrementation of i for the next level of the task dependency graph. Otherwise, the function ends in the operation 750 with the call of the caliper-solver 12, as with the operation 418.
  • the filtering condition is expressed by the formula (20) of Appendix A.
  • This formula is a mathematical expression that the initial matrix A and the preconditioner M satisfy a stability condition which is based on the comparing their product with a vector.
  • mini-blocks are mentioned are that they must not be separated during the optional operation 320, and when the matrix A is cut into blocks. A given mini-block must always be contained in a single block of matrix A or matrix B.
  • the formula (100) of the Annex A is adapted to a matrix t with only one column, that is to say a vector, and the mini-blocks are thus of size once a, it is to say scalars.
  • the Applicants thus generalized the formula (.100) in the form of the formula (180), in which Diag () designates a function which creates a diagonal matrix whose elements are designated as arguments of this function, and in which the surgery "/.” refers to the term division of the dies.
  • A1 / .A2 is a matrix A3 where each term A3 (ij) is equal to the quotient of ⁇ (i, j) by A2 (i, j).
  • the driver 14 can be integrated with the adapter 10 and the computer 12, that is to say that they are arranged to interact, instead of being separate elements ordered who ignore each other.
  • the presentation of the elements of system 2 is mainly functional.
  • these elements can be separated physically and connected by communication links, or implemented in a distant way in time, or put in place. implemented on the same equipment with the driver 14 defined by the intrinsic links between these elements and a user interface.
  • the diseretizer 8, the adapter 10, the computer 12 and the driver 14 can be implemented in the form of analog elements, such as integrated circuits or daughterboards, or in the form of digital elements, that is, in the form of programs implemented by a computer, possibly remote and / or distributed.
  • analog elements such as integrated circuits or daughterboards
  • digital elements that is, in the form of programs implemented by a computer, possibly remote and / or distributed.
  • Matrix or matrix representation therefore means any digital data structure that allows the matrix to be processed within the scope of the invention.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

Un dispositif informatique de calcul polyvalent du type comprenant: un calculateur- sol veux (12), agencé pour recevoir une représentation matricielle de travail correspondant à un système d'équations, ainsi que des données de résidus, et pour fournir une solution du système d'équations à partir des données de résidus, un adaptateur (10), agencé pour recevoir une représentation matricielle initiale correspondant à un système d'équations à traiter, ainsi qu'une représentation matricielle filtrante pour ce système d'équations, et agencé pour calculer une représentation matricielle de travail correspondant à un système d'équations soluble par le calculateur- solveur, la représentation matricielle de travail étant contrainte à vérifier avec la représentation matricielle initiale une condition de stabilité comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée, et comportant respectivement la représentation matricielle initiale, et la représentation matricielle, de travail, l'adaptateur est agencé pour calculer itérativement par blocs une matrice intermédiaire à partir de la représentation matricielle initiale et de ladite représentation numérique de représentation matricielle filtrante, tandis que le calculateur- solveur est agencé pour travailler sur cette matrice intermédiaire, par blocs, de façon à fournir une solution du système d'équations de la représentation matricielle initiale, sans inversion complète de celle-ci.

Description

Dispositif informatique de calcul polyvalent
L'invention concerne la modélisation et la simulation des systèmes physiques complexes,
Dans de nombreux domaines de la physique moderne, les équations qui régissent un phénomène physique ne peuvent pas être résolues de manière théorique. C'est notamment le cas de tous les problèmes qui ont hait à la mécanique des fluides, par exemple dans la modélisation de l'exploitation d'un champ pétrolifère.
Dans ces situations, les équations différentielles sont résolues de manière numérique, c'est-à-dire par discrétisation des équations générales, en fonction des paramètres particuliers de la simulation. Ces systèmes discrets sont résolus par l'utilisation de matrices de très grande taille, dans lesquelles les équations discrétisées forment les bases des systèmes. Mais ces matrices de base sont difficilement inversibles.
Pour résoudre ce problème, des méthodes itératives sont largement utilisées aujourd'hui. Ces méthodes, par exemple celle dite « GMRES », sont basées sur des sous- espaces de Krylov. Dans le but d'accélérer la convergence des méthodes itératives, des "préconditionnetirs" ont été créés. Ce sont des éléments qui calculent une matrice proche de la matrice de base, et dont l'inverse peut être appliqué de manière efficace à un vecteur arbitraire.
Les préconditionneurs sont intéressants, mais posent des problèmes avec certains vecteurs pour lesquels ils ne reproduisent pas fidèlement la matrice de base. Pour répondre à ce problème, des préconditionneurs "satisfaisant une propriété de filtrage" ont été développés, qui ont la particularité d'être fidèles à la matrice de base pour un vecteur choisi particulier.
À ce jour, les méthodes pour produire des préconditionneurs satisfaisant une condition de filtrage nécessitent des matrices de base présentant une forme très spécifique, ce qui limite fortement l'usage et l'utilité des préconditionneurs satisfaisant une propriété de filtrage.
L'invention vient améliorer la situation.
A cet effet, l'invention propose un dispositif informatique de calcul polyvalent du type comprenant :
~ un calculateur-solveur, agencé pour recevoir une représentation matricielle de travail et une représentation matricielle initiale correspondant à un système d'équations, ainsi que des données de résidus, et pour fournir une solution du système d'équations à partir des données de résidus,
- un adaptateur, agencé pour recevoir une représentation matricielle initiale correspondant à un système d'équations à traiter, ainsi qu'une représentation matricielle filtrante formant représentation matricielle filtrante pour ce système d'équations, et agencé pour calculer une représentation matricielle de travail correspondant à un système d'équations solubie par le calculateur-solveur,
îa représentation matricielle de travail étant contrainte à vérifier avec la représentation matricielle initiale une condition de stabilité comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée, et comportant respectivement la représentation matricielle initiale, et la représentation matricielle de travail,
L'adaptateur est agencé pour calculer itérativement par blocs une matrice intermédiaire à partir de la représentation matricielle initiale et de ladite représentation matricielle filtrante, tandis que le calculateur-solveur est agencé pour travailler sur cette matrice intermédiaire, par blocs, de façon à fournir une solution du système d'équations de la représentation matricielle initiale, sans inversion complète de celle-ci, tandis que ledit calcul itératif de l'adaptateur obéit à une règle de calcul où un bloc courant (i,j) de la matrice intermédiaire est défini par la différence entre le bloc (i.j) correspondant de la représentation matricielle initiale et une somme de blocs chacun défini par un produit faisant intervenir deux blocs déjà calculés de la matrice intermédiaire, et un bloc auxiliaire d'une matrice d'approche qui est contraint à vérifier avec un bloc diagonal déjà calculé de la matrice intermédiaire une condition d'équivalence, comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée et un bloc déjà calculé de la matrice intermédiaire, et comportant respectivement l'inverse dudit bloc diagonal déjà calculé de la matrice intermédiaire, et ledit bloc auxiliaire de la matrice d'approche. Un tel dispositif est particulièrement avantageux, car il permet de mettre en œuvre un préconditionneur satisfaisant une condition de filtrage, et ce quelque soit la matrice de base qui définit le système que l'on cherche à résoudre.
L'invention concerne également un procédé comprenant :
a) recevoir une représentation matricielle initiale correspondant à un système d'équations à traiter et une représentation matricielle filtrante,
b) calculer une représentation matricielle de travail vérifiant avec la représentation matricielle initiale une condition de stabilité comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée, et comportant respectivement la représentation matricielle initiale, et la représentation .matricielle de travail,
c) recevoir des données de résidus, et résoudre le système d'équations défini par la représentation matricielle initiale, à partir des données de résidus, de la représentation matricielle de travail et de la représentation matricielle initiale.
L'opération b) comprend calculer itérativement par blocs une matrice intermédiaire à partir de ia représentation matricielle initiale et de ladite représentation matricielle filtrante en répétant itérativement, pour chaque bloc d'indice courant (i,j) de la matrice intermédiaire :
bl) calculer une somme de blocs chacun défini par un produit faisant intervenir deux blocs déjà calculés de la matrice intermédiaire, et un bloc auxiliaire d'une matrice d'approche qui est contraint à vérifier avec un bloc diagonal déjà calculé de la matrice intermédiaire une condition d'équivalence, comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée et un bloc déjà calculé de ia matrice intermédiaire, et comportant respectivement l'inverse dudit bloc diagonal déjà calculé de la matrice intermédiaire, et ledit bloc auxiliaire de la matrice d'approche
b!2) calculer la différence entre le bloc (i.j) de la matrice de la représentation matricielle initiale et la somme issue de l'opération bl). L'opération c) comprend travailler sur la matrice intermédiaire, par blocs, de façon à fournir une solution du système d'équations de la représentation matricielle initiale, sans inversion complète de celle-ci,
D'autres caractéristiques et avantages de l'invention apparaîtront mieux à la lecture de la description qui suit, tirée d'exemples donnés à titre iliustratif et non limitatif, tirés des dessins sur lesquels :
- la figure 1 représente une vue schématique d'un système de modélisation et de simulation selon l'invention,
- la figure 2 représente un diagramme de flux simplifié d'une opération de modélisation et de simulation au moyen du système de la figure 1,
- la figure 3 représente un diagramme de flux simplifié d'une opération de la figure 2,
~ la figure 4 représente un exemple d'une fonction de calcul d'un préconditionneur selon une opération de îa figure 3,
- la figure 5 représente un exemple de matrice obtenue après une première opération de réordonnancement selon une opération optionnelle de la figure 3,
- la figure 6 représente un exemple de matrice obtenue après une deuxième opération de réordonnaneement selon une opération optionnelle de la figure 3, et
- la figure 7 représente un exemple de modification de la fonction de la figure 4 pour mettre en place une parallélisation des calculs à partir d'une matrice réordonnée du type de celle représentée sur les figures 5 et 6.
Les dessins et la description ci-après contiennent, pour l'essentiel, des éléments de caractère certain. îls pourront donc non seulement servir à mieux faire comprendre la présente invention, mais aussi contribuer à sa définition, le cas échéant.
En outre, la description détaillée est augmentée de l'annexe A, qui dorme la formulation de certaines formules matliématiques mises en œuvre dans le cadre de l'invention. Cette Annexe est mise à part dans ira but de clarification, et pour faciliter les renvois, Elle est partie intégrante de la description, et pourra donc non seulement servir à mieux faire comprendre la présente invention, mais aussi contribuer à sa définition, le cas échéant. La modélisation et la simulation des systèmes physiques sont devenues des enjeux de taille. Par exemple, dans l'exploitation d'un forage d'hydrocarbure, il y a une première phase pendant laquelle le pétrole sort naturellement. Puis, au fur et à mesure que la pression baisse, il devient nécessaire d'agir pour récupérer le pétrole.
Pour cela, il est possible par exemple d'utiliser un flux d'eau, qui est introduit dans le puits pour faire remonter la pression et faire rejaillir le pétrole. Mais ces opérations périlleuses nécessitent une connaissance poussée du puits et des réactions de celui-ci dans ces circonstances.
Les équations qui déterminent ce problème physique sont très complexes, et pour la plupart n'admettent que des solutions par discrétisation et méthode numérique de type différences finies ou volumes finis. Les problèmes ainsi discrétisés peuvent alors être résumés à la formule (10) de l'Annexe A, dans laquelle A est la matrice de base qui définit le système d'équations discrétisé, x est le vecteur qui est recherché, et y est le vecteur- résultat connu.
Ce type de problème est bien connu en algèbre, et il s'agit de trouver la matrice inverse de A pour calculer x. Mais l'inversion des matrices est un problème complexe, qui monopolise des puissances de calcul qui croissent de manière exponentielle avec la taille de la matrice à inverser.
Pour cela, des méthodes itératives basées sur des sous-espaces de Krylov, comme GMRES, sont largement utilisées aujourd'hui. Pour accélérer la convergence de ces méthodes, des "préconditionneurs" ont été proposés. Les préconditionneurs sont des matrices qui permettent d'approcher rapidement la matrice inverse de A. En utilisant le préconditionneur M, la méthode itérative résoud le système linéaire
M-l A x = M~l b. Dans ce mode de résolution, on calcule des opérations de type M-l v et A v, où v est un vecteur, sans calculer de manière explicite l'inverse de M. Comme cela a été expliqué plus haut, il existe une classe particulière de préconditionneurs, les préconditionneurs satisfaisant une propriété de filtrage.
Les préconditionneurs satisfaisant une propriété de filtrage présentent l'avantage supplémentaire de se comporter de manière identique à ia matrice A pour un vecteur choisi, comme cela est explicité dans la formule (20) de l'annexe A, dans laquelle M est le préconditionneur, et t le vecteur choisi.
A ce jour, seules les matrices A qui sont tridiagonales par blocs sont utilisées pour produire un préconditionneur satisfaisant une propriété de filtrage. Cela restreint considérablement leur champ d'application.
En outre, les méthodes de calcul de ces préconditionneivrs sont pour ia plupart séquentielles, ce qui les rend rapidement prohibitives en termes de coût de calcul, et donc peu utilisables en pratiques. En effet, seules les matrices issues d'un maillage structuré peuvent être traitées en parallèle, ce qui restreint considérablement leur champ d'application.
La figure 1 représente un dispositif informatique de calcul polyvalent 2 selon l'invention. Le dispositif 2 comprend un ensemble de capteurs 4, un numériseur 6, un discrétiseur 8, un adaptateur 10, un calculateur-soiveur 12, et un pilote 14 qui les commande.
Dans l'exemple décrit ici, l'ensemble de capteurs 4 servent à obtenir les données qui contraignent le système physique à modéliser, et le numériseur 6 sert à transformer ces données analogiques pour les injecter dans des équations théoriques.
Ces éléments sont pour ainsi dire indifférents au problème résolu par l'invention : ils servent à définir le cadre pour son application pratique. Aussi, leur réalisation pourra être très variée.
Le discrétiseur 8 est appelé par le pilote 14 pour discrétiser les équations théoriques particularisées avec les données réelles, et pour en. tirer un système d'équations linéaires. Ce système présente en général une taille très importante, et ses lignes forment la matrice A. Là encore, cet élément peut être réalisé de nombreuses manières différentes.
Enfin, l'adaptateur 10 et le calculateur 12 sont appelés par le pilote 14 pour calculer le préconditionneur et pour tirer une solution correspondant à une situation particulière que l'on cherche à modéliser. Dans ce cas, les données de "second membre" de l'équation faisant intervenir la matrice A sont également appelées données de résidus, en référence aux méthodes de Newton, Le pilote 14 peut appeler l'adaptateur 10 et le calculateur 12 pour faire évoluer cette solution particulière par pas de temps successifs, et pour donner ainsi line simulation de l'évolution du système physique modélisé.
Selon une première variante, l'adaptateur 10 est appelé une unique fois par le pilote 14 pour' calculer un préconditionneur qui est utilisé pour toute la durée de la simulation, et le résultat calculé par le calculateur 12 pour un pas de temps donné est utilisé comme entrée au pas de temps suivant.
Selon une deuxième variante, l'adaptateur 10 peut être sélectivement appelé par le pilote 14, en fonction de l'évolution de la simulation, notamment si celle-ci tend à modifier le système à l'origine de la matrice A. ïci encore, les techniques de simulation sont variées. La figure 2 montre un diagramme de flux simplifié montrant les opérations résumées plus haut:
- dans une opération 200, l'ensemble de capteurs 4 est appelé pour mesurer tous les paramètres nécessaires à la simulation,
- dans une opération 220, le numériseur 6 et le discrétiseur 8 sont appeiés pour modéliser le système de manière numérique, avec les mesures tirées de l'opération 200,
- dans une opération 240, l'adaptateur 10 et le calculateur 12 sont appelés pour effectuer la simulation en tant que telle.
La figure 3 montre un diagramme de flux simplifié de l'opération 240. Dans une opération 300, le pilote 14 transmet la matrice À tirée de l'opération 220 à l'adaptateur 10. Dans une opération 320 optionnelle, l'adaptateur 10 réordonne les éléments de la matrice A pour permettre un traitement ultérieur en parallèle. Cette opération peut être réalisée de plusieurs manières, par exemple par une dissection emboîtée ou par une partition en plusieurs domaines indépendants qui peuvent également se recouvrir et avoir une sous-partition récursive.
Un tel recouvrement augmente légèrement les coûts de calcul mais offre de meilleurs taux de convergence et une robustesse supérieure, comme cela est réalisé dans la méthode de Schwarz. L'adaptateur 10 permet ainsi d'obtenir une matrice B réordonnée qui comprend des blocs nuls.
Ensuite, dans une opération 340, l'adaptateur 10 traite la matrice A, réordonnée ou pas, pour en. tirer une représentation d'une préconditionneur M satisfaisant une propriété de filtrage. Enfin, dans une opération 360, îe pilote 14 appelle le calculateur 12 avec ia représentation du préconditionneur M pour réaliser la simulation.
Le dispositif formé par l'adaptateur 10, le calculateur 12 et le pilote 14, permet donc :
- de calculer une représentation d'un préconditionneur vérifiant une propriété de filtrage pour une quelconque matrice A en entrée, et
- paralléîiser de manière massive les calculs liés au préconditionneur lorsque l'adaptateur 10 est appelé pour réordonner la matrice À.
La figure 4 représente un diagramme de flux de calcul du précondiîionneur selon l'invention. Ce calcul repose sur ia décomposition du préconditionneur M sous la forme de la formule (30) de l'annexe A. Cette décomposition en matrice LDU est connue en principe, mais le calcul des éléments est différent. Les formules (40), (50) et (60) de l'annexe A donnent la composition de ces éléments respectifs. La décomposition du préconditionneur sous la forme LDU est très avantageuse, car elle permet de résoudre le système réel sans avoir à inverser la matrice M. Plus précisément, la technique contient de nombreux algorithmes qui permettent une résolution simplifiée d'une équation matricielle lorsque La décomposition LDU est utilisée.
Pour cette raison, le calcul du préconditionneur M en tant que tel n'est j amais réalisé, et seules ses composantes LDU sont calculées et stockées. Ensuite, le caleulateur- solveur 12 les appelle sélectivement pour résoudre le système. II serait néanmoins possible de calculer le préconditionneur M, en appliquant ia formule (30).
Pour calculer les éléments de la décomposition LDU, les Demandeurs ont découvert une formulation basée sur le calcul d'une matrice C qui correspond à la somme des matrices L, D et U (formule (70) de l'Annexe A).
Du fait de la forme respective des matrices L, D et U, il apparaît que chaque élément de la matrice C correspond à un unique élément de chacune des matrices L, D, U. Ainsi :
- Dii = Cii, et Djj - 0 pour i différent de j,
- Lij = Cjj pour i>1 et i strictement supérieur à j, et Lij = 0 pour i inférieur ou égal à j ,
- Uij - Oj pour j>1 et j strictement supérieur à i, et Uij = 0 pour j inférieur ou égal à i. Les Demandeurs ont découvert que la matrice C peut être établie conformément à la formule (80) de l'annexe A, dans laquelle le terme Fkj satisfait la formule (90) de l'annexe A.
Dans la formule (80), la première ligne ne représente qu'une initialisation. Conccptuellernent, cette première ligne est équivalente à la deuxième ligne. En effet, pour i=1 ou j=1, alors rnin(ij)-1 vaut zéro, ce qui signifie que la somme de cette deuxième ligne ne comprend aucun terme, et donne un résultat identique à la première ligne. La formule (90) exprime le fait que la matrice F qui regroupe les termes est une matrice qui approche le bloc diagonal Dkk pour la condition portant sur l'indice j. Ainsi, le fait que satisfait la formule (90) peut être vu comme une condition d'équivalence avec l'inverse du bloc diagonal Dkk pour la condition portant sur l'indice j .
Lorsque le vecteur Ukjtj n'a aucune composante nulle, le calcul de Fy est aisé. Cependant, pour tenir compte des cas où ce vecteur présente des composantes nulles, il est possible de modifier les approches de ia littérature conformément à la formule (100) de l'annexe A,
Les Demandeurs ont également découvert un autre calcul pour la matrice C, dans laquelle, dans la formule (80), le terme Fkj est remplacé par le terme Gy, qui est défini avec la formule (110) de l'annexe A. Pour calculer le terme Gkj, on calcule d'abord le terme Fy correspondant grâce à la formule (100), puis on applique la formule (110). ïl est aussi possible de définir le terme Fkj satisfaisant la condition de filtrage (90) par la formule (120) combinée à la formule (140) de l'Annexe A, qui est issue des méthodes de déflation de l'algèbre linéaire. Dans le cas où les blocs sont symétriques, on peut utiliser la formule (130) de l'Annexe A, qui est une version simplifiée de la formule (120).
Dans l'état actuel des recherches des Demandeurs, la formule (100) est préférée à la formule (1 10), pour des raisons de stabilité, La figure 4 représente un exemple de fonction permettant de calculer les matrices L, D et U. Pour cela, chaque terme de la matrice est calculé, et affecté à la matrice L, D et U comme il convient.
En variante, les termes ne sont, pas affectés à chaque matrice L, D, U, mais à la seule matrice C. qui est directement utilisée par le caleulateur-solveur 12, En fait, on a vu plus haut que la matrice C est équivalente aux matrices L, D, et U, et ces deux variantes représentent seulement des manières différentes d'exprimer le précondîtionneur M. Au vu de ce qui précède, le premier élément de la matrice D peut être calculé directement, puisqu'il correspond au premier terme de la diagonale de la matrice À (ou B si l'opération 320 est exécutée), De même tous les termes non nuls de respectivement la première colonne de la matrice L et la première ligne de la matrice U sont initialisés avec le terme correspondant de 3a matrice A, Un indice i est mitialisé à 1. Cela est réalisé dans une opération 400. Par premier élément ou premier terme, on entend bloc. En effet, la matrice A (ou B si l'opération 320 est exécutée) est prédécoupée en blocs rectangulaires dont les côtés sont des paramètres qui peuvent être choisis librement. Seuls les blocs diagonaux de A doivent être carrés. À ce jour, les Demandeurs utilisent des valeurs de côtés telles que le produit de ces valeurs est égal à la taille de la mémoire tampon, c'est-à-dire qu'un bloc donné de la matrice A peut être stocké dans le tampon mémoire.
En variante, les Demandeurs utilisent également de valeurs de côtés telles que leur produit est inférieur à la taille de la mémoire tampon. Si l'opération 320 est exécutée, la taille des blocs de la matrice B est déterminée par cette opération.
Ensuite, on exécute une boucle dite globale qui va calculer tous les autres termes de la matrice C, et donc tous les termes des matrices L, D, ei U, ce qui permet de définir le préconditionneur M.
La boucle globale consiste, à chaque itération, à calculer d'abord le terme diagonal, puis à calculer au moyen d'une boucle locale les autres termes, par indice croissant de ligne et de colonne.
D'abord, l'indice i de la boucle globale est Încrémenté dans une opération 402, puis dans une opération 404? le terme diagonal DÛ est calculé conformément à la formule (80). Le ternie Dii correspond comme on l'a vu précédemment au terme C«. L'opération 404 inclut le calcul de la matrice F tel que la formule (90) soit satisfaite, conformément à la formule (100). Si le deuxième mode de calcul est retenu, le calcul de G est également réalisé, conformément à la formule (110).
Ensuite, l'indice j de la boucle locale est initiaîisé à 2 dans une opération 406, Cela est suivi d'une opération 408 de fin de boucle locale laquelle il est testé si j est égal à i. Quand i vaut 2, cela permet de passer directement à l'itération de boucle globale suivante, comme Cu et C21 sont connus,
La boucle locale est alors exécutée, avec les calculs de Ljj qui correspond à Cy dans une opération 410, et de Uji qui correspond à Cjidans une opération 412.
On notera que les opérations 410 et 4! 2 peuvent être réalisées en parallèle, ce qui est avantageux. En effet, de par la formule (80), le calcul des termes Cy et ¾ est indépendant.
Ensuite, l'indice de la boucle locale] est Încrémenté dans tme opération 414, et la boucle locale reprend avec Le test de l'opération 408. Lorsque la boucle locale est terminée, c'est-à-dire lorsque tous les termes de la ligne de L et tous les termes de la colonne de U ont été calculés, un test vérifie dans une opération 416 si i est égal au nombre de blocs N de la matrice M.
Si c'est îe cas, alors la boucle globale est terminée, et l'opération se termine dans une opération 418 l'appel du calculateur -solveur (12). Sinon, la boucle globale reprend avec l'incrémentation de l'indice i de la boucle globale dans l'opération 402.
Comme mentionné plus haut, la fonction de la figure 4 permet d'obtenir une décomposition de M sous la forme LDU, qui sert de base aux méthodes de résolutions connues en algèbre. Ainsi, le préconditionneur M n'est pas explicitement calculé. Cependant il serait possible de calculer explicitement le précondi donneur, en appliquant la formule (30). Les Demandeurs ont développé un dispositif mettant en œuvre un préconditionneur et une méthode de calcul d'une représentation de ce préconditiormeur qui sont particulièrement adaptés à la parailélisation des calculs. Pour mieux comprendre cela, il convient d'expliquer plus en détail un exemple de réalisation de l'opération 320. Cet exemple sera ici basé sur le cas d'une dissection emboîtée à deux niveaux.
Dans ce type d'opération, la matrice A est modifiée pour lui donner la forme d'une matrice B en forme de "flèche". Pour cela, la matrice A est "réordonnée" une première fois pour- lui donner la forme de la matrice représentée sur la figure 5, puis les sous- matrices de la matrice de la figure 5 sont elles-mêmes réordonnées de la même manière.
La matrice de la figure 5 comporte trois blocs diagonaux Bl, B2 et B3, deux blocs B4 et B5 respectivement le long des bords bas et droit de la matrice B, et est nulle ailleurs. Ce réordonnancement de la matrice A est possible du fait de sa faible densité. Le même réordonnancement est réalisé sur le vecteur t.
Une fois que la matrice de figure 5 a été calculée» il est possible de réappliquer ce même réordormancement aux blocs Bl et B2, ce qui aboutit à la matrice B de Sa figure 6. On remarquera sur cette figure que le bloc B3 de la figure 5 est le bioc B77 de la figure 6, et que les blocs B4 et B5 de la figure 5 correspondent respectivement aux blocs B71 à B76 et aux blocs B17 à B67 de la figure 6. Si on analyse la formule (80) de l'Annexe A, il apparaît par la première ligne de cette formule que les blocs Cl 1 à C17 et C21 à C71 sont connus directement à partir de B. 11 apparaît également que, dans l'application de la deuxième ligne de cette formule, de nombreux blocs sont nuls ou connus, ce qui rend possible de calculer certains blocs Qj de manière indépendante les uns des antres, et donc que leur calcul peut être réalisé en parallèle.
Ainsi, Bl 1, B22, B44 et B55. Si l'on calcule ces blocs en parallèle, ia même situation se répète, et de nouveaux blocs B33 et B66 peuvent à leur tour être calculés en parallèle. Et ainsi de suite. D'une manière générale, les Demandeurs ont donc découvert que l'opération 320 peut être utilisée pour produire une matrice B équivalente à la matrice A, et qui comporte des domaines séparés.
À partir de ces domaines, il est possible de créer un graphe de dépendance de tâches, dans lequel pour un niveau donné, tous les nœuds représentent des blocs pouvant être calculés en parallèle pour l'application de la formule (80) de l'Annexe A. Une fois tous les blocs liés aux nœuds d'un niveau donné calculés, il devient possible de calculer les blocs liés aux nœuds du prochain niveau dans le graphe, là encore en parallèle. Pour calculer et ordonnancer le graphe de dépendance de tâches, plusieurs techniques peuvent être utilisées, comme cela est connu dans la théorie des graphes. Le graphe de dépendances de tâches est calculé en parcourant la structure de la matrice B. Les nceuds de ce graphe représentent des tâches, c'est à dire des calculs de blocs Cjj de la formule (80) de l'Annexe A.
Les dépendances dans ce graphe représentent l'ordre des calculs imposé par la formule (80) de l'Annexe A. Ce graphe peut être ensuite ordonnancé en utilisant un ordonnancement statique ou un ordonnancement dynamique des tâches sur les processeurs.
Par exemple, un ordonnancement dynamique affecte au cours de l'exécution parallèle les tâches prêtes à être exécutées sur les processeurs disponibles. Un ordonnancement statique établit dans une première phase l'ordre des tâches à exécuter sur les différentes unités de calcul parallèles dans le but de minimiser le temps de calcul parallèle.
Dans une deuxième phase leur exécution a lieu. Une distribution statique des données peut utilisée sur les processeurs de type sous-arbre vers sous-cube ou de type bi- dimensionneL
C'est cette observation qui a amené les modifications du diagramme de la figure 4, qui sont présentées avec la figure 7. Dans cette fonction la première opération 700 correspond à l'opération 400 de la figure 4, Une différence avec ropération 400 est que l'opération 700 comprend le calcul du graphe de dépendance de tâches décrit plus haut
Ensuite l'indice i, qui correspond ait niveau du graphe de dépendance de tâches qui est actuellement parcouru, est incrémenté dans une opération 710.
Une fois l'indice i incrémenté, le pilote Î4 appelle l'adaptateur 10 dans une opération 720 pour récupérer tous les nœuds du niveau i du graphe de dépendance de tâches, au moyen d'une fonction Dep Gr().
Le résultat de cette fonction est stocké dans une liste List, qui est une variable locale qui contient à chaque itération la liste des couples (ki) qui identifient les blocs indépendants de même niveau du graphe de dépendance de tâches.
Ensuite, l'adaptateur calcule les blocs CM pour tous les couples de la iîst List dans une opération 730. Ce calcul est exécuté en parallèle, comme tous les blocs sont indépendants par rapport à la formule (80) de l'Annexe A, et sont répartis sur tous les processeurs et les cœurs de processeurs disponibles.
On remarquera ici que le calcul des blocs par l'opération 730 est différent de celui des opérations 404, 410 et 412. En effet, si la formule utilisée pour le calcul est la même, les indices des blocs sont totalement indépendants.
Là où la fonction de la figure 4 opère en calculant d'abord le terme diagonal, puis les termes de la ligne et la colonne correspondante, ici c'est le graphe de dépendance de tâches qui détermine quels sont les blocs calculés. On notera que dans l'exemple présent, la nature de la matrice B crée une certaine symétrie des indices de chaque niveau dans le graphe de dépendance de tâches. Cependant, l'application d'une autre méthode que la dissection emboîtée pourra limiter cette symétrie, et les blocs pourront être calculés dans un ordre apparemment arbitraire. Enfin dans une opération 740, l'adaptateur 10 vérifie si l'indice i est inférieur à N, le nombre de niveaux du graphe. Si c'est le cas, alors la fonction reprend en 7Î0 avec l'incrémentation de i pour le niveau suivant du graphe de dépendance de tâches. Sinon, la fonction se termine dans l'opération 750 avec l'appel du caleuîateur-solveur 12, comme avec l'opération 418.
Dans ce qui précède, la condition de filtrage est exprimée par la formule (20) de l'Annexe A. Cette formule est une expression mathématique du fait que la matrice initiale A et le préconditiormeur M vérifient une condition de stabilité qui est basée sur la comparaison de leur produit avec un vecteur.
Cependant, la condition de stabilité ne doit pas être limitée à la seule formule (20). Ainsi, les Demandeurs ont également utilisé avec succès la formule (150) de l'Annexe A.
Comme la formule (150) est presque la transposée de la formule (20), l'utilisation de la formule (150) comme condition de stabilité ne change rien au mode de fonctionnement de l'invention. En conséquence» les formules (80) et (90) n'ont qu'à être légèrement modifiées, comme présenté avec les formules (160) et (170). Les formules (100) à (1.40) pourront être adaptées de manière similaire.
En outre, tous les exemples qui précèdent ont été réalisés pour une condition de stabilité utilisant un vecteur t. Cependant, lorsqu'un système physique est modéiisé, de nombreuses grandeurs sont utilisées. Les expériences des Demandeurs montrent qu'il est avantageux d'utiliser une condition de stabilité utilisant une matrice dont chaque colonne concerne une grandeur physique. Ainsi, si deax grandeurs physiques caractérisent une mise en équation donnée, il est avantageux d'utiliser comme élément de filtrage t une matrice ayant deux colonnes. Dans la pratique, cela ne change pas la philosophie de l'invention, et les calculs présentés précédemment s'en trouvent peu ou pas modifiés. En effet, dans ce type de situation, les équations représentées dans la matrice A vont être associées par "mini-blocs" carrés dont ie côté est égal au nombre de colonne de la matrice t. Donc, dans le cas décrit au paragraphe précédent, chaque mini-bloc serait un bloc carré de deux fois deux termes de la matrice A initiale.
La raison pour laquelle ces mini-blocs sont mentionnés est qu'ils ne doivent pas être séparés îors de l'opération optionnelle 320, et lorsque la matrice A est découpée en blocs. Un mini-bloc donné doit toujours être contenu dans un unique bloc de la matrice A ou de la matrice B.
Le seul élément qui change légèrement est le calcul de !¾. En effet, la formule (100) de l'Annexe À est adaptée à une matrice t à une seule colonne, c'est-à-dire un vecteur, et les mini-blocs sont donc de taille un fois un, c'est-à-dire des scalaires. Les Demandeurs ont donc généralisé la formule (.100) sous la forme de la formule (180), dans laquelle Diag() désigne une fonction qui crée une matrice diagonale dont les éléments sont désignés en arguments de cette fonction, et dans laquelle l'opération "/." désigne la division terme à ternie des matrices. Ainsi A1/.A2 est une matrice A3 dont chaque terme A3(ij) est égal au quotient de Àî(i,j) par A2(i,j).
Une autre manière de voir cette modification est de remarquer que la formule (100) peut être vue comme un cas particulier de la formule (180), dans laquelle t a une seule colonne. Dans ce qui précède, l'adaptateur 10, le calculateur 12 et ie pilote 14 peuvent être réalisés de plusieurs manières.
Tout d'abord, le pilote 14 peut être réalisé de manière intégrée avec l'adaptateur 10 et le calculateur 12, c'est-à-dire que ceux-ci sont agencés pour savoir interagir, au lieu d'être des éléments séparés commandés qui s'ignorent.
De plus, la présentation des éléments du système 2 est principalement fonctionnelle. Ainsi, ces éléments peuvent être séparés physiquement et reliés par des liens de communication, ou mis en œuvre de manière éloignée dans le temps, ou encore mis en œuvre sur un même équipement avec le pilote 14 défini par les liaisons intrinsèques entre ces éléments et une interface utilisateur.
En outre, le diserétiseur 8, l'adaptateur 10, le calculateur 12 et le pilote 14 peuvent être mis en œuvre sous la forme d'éléments analogiques, comme des circuits intégrés ou des cartes filles, ou sous la forme d'éléments numériques, c'est-à-dire sous la forme de programmes mis en œuvre par un ordinateur, éventuellement de manière éloignée et/ou répartie. On notera également que, dans ce qui précède, il est souvent indifféremment fait référence à des matrices ou à leur représentation, il va de soi qu'un ordinateur- ne sait pas ce qu'est une matrice, et que c'est bien là représentation numérique de ces matrices, c'est-à-dire les données qui définissent ces matrices qui sont visées. Par matrice ou par représentation matricielle, on entend donc toute structure de données numériques qui pennet de traiter la matrice dans le cadre de l'invention.
On noiera enfin la visée particulièrement pratique du dispositif de l'invention, qui pennet la simulation et la résolution de nombreux problèmes physiques qui n'étaient pas semblés précédemment, par exemple dans l'industrie pétrolière.
Figure imgf000021_0001
Figure imgf000022_0001

Claims

Revendications
1. Dispositif informatique de calcul polyvalent du type comprenant :
- un calcuiateur-solveur (12), agencé pour recevoir une représentation matricielle de travail (M) et une représentation matricielle initiale (A) correspondant à un système d'équations, ainsi que des données de résidus, et pour fournir une solution du système d'équations à partir des données de résidus,
- un adaptateur (10), agencé pour recevoir une représentation matricielle initiale (A) correspondant à un système d'équations à traiter, ainsi qu'une représentation matricielle filtrante (t) pour ce système d'équations, et agencé pour calculer une représentation matricielle de travail (M) correspondant à un système d'équations soluble par le calcuiateur-solveur,
la représentation matricielle de travail (M) étant contrainte à vérifier avec la représentation matricielle initiale (A) une condition de stabilité ((20), (150)) comprenant une comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante (t) ou sa transposée, et comportant respectivement la représentation matricielle initiale (A), et la représentation matricielle de travail (M),
caractérisé en ce que l'adaptateur est agencé pour calculer itérât! vement par blocs une matrice intermédiaire (C) à partir de la représentation matricielle initiale (A) et de ladite représentation numérique de représentation matricielle filtrante (t), tandis que le calcuiateur-solveur est agencé pour travailler sur cette matrice intermédiaire, par blocs, de façon à. fournir une solution du système d'équations de Sa représentation matricielle initiale (A), sans inversion complète de celle-ci,
tandis que ledit calcul itératif de l'adaptateur obéit à une règle de calcul ((80) ; (160)) où un bloc courant (i,j) de la matrice intermédiaire (C) est défini par la différence entre le bloc (i j) correspondant de la représentation matricielle initiale (A) et une somme de blocs chacun, défini par un produit faisant intervenir deux blocs déjà calculés de la matrice intermédiaire (C), et un bloc auxiliaire d'une matrice d'approche (F) qui est contraint à vérifier avec un bloc diagonal déjà calculé de îa matrice intermédiaire une condition d'équivalence ((90), (170)), comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante (t) ou sa transposée et un bloc déjà calculé de la matrice intermédiaire, et comportant respectivement l'inverse dudit bîoc diagonal déjà calculé de la matrice intermédiaire (C), et ledit bîoc auxiliaire de la matrice d'approche (F),
2. Dispositif selon la revendication 1, dans lequel la condition de stabilité (20) comprend une comparaison de deux produits matriciels, tous deux entre respectivement la représentation matricielle initiale (A), et la représentation matricielle de travail (M), et ladite représentation matricielle filtrante (t), et dans lequel la somme de blocs est réalisée, pour un indice k non nul et strictement inférieur au minimum entre i et j, chaque bloc de cette somme étant définie par le produit matriciel de la forme PQR,
* où P est le bloc (ijk) de la matrice intermédiaire (C),
* où Q est le bloc auxiliaire (kj) de la matrice d'approche (F) d'indice, la condition d'équivalence (90) comprenant la comparaison de deux produits matriciels de respectivement ledit bloc auxiliaire (kj) de la maixice d'approche (F) et l'inverse du bloc diagonal (k,k) déjà calculé de la matrice intermédiaire (C), avec îe bloc (kj) déjà calculé de la représentation matricielle intermédiaire (C) et avec le bloc j de ladite représentation matricielle filtrante (t), et
* où R est le bloc (k j) de la matrice intermédiaire (C).
3. Dispositif selon la revendication 1 , dans lequel la condition de stabilité (150) comprend une comparaison de deux produits matriciels, tous deux entre la transposée de la représentation matricielle filtrante (t) et respectivement la représentation matricielle initiale (A), et la représentation matricielle de travail (M), et dans lequel la somme de blocs est réalisée, pour un indice k non nul et strictement inférieur au minimum entre i et j, chaque bloc de cette somme étant définie par le produit matriciel de la forme PQR,
* où P est le bloc (i,k) de la matrice intermédiaire (C),
* où Q est le bloc auxiliaire (i.k) de la matrice d'approche (F), la condition d'équivalence (170) comprenant la comparaison de deux produits matriciels du bloc i de ladite représentation matricielle filtrante (t) avec le bloc (i,k) déjà calculé de la représentation matricielle intermédiaire (C), et avec respectivement ledit bloc auxiliaire (î,k) de la matrice d'approche (F) et J'inverse du bloc diagonal (k.k) déjà calculé de la matrice intermédiaire (C), et
* où R est le bloc (kj) de la matrice intermédiaire (C),
4. Dispositif selon l'une des revendications précédentes, dans lequel le bloc auxiliaire de la matrice d'approche (F) est calculé à partir d'une division terme à terme faisant intervenir le bloc déjà calculé de la matrice 'intermédiaire (C) d'indice (k,k), et soit le bloc j de ladite représentation matricielle filtrante (t) et le bloc (kj) déjà calculé de la matrice intermédiaire (C), soit le bloc i ladite représentation matricielle filtrante (t) et le bloc (i,k) déjà calculé de la matrice intermédiaire (C).
5. Dispositif selon l'une des revendications 1 à 3, dans lequel la matrice d'approche est calculée par une méthode de déflation, dans laquelle un premier terme (Z) fait intervenir soit le bloc j de ladite représentation matricielle filtrante (t) et le bloc (k j) déjà calculé de la matrice intermédiaire (C), soit le bloc i de ladite représentation matricielle filtrante (t) et le bloc (isk) déjà calculé de la matrice intermédiaire (C), dans laquelle un second terme (Q) fait intervenir le bloc diagonal (k,k) déjà calculé de la matrice intermédiaire (C) et le premier terme, et un troisième terme (P) fait intervenir les premier et seconds termes ainsi que la matrice de la représentation matricielle initiale (A).
6. Dispositif selon l'une des revendications précédentes, dans lequel la représentation matricielle filtrante (t) est un vecteur colonne.
7. Dispositif selon Tune des revendications précédentes, dans lequel l'adaptateur (10) est agencé pour réordonner la représentation matricielle initiale (A) pour produire une représentation matricielle modifiée (B) selon une règle d'ordonnancement agencée pour associer des blocs de la matrice de la représentation matricielle initiale (A) en fonction d'une condition de dépendance ((80),(160)), et pour calculer la représentation matricielle de travail (M) en remplaçant représentation matricielle initiale (A) par la représentation matricielle modifiée (B).
8. Dispositif selon la revendication 7, dans lequel l'adaptateur (10) est agencé pour calculer un graphe à partir de la représentation matricielle modifiée (B) dans lequel pour un niveau donné les calculs des blocs de la matrice intermédiaire (C) peuvent être réalisés de manière indépendante, et pour calculer les blocs d'un même niveau du graphe en parallèle.
9. Dispositif selon l'une des revendications précédentes, comprenant en outre un ensemble de capteurs 4, un numériseur 6, un discrétiseur 8 et un pilote 14, dans lequel le pilote 14 est agencé pour appeler le discrétiseur 8 avec des données tirées du numériseur 6 qui opère sur des données tirées de l'ensemble de capteur 4, pour produire la représentation matricielle initiale (À) et les données de résidus, et agencé pour commander l'adaptateur (10) et le caîculateur-soiveur (12) en conséquence.
10. Dispositif selon l'une des revendications précédentes, dans lequel le système d'équations est représentatif d'un système physique complexe du monde réel, tel qu'un champ pétrolifère.
11. Procédé de calcul polyvalent du type comprenant :
a) recevoir une représentation matricielle initiale correspondant à un système d'équations à traiter et une représentation matricielle filtrante,
b) calculer une représentation matricielle de travail vérifiant avec la représentation matricielle initiale une condition de stabilité comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée, et comportant respectivement fa représentation matricielle initiale, et la représentation matricielle de travail,
c) recevoir des données de résidus, et résoudre le système d'équations défini par la représentation matricielle initiale, à partir des données de résidus, de la représentation matricielle de travail et de la représentation matricielle initiale,
caractérisé en ce que l'opération b) comprend calculer itérativement par blocs une matrice intermédiaire à partir de la représentation matricielle initiale et de ladite représentation numérique de représentation matricielle filtrante en répétant itérativement, pour chaque bloc courant (i j) de la matrice intermédiaire, :
bl) calculer une somme de blocs chacun défini par un produit faisant intervenir deux blocs déjà calculés de la matrice intermédiaire, et un bloc auxiliaire d'une matrice d'approche qui est contraint à vérifier avec un bloc diagonal déjà calculé de ia matrice intermédiaire une condition d'équivalence, comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée et un bloc déjà calculé de la matrice intermédiaire, et comportant respectivement l'inverse dudit bloc diagonal déjà calculé de la matrice intermédiaire, et ledit bloc auxiliaire de la matrice d'approche,
b2) calculer la différence entre le bloc (i,j) de la matrice de la représentation matricielle initiale et la somme issue de l'opération b l ), et ce que l'opération c) comprend travailler sur la matrice intermédiaire, par blocs, de façon à fournir une solution du système d'équations de la représentation matricielle initiale, sans inversion complète de celle-ci.
PCT/FR2011/051859 2010-08-03 2011-08-02 Dispositif informatique de calcul polyvalent WO2012017177A2 (fr)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US13/813,837 US20130226980A1 (en) 2010-08-03 2011-08-02 Multipurpose calculation computing device

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
FR1003248 2010-08-03
FR1003248A FR2963692A1 (fr) 2010-08-03 2010-08-03 Dispositif informatique de calcul polyvalent

Publications (2)

Publication Number Publication Date
WO2012017177A2 true WO2012017177A2 (fr) 2012-02-09
WO2012017177A3 WO2012017177A3 (fr) 2012-04-05

Family

ID=44227852

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/FR2011/051859 WO2012017177A2 (fr) 2010-08-03 2011-08-02 Dispositif informatique de calcul polyvalent

Country Status (3)

Country Link
US (1) US20130226980A1 (fr)
FR (1) FR2963692A1 (fr)
WO (1) WO2012017177A2 (fr)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9760537B2 (en) 2014-10-28 2017-09-12 International Business Machines Corporation Finding a CUR decomposition

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6081268A (en) * 1997-09-03 2000-06-27 Digital Equipment Corporation Method for ignoring redundant constraints in a graphic editor
GB0208329D0 (en) * 2002-04-11 2002-05-22 Univ York Data processing particularly in communication systems
WO2005121840A1 (fr) * 2004-06-07 2005-12-22 Exxonmobil Upstream Research Company Procede pour resoudre une equation matricielle de simulation de reservoir implicite
GB2455077A (en) * 2007-11-27 2009-06-03 Polyhedron Software Ltd Estimating the state of a physical system using generalized nested factorisation

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
ACHDOU Y ET AL: "Low frequency tangential filtering decomposition", NUMERICAL LINEAR ALGEBRA WITH APPLICATIONS, vol. 14, no. 2, mars 2007 (2007-03), pages 129-147, XP055002531, ISSN: 1070-5325, DOI: 10.1002/nla.512 *
Grigori L et al: "A class of multilevel parallel preconditioning strategies", INRIA Rapport de Recherche No. RR-7410, HAL: inria-00524110, version 1, 6 octobre 2010 (2010-10-06), XP055002490, Extrait de l'Internet: URL:http://hal.inria.fr/docs/00/52/41/10/PDF/Paper.pdf [extrait le 2011-07-11] *
Grigori L ET AL: "Generalized Filtering Decomposition", INRIA Rapport de Recherche No. RR-7569, HAL : inria-00576894, version 1, 15 mars 2011 (2011-03-15), XP055002495, Extrait de l'Internet: URL:http://hal.inria.fr/docs/00/57/68/94/PDF/RR-7569.pdf [extrait le 2011-07-11] *
GRIGORI L ET AL: "Two sides tangential filtering decomposition", JOURNAL OF COMPUTATIONAL AND APPLIED MATHEMATICS, vol. 235, no. 8, 15 février 2011 (2011-02-15), pages 2647-2661, XP055002485, ISSN: 0377-0427, DOI: 10.1016/j.cam.2010.11.016 *
Kumar P et al: "Combinative preconditioning based on Relaxed Nested Factorization and Tangential Filtering preconditioner", INRIA Rapport de Recherche No. RR-6955, HAL: inria-00392881, version 1, 9 juin 2009 (2009-06-09), XP055002493, Extrait de l'Internet: URL:http://hal.inria.fr/docs/00/39/28/81/PDF/RR-6955.pdf [extrait le 2011-07-11] *
Wang R et al: "A twisted block tangential filtering decomposition preconditioner", Mathematical Problems in Engineering, Volume 2009, Article ID 282307, 2009, XP055002573, DOI: 10.1155/2009/282307 Extrait de l'Internet: URL:http://www.emis.de/journals/HOA/MPE/Volume2009/282307.pdf [extrait le 2011-07-12] *

Also Published As

Publication number Publication date
FR2963692A1 (fr) 2012-02-10
WO2012017177A3 (fr) 2012-04-05
US20130226980A1 (en) 2013-08-29

Similar Documents

Publication Publication Date Title
EP0198729B1 (fr) Système de simulation d'un circuit électronique
EP0953168B1 (fr) Dispositif electronique de traitement de donnees-image, pour la simulation du comportement deformable d'un objet
EP3472736B1 (fr) Procédé d'estimation du facteur d'intensité des contraintes et procédé de calcul de durée de vie associé
FR3028031B1 (fr) Procede d'estimation d'un etat de navigation contraint en observabilite
WO2006021713A1 (fr) Procede d'ordonnancement de traitement de tâches et dispositif pour mettre en œuvre le procede
WO2016177959A1 (fr) Procede des simplification de modele de geometrie
FR2863052A1 (fr) Methode pour determiner les composantes d'un tenseur de permeabilite effectif d'une roche poreuse
FR3088464A1 (fr) Procédé de construction de réseau de neurones pour la simulation de systèmes physiques
FR3019344A1 (fr) Procede de construction d'un maillage optimise pour la simulation de reservoir dans une formation souterraine
EP1926033A1 (fr) Procédé, programme et système informatique de simulation de chenaux
Denoyelle Theoretical and numerical analysis of super-resolution without grid
EP1462605A1 (fr) Méthode de pseudoisation et d'éclatement pour décrire des fluides hydrocarbones
FR2966951A1 (fr) Procede de simulation pour determiner des coefficients aerodynamiques d'un aeronef
WO2012017177A2 (fr) Dispositif informatique de calcul polyvalent
FR2903794A1 (fr) Procede de modelisation de l'activite de commutation d'un circuit numerique
Almeida et al. Hybrid Evolvable Hardware for automatic generation of image filters
WO2012035272A1 (fr) Dispositif informatique de calcul polyvalent
EP3443369B1 (fr) Systeme et procede de test d'un circuit integre
EP0855038B1 (fr) Diagnostic de réseaux de composants, avec modélisation par bandes
FR3131413A1 (fr) Réseau de neurones avec génération à la volée des paramètres du réseau
WO2018229392A1 (fr) Dispositif de caracterisation et/ou de modelisation de temps d'execution pire-cas
EP2901411A1 (fr) Dispositif de decomposition d'images par transformee en ondelettes
Kozhanova High-Order Numerical Methods for Shock-Bubble Interaction Computations
Clavier Analytical and Geometric approches of non-perturbative Quantum Field Theories
Pinto Wavelet-based multiscale simulation of incompressible flows

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 11755394

Country of ref document: EP

Kind code of ref document: A2

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 13813837

Country of ref document: US

122 Ep: pct application non-entry in european phase

Ref document number: 11755394

Country of ref document: EP

Kind code of ref document: A2