No 3D Matrices: A Unified Tensor-Product View of Matrix-Free Cartesian PDE Solvers

Résumé Technique : Pas de matrices 3D : Une vue unifiée par produit tensoriel des solveurs d'EDP cartésiens sans matrice

Énoncé du problème
Les codes de mécanique des fluides numérique (CFD) de production et les solveurs d'éléments spectraux reposent depuis longtemps sur un secret structurel : les opérateurs tridimensionnels régissant les EDP cartésiennes (tels que le Laplacien, les matrices de masse, et les schémas de dérivées compactes) ne sont jamais assemblés, factorisés ou stockés sous forme de matrices monolithiques $N \times N$, où $N = N_x N_y N_z$. Bien que les manuels présentent souvent ces opérateurs comme des matrices 3D denses ou creuses, le faire est informatiquement prohibitif. Pour une grille de 200 points par côté, un Laplacien dense nécessiterait environ un demi-pétaoctet de stockage, et même une factorisation directe par stencil à 7 points d'une matrice creuse épuiserait la mémoire d'une station de travail. Le problème traité est le décalage entre la représentation théorique de ces opérateurs et l'implémentation efficace « sans matrice » utilisée en pratique. L'article cherche à unifier la justification algébrique de cet « effondrement » des opérateurs 3D en opérations 1D et à codifier les astuces d'implémentation spécifiques qui rendent cette approche efficace sur le matériel.

Méthodologie
L'article dérive un cadre unifié basé sur l'algèbre des produits de Kronecker. L'idée centrale est que tout opérateur 3D cartésien $L$ peut être exprimé comme une somme de produits tensoriels d'opérateurs monodimensionnels :
$$L = \sum_{r} C_z^{(r)} \otimes C_y^{(r)} \otimes C_x^{(r)}$$
où $C_\xi$ sont de petites matrices de bande ou diagonales agissant le long des axes de coordonnées. La méthodologie repose sur deux identités fondamentales :

1. La règle du produit mixte : $(A \otimes B)(C \otimes D) = (AC) \otimes (BD)$, ce qui permet la composition d'opérateurs directionnels sans former leur produit.
2. La règle de l'inverse : $(A \otimes B)^{-1} = A^{-1} \otimes B^{-1}$, ce qui explique pourquoi l'inversion d'une matrice de masse ou d'un facteur compact produit tensoriel se réduit à des résolutions 1D indépendantes.

L'article applique systématiquement cette algèbre à :

• Les différences finies : Différences centrales et stencils d'ordre supérieur.
• Les schémas compacts (Padé) : Relations implicites $A_x u' = R_x u$, évitant la formation du produit dense $A_x^{-1}R_x$.
• Les méthodes de Galerkin et spectrales : Incluant les B-splines, l'analyse isogéométrique (IGA) et le colocalisation, montrant que les matrices de masse et de rigidité héritent de la structure de Kronecker des bases en produit tensoriel.
• Le pas de temps : Les méthodes de direction alternée implicite (ADI) (Peaceman–Rachford, Douglas–Gunn) sont montrées comme décomposant une étape implicite 3D couplée en trois résolutions 1D séquentielles de bandes.
• Les solveurs directs : La diagonalisation rapide (Lynch–Rice–Thomas) pour les équations de Poisson et de Helmholtz, où le problème 3D est découplé en divisions scalaires dans le sous-espace propre des opérateurs 1D.

Contributions principales
L'article apporte trois contributions majeures, allant au-delà de la réduction algébrique pour atteindre les « astuces de production » nécessaires au calcul de haute performance :

1. Dérivation algébrique unifiée : Il fournit une dérivation unique et limpide montrant que la réduction aux opérations 1D s'applique universellement aux différences finies, aux schémas compacts, aux méthodes de Galerkin, aux B-splines et à la colocalisation. Il clarifie que l'« opérateur 3D » est une abstraction mathématique qui n'a jamais besoin d'être instanciée en mémoire.
2. Trois astuces d'implémentation de production :
   • Remodelage Multi-RHS (Multi-Right-Hand-Side) : Au lieu de boucler sur $N_y N_z$ lignes 1D indépendantes, le champ est remodelé en une matrice de taille $N_x \times (N_y N_z)$. Cela présente l'opération comme un noyau groupé unique, permettant l'utilisation d'appels BLAS-3 GEMM optimisés pour les facteurs denses/locaux à l'élément, ou de solveurs de bandes groupés pour les stencils.
   • Factorisation de somme (Sum Factorization) : Pour les méthodes d'éléments spectraux d'ordre élevé, l'article détaille comment réduire le coût de la quadrature de $O(n^{2d})$ à $O(d n^{d+1})$ (où $n$ est le nombre de points par direction) en effectuant des contractions 1D séquentielles plutôt qu'un produit tensoriel 3D complet.
   • Décomposition en crayons (Pencil Decomposition) : Pour la mémoire distribuée (MPI), la grille est décomposée en « crayons » (plaques) alignés avec une direction. Des transpositions collectives sont utilisées pour faire pivoter l'orientation du crayon, garantissant que chaque direction de balayage est contiguë en mémoire pour les rangs locaux.
3. Analyse de complexité : L'article établit rigoureusement que pour une largeur de stencil ou un degré polynomial fixe, le coût de calcul est $O(N)$ et le stockage de l'opérateur est $O(N_x + N_y + N_z)$, contrastant nettement avec les coûts de stockage $O(N)$ et de factorisation $O(N^2)$ ou $O(N^{4/3})$ des matrices 3D assemblées.

Résultats
L'article valide ses affirmations par une comparaison numérique d'une résolution de Poisson 3D sur un cube unité :

• Stockage : Pour $N \approx 1,7 \times 10^7$ inconnues ($256^3$), le solveur direct creux assemblé nécessite des gigaoctets de remplissage (et échoue à s'exécuter au-delà de $48^3$ sur une station de travail standard), tandis que la méthode de diagonalisation rapide sans matrice ne nécessite que des kilo-octets de stockage (les valeurs propres 1D).
• Temps : Le solveur sans matrice termine le problème $256^3$ en environ 2 secondes, tandis que le solveur direct assemblé prend près d'une minute pour $48^3$ et devient infaisable pour des grilles plus grandes.
• Précision : Les deux méthodes produisent des solutions discrètes identiques (accord à $\approx 10^{-14}$), confirmant que l'approche sans matrice n'introduit aucune perte de précision.

Signification
L'article affirme que sa signification réside dans le fait de « mettre en lumière » le secret structurel qui est resté une « habitude non écrite » des spécialistes de code depuis un demi-siècle. En présentant le produit de Kronecker comme le principe organisateur, il unifie des méthodes de discrétisation disparates (différences finies, splines, éléments spectraux) sous un seul paradigme algorithmique. Il soutient que la « matrice 3D » est un artefact de notation, et non une composante nécessaire du calcul. L'article sert de manuel pour implémenter ces réductions, fournissant une base de référence $O(N)$ précise contre laquelle les futurs accélérateurs sans matrice, à faible rang ou appris, pourront être jugés. Il souligne que l'efficacité des solveurs cartésiens modernes ne provient pas d'algorithmes nouveaux et complexes, mais de l'exploitation correcte de la structure séparable de la grille et de l'opérateur via des stratégies spécifiques de disposition des données et de noyaux.