LOBPCG
Le Gradient conjugué préconditionné par bloc localement optimal (acronyme anglophone LOBPCG) est une méthode utilisée pour trouver les valeurs propres les plus grandes (ou les plus petites) et les vecteurs propres correspondants d'un problème de valeurs propres généralisé symétrique
pour une paire de matrices hermitiennes complexes ou réelles symétriques, où la matrice est également supposée définie positive.
Caractéristiques principales
- La méthode est dite matrix-free c'est-à-dire qu'elle ne nécessite pas de stocker explicitement les coefficients de la matrice, mais on peut accéder à la matrice en évaluant les produits matrice-vecteur[1].
- La méthode ne nécessite pas de factorisation de matrice, même pour un problème aux valeurs propres généralisé.
- Le coût de calcul par itération ainsi que l'occupation mémoire sont compétitifs par rapport à la méthode de Lanczos, calculant une seule paire vecteur propre, valeur propre d'une matrice symétrique.
- La convergence linéaire est théoriquement garantie et pratiquement observée.
- On observe une convergence accélérée due au préconditionnement direct, contrairement à la méthode de Lanczos, la matrice du préconditionneur peut être fixe ou varier d'une itération à l'autre mais rester symétrique et définie positive.
- La méthode peut tirer parti de l'utilisation de techniques de décomposition de domaine et d'un préconditionneur multigrille.
- La méthode peut être redémarrée à chaud (c'est-à-dire que l'on utilise des vecteurs propres approchés venant d'une autre méthode pour initier le calcul) et elle calcule une approximation des vecteurs propres, valeurs propres à chaque itération.
- La méthode est plus stable numériquement que la méthode de Lanczos, et peut fonctionner en arithmétique de faible précision.
- Elle est facile à mettre en œuvre, avec de nombreuses implantations existantes.
- La méthode repose sur la notion de bloc et permet ainsi d'utiliser des opérations matrice-matrice très efficaces, par exemple BLAS 3.
- La taille des blocs peut être ajustée pour équilibrer la vitesse de convergence par rapport aux coûts de calcul des orthogonalisations et au coût de la méthode Rayleigh-Ritz à chaque itération.
Algorithme
Version à un vecteur
On commence par décrire un algorithme qui recherché un seul couple vecteur propre / valeur propre, par exemple celui correspondant à la plus petite valeur propre.
Préliminaires : rappel sur la méthode de montée ou descente de gradient pour les problèmes de valeurs propres
La méthode effectue une maximisation (ou minimisation) itérative du quotient de Rayleigh généralisé
ce qui permet de trouver les paires vecteurs propres, valeurs propres les plus grandes (ou les plus petites) de
La direction utilisée par l'algorithme du gradient est précisément le gradient du quotient de Rayleigh généralisé, direction proportionnelle au vecteur
appelé vecteur propre résiduel. Si un préconditionneur est disponible, il est appliqué au résidu et donne le vecteur
appelé résidu préconditionné. Sans préconditionnement, on pose et donc . La méthode itérative s'écrit
ou, alternativement de manière moins concise,
Cette méthode est connue sous le nom de montée (ou descente) de gradient préconditionnée, où le scalaire est le pas d'itération. La taille optimale du pas peut être déterminée en maximisant le quotient de Rayleigh, c'est-à-dire,
(ou en cas de minimisation), auquel cas la méthode est dite localement optimale.
On note que le maximum (ou minimum) est pris pour appartenant à l'espace vectoriel de dimension deux engendré par les deux vecteurs et .
Formule de récurrence à trois termes
Pour accélérer considérablement la convergence de la montée (ou descente) de gradient préconditionnée localement optimale, un vecteur supplémentaire peut être ajouté à la formule de récurrence à deux termes pour la rendre à trois termes (devient une fonction de , et )
La maximisation/minimisation du quotient de Rayleigh se fait maintenant dans un espace vectoriel de dimension trois et peut être effectuée numériquement par la méthode Rayleigh-Ritz. L'ajout de vecteurs supplémentaires, voir, par exemple, l'extrapolation de Richardson, n'entraîne pas d'accélération significative[2] mais augmente les coûts de calcul, et n'est donc généralement pas recommandé.
Améliorations de la stabilité numérique
Au fur et à mesure que les itérations convergent, les vecteurs et deviennent presque linéairement dépendants, ce qui entraîne une perte de précision et rend la méthode Rayleigh-Ritz numériquement instable en présence d'erreurs d'arrondi. La perte de précision peut être évitée en substituant le vecteur avec un vecteur , qui peut être plus éloigné de , dans la base du sous-espace tridimensionnel , tout en gardant le sous-espace inchangé et en évitant l'orthonormalisation ou toute autre opération supplémentaire[2]. De plus, l'orthonormalisation de la base du sous-espace tridimensionnel peut être nécessaire pour les problèmes de valeurs propres mal conditionnés afin d'améliorer la stabilité et la précision atteignable.
Analogues du sous-espace de Krylov
Il s'agit d'une version à vecteur unique de la méthode LOBPCG - une possible généralisation des solveurs linéaires à gradient conjugué préconditionnés au cas des problèmes de valeurs propres symétriques.[2] Même dans le cas trivial et l'approximation résultante avec sera différent de celui obtenu par l'algorithme de Lanczos, bien que les deux approximations appartiendront au même sous-espace de Krylov.
Scénarios d'utilisation pratiques
L'extrême simplicité et la haute efficacité de la version à vecteur unique de LOBPCG la rendent attrayante pour les applications de calcul de valeurs propres en environnement contraint en ressources de calcul, par exemple pour la détection d'anomalies en temps réel basée sur le partitionnement spectral. On utilise aussi pour le partitionnement de graphe sur ASIC ou FPGA, ou encore pour la modélisation de phénomènes physiques sur les supercalculateurs exascale du classement TOP500.
version par bloc de l'algorithme
Résumé
On cherche maintenant à calculer plusieurs couples propres simultanément. On pourrait appliquer la méthode LOBPCG à un seul vecteur (décrite ci-dessus) pour avoir le couple correspondant à la plus grande valeur propre (en module), noté et enchaîner avec une déflation orthogonale pour trouver les couples suivant. Si on appelle la matrice (supposée symétrique, mais on peut généraliser on cas non-symétrique) dont on cherche les couples propres, rappelons que la déflation orthogonale consiste à introduire une matrice dont on peut montrer qu'elle a les mêmes couples propres que sauf . Il suffit d'appliquer à nouveau la méthode LOBPCG à un seul vecteur pour trouver un deuxième couple propre et de répéter cette opération autant de fois que nécessaire.
Cependant, avec approche, chaque nouveau couple propre calculé de manière approchée dépend du calcul de tous les autres couples propres précédents et ainsi les erreurs de calcul s'accumulent et les estimations des couples propres deviennent de moins en moins bonnes. Pour améliorer cette situation, la méthode LOBPCG par bloc se propose de calculer plusieurs couples propres ensemble d'un seul bloc[2], de manière rapide, précise et robuste des vecteurs propres. Elle améliore également les cas où l'on a plusieurs plusieurs couples propres très proches les uns des autres qui seraient très mal estimés par la déflation orthogonale (convergence très lente). La taille du bloc peut être ajustée pour équilibrer la stabilité numérique par rapport à la vitesse de convergence et aux coûts de calcul des orthonormalisations et de la méthode Rayleigh-Ritz à chaque itération.
Principe de la méthode
L'approche par blocs dans LOBPCG remplace les vecteurs uniques et avec des vecteurs-blocs, c'est-à-dire des matrices et , où, par exemple, chaque colonne de se rapproche de l'un des vecteurs propres. Toutes les colonnes sont itérées simultanément, et la prochaine matrice de vecteurs propres approximatifs est déterminé par la méthode Rayleigh-Ritz sur le sous-espace couvert par toutes les colonnes de matrices et . Chaque colonne de est calculée simplement comme le résidu préconditionné pour chaque colonne de La matrice est déterminée de telle sorte que les sous-espaces engendrés par les colonnes de et de sont identiques.
Compromis entre stabilité numérique et efficacité
La méthode Rayleigh-Ritz est appliquée sur le sous-espace engendré par toutes les colonnes de matrices et , où en principe une base peut être choisie arbitrairement. Cependant, la méthode Rayleigh-Ritz peut en pratique devenir numériquement instable si certains des vecteurs de base ne sont qu'approximativement linéairement indépendant.
Le choix de la base de sous-espace devient crucial autant pour d'assurer la stabilité numérique de la méthode Rayleigh-Ritz que pour réduire à des coûts de calcul. L'approche sans doute la plus stable consiste à rendre les vecteurs de base orthogonaux, par exemple, par le processus d'orthonormalisation de Gram-Schmidt, opération coûteuse en temps de calcul. Par exemple, les implémentations LOBPCG[3],[4], utilise une factorisation de Cholesky instable mais efficace de la matrice normale, qui est effectuée uniquement sur les matrices individuelles et , plutôt que sur l'ensemble du sous-espace. Sur une architecture de calcul moderne, la quantité de mémoire disponible est suffisante pour se permettre d'utiliser des tailles de bloc de l'ordre de , néanmoins le temps de calcul consacré aux orthogonalisations et à la méthode Rayleigh-Ritz commence à dominer.
Verrouillage des vecteurs propres précédemment convergés
Quand on utilise une méthode par bloc itérative pour approcher les vecteurs propres, il arrive très souvent qu'à l'intérieur d'un bloc, certains vecteurs propres convergent plus rapidement que d'autres. Le verrouillage des vecteurs propres déjà convergés consiste à les retirer de la boucle itérative, afin d'éliminer des calculs inutiles et d'améliorer la stabilité numérique. Une simple suppression d'un vecteur propre peut probablement entraîner la formation de son double dans des vecteurs toujours itératifs. Le fait que les vecteurs propres d'une matrice symétrique soient orthogonaux par paire suggère de garder tous les vecteurs itératifs orthogonaux aux vecteurs verrouillés.
Le verrouillage peut être mis en œuvre différemment en maintenant la précision et la stabilité numériques tout en minimisant les coûts de calcul. Par exemple, les implémentations LOBPCG[3],[4],[2],[5], séparent le verrouillage dur, c'est-à-dire une déflation par restriction, où les vecteurs propres verrouillés servent d'entrée de code et ne changent pas, du verrouillage doux, où les vecteurs verrouillés ne participent pas à l'étape itérative généralement la plus coûteuse en calcul des résidus, mais cependant, participent pleinement à la méthode Rayleigh-Ritz et peuvent donc être modifiés par la méthode Rayleigh-Ritz.
Théorie et pratique de la convergence
La méthode LOBPCG, par construction, assure[2] une minimisation du quotient de Rayleigh au moins aussi rapide que la descente de gradient appliqué à un du bloc, dont on connaît la convergence théorique. Chaque vecteur propre est un point stationnaire du quotient de Rayleigh, où le gradient s'annule. Ainsi, la descente du gradient peut ralentir à proximité de n'importe quel vecteur propre, cependant, soit la méthode converge vers le vecteur propre avec un taux de convergence linéaire ou, si ce vecteur propre est un point col, le quotient de Rayleigh itératif est susceptible de converger vers la deuxième plus basse valeur propre. La convergence dans le pire cas a été déterminée[2] et dépend de l'écart relatif entre la valeur propre et le reste du spectre de la matrice et de la qualité du préconditionneur, s'il est présent.
Pour une matrice quelconque, il n'y a évidemment aucun moyen de prédire les vecteurs propres et donc de générer les premières approximations qui fonctionnent toujours bien. La solution itérative de LOBPCG peut être sensible aux valeurs initiales des vecteurs propres, et par exemple, mettre plus de temps à converger et ralentir lorsque les paires propres suivantes doivent être calculées. De plus, en théorie, on ne peut pas nécessairement garantir la convergence vers le plus petit couple propre, bien que la probabilité de l'échec soit nulle. Une fonction gaussienne aléatoire de bonne qualité avec une moyenne nulle est généralement prise pour générer les approximations initiales de la méthode LOBPCG. Pour fixer les valeurs approchées initiales, on peut choisir une graine fixée pour le générateur de nombres aléatoires.
Contrairement à la méthode de Lanczos, LOBPCG présente rarement une convergence superlinéaire asymptotique dans la pratique.
Application à l'analyse en composantes principales (ACP) et à la décomposition en valeur singulière (SVD)
la méthode LOBPCG peut être facilement adaptée au calcul de plusieurs plus grandes valeurs singulières et les vecteurs singuliers correspondants et ainsi être utilisée pour l'analyse en composante principale (ACP). Elle est implantée dans le module Python SciPy depuis la version 1.4.0[6].
Implémentations logicielles
Andrew Knyazev, l'auteur original de la méthode LOBPCG, a publié une implantation logicielle de référence appelée Block Locally Optimal Eigenvalue Xolvers (BLOPEX)[7],[8] avec des interfaces vers les bibliothèques d'algèbre linéaire parallèle PETSc, hypre et Parallel Hierarchical Adaptive MultiLevel method (PHAML)[9]. D'autres implantations sont disponibles dans les logiciels GNU Octave[10], MATLAB[3], Java[11], Anasazi (Trilinos)[12], SLEPc[13],[14], SciPy[4], Julia[15], MAGMA[16], PyTorch[17], Rust[18], OpenMP et OpenACC[19], CuPy (Une bibliothèque de tableaux compatible NumPy pour les processeurs graphiques)[20], et NVIDIA AMGX[21]. Il existe également une implantation de LOBPCG compatible avec TensorFlow[22],
Applications
Les progiciels scikit-learn et Megaman[23] utilisent LOBPCG pour adapter le regroupement spectral[24] et l'apprentissage automatique en grande dimension[25] pour les problèmes de réduction de la dimension. NVIDIA a implanté[26] la méthode LOBPCG dans sa bibliothèque nvGRAPH introduite dans CUDA 8. le logiciel Sphynx[27], implante un algorithme de partitionnement de graphe parallèle hybride (CPU et GPU) à mémoire distribuée et partagée, il utilise le regroupement spectral pour le partitionnement de graphe, en calculant les vecteurs propres sur la matrice laplacienne du graphe en utilisant la méthode LOBPCG du package Anasazi de la bibliothèques Trilinos.
LOBPCG est implémenté dans ABINIT[28] (y compris la version CUDA) et Octopus[29]. Il a été utilisé pour des matrices de plusieurs milliards de dollars par les finalistes du prix Gordon Bell, sur le supercalculateur Earth Simulator au Japon[30],[31]. Le modèle Hubbard pour les systèmes électroniques fortement corrélés pour comprendre le mécanisme derrière la supraconductivité utilise LOBPCG pour calculer l'état fondamental de l'hamiltonien sur l'ordinateur K.[32] Il existe des versions MATLAB[33] et Julia[34],[35],[36] de LOBPCG pour les équations de Kohn-Sham et la théorie de la fonctionnelle de la densité (DFT) utilisant la base des ondes simples. Les implémentations récentes incluent TTPY[37], Platypus‐QM[38], MFDn[39], ACE-Molecule[40], LACONIC.
L'implantation de LOBPCG de la bibliothèque BLOPEX est utilisée pour la configuration du préconditionneur dans la bibliothèque de solveurs BDDCML (Multilevel Balancing Domain Decomposition by Constraints), qui est incorporé dans OpenFTL (Open Finite Element Template Library) et utilisé par le simulateur Flow123d pour étudier les écoulements d'eau souterraine, de solutés et de transport de chaleur dans les milieux poreux fracturés. LOBPCG a été implémenté[41] dans LS-DYNA.
LOBPCG est l'un des principales méthodes de recherche de valeurs propres dans PYFEMax et dans le logiciel d'éléments finis multiphysique haute performance Netgen/NGSolve. L'implantation de LOBPCG fournie par la bibliothèque hypre est intégré à la bibliothèque C++ MFEM (simulation par la méthode des éléments finis), qui est utilisée dans de nombreux projets, notamment BLAST, XBraid, VisIt, xSDK, et par l'institut FASTMath du SciDAC pour la conception de discrétisations exascales efficaces (CEED).
Un filtre passe-bas approché basé sur les itérations de la méthode LOBPCG peut être utilisé pour le débruitage ; voir[42], par exemple, en traitement d'images.
La segmentation d'images par regroupement spectral effectue un plongement de faible dimension à l'aide d'une matrice d'affinité entre les pixels, suivi d'un regroupement des composants des vecteurs propres dans l'espace de faible dimension, par exemple en utilisant le graphe Laplacien pour le filtrage bilatéral. La segmentation d'images via un partitionnement de graphe spectral par LOBPCG avec préconditionnement multigrille a été proposée pour la première fois dans[43] et actuellement testée dans[44],[45]. Cette dernière approche a ensuite été implémentée dans Python Scikit-learn[46] qui utilise l'implantation de LOBPCG fournie par SciPy avec un préconditionnement algébrique multigrille pour résoudre le problème des valeurs propres pour le graphe Laplacien.
Références
- ↑ (en) Auteur inconnu, « Recent implementations, applications, and extensions of the Locally Optimal Block Preconditioned Conjugate Gradient method (LOBPCG) », ..
- ↑ a b c d e f et g Knyazev, « Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method », SIAM Journal on Scientific Computing, vol. 23, no 2, , p. 517–541 (DOI 10.1137/S1064827500366124, lire en ligne).
- ↑ a b et c MATLAB File Exchange function LOBPCG.
- ↑ a b et c SciPy sparse linear algebra function lobpcg.
- ↑ A. Knyazev « Hard and soft locking in iterative methods for symmetric eigenvalue problems » () (DOI 10.13140/RG.2.2.11794.48327, lire en ligne)
—Eighth Copper Mountain Conference on Iterative Methods March 28 - April 2, 2004. - ↑ LOBPCG for SVDS in SciPy.
- ↑ GitHub BLOPEX.
- ↑ Knyazev, Argentati, Lashuk et Ovtchinnikov, « Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in Hypre and PETSc », SIAM Journal on Scientific Computing, vol. 29, no 5, , p. 2224 (DOI 10.1137/060661624, Bibcode 2007arXiv0705.2626K, arXiv 0705.2626).
- ↑ PHAML BLOPEX interface to LOBPCG.
- ↑ Octave linear-algebra function lobpcg.
- ↑ Java LOBPCG at Google Code.
- ↑ Anasazi Trilinos LOBPCG at GitHub.
- ↑ Native SLEPc LOBPCG.
- ↑ SLEPc BLOPEX interface to LOBPCG.
- ↑ Julia LOBPCG at GitHub.
- ↑ Anzt, Tomov et Dongarra, « Accelerating the LOBPCG method on GPUs using a blocked sparse matrix vector product », Proceedings of the Symposium on High Performance Computing (HPC '15). Society for Computer Simulation International, San Diego, CA, USA, hPC '15, , p. 75–82 (ISBN 9781510801011, lire en ligne).
- ↑ PyTorch LOBPCG at GitHub.
- ↑ Rust LOBPCG at GitHub.
- ↑ Fazlay Rabbi, Christopher S. Daley, Hasan M. Aktulga et Nicholas J. Wright « Evaluation of Directive-based GPU Programming Models on a Block Eigensolver with Consideration of Large Sparse Matrices » () (lire en ligne).
- ↑ CuPy: A NumPy-compatible array library accelerated by CUDA LOBPCG at GitHub.
- ↑ NVIDIA AMGX LOBPCG at GitHub.
- ↑ Rakhuba, Novikov et Osedelets, « Low-rank Riemannian eigensolver for high-dimensional Hamiltonians », Journal of Computational Physics, vol. 396, , p. 718–737 (DOI 10.1016/j.jcp.2019.07.003, Bibcode 2019JCoPh.396..718R, arXiv 1811.11049, lire en ligne).
- ↑ McQueen et al., « Megaman: Scalable Manifold Learning in Python », Journal of Machine Learning Research, vol. 17, no 148, , p. 1–5 (Bibcode 2016JMLR...17..148M, lire en ligne).
- ↑ « Sklearn.cluster.SpectralClustering — scikit-learn 0.22.1 documentation ».
- ↑ « Sklearn.manifold.spectral_embedding — scikit-learn 0.22.1 documentation ».
- ↑ Naumov, « Fast Spectral Graph Partitioning on GPUs », NVIDIA Developer Blog, (lire en ligne).
- ↑ « SGraph partitioning with Sphynx ».
- ↑ ABINIT Docs: WaveFunction OPTimisation ALGorithm.
- ↑ Octopus Developers Manual:LOBPCG.
- ↑ S. Yamada, T. Imamura et M. Machida « 16.447 TFlops and 159-Billion-dimensional Exact-diagonalization for Trapped Fermion-Hubbard Model on the Earth Simulator » () (DOI 10.1109/SC.2005.1).
- ↑ S. Yamada, T. Imamura, T. Kano et M. Machida « Gordon Bell finalists I—High-performance computing for exact numerical approaches to quantum many-body problems on the earth simulator » () (DOI 10.1145/1188455.1188504)
—Proc. ACM/IEEE conference on Supercomputing (SC '06). - ↑ S. Yamada, T. Imamura et M. Machida « High Performance LOBPCG Method for Solving Multiple Eigenvalues of Hubbard Model: Efficiency of Communication Avoiding Neumann Expansion Preconditioner » () (DOI 10.1007/978-3-319-69953-0_14)
—Asian Conference on Supercomputing Frontiers. - ↑ Yang, Meza, Lee et Wang, « KSSOLV - a MATLAB toolbox for solving the Kohn-Sham equations », ACM Trans. Math. Softw., vol. 36, no 2, , p. 1–35 (DOI 10.1145/1499096.1499099).
- ↑ Fathurrahman, Agusta, Saputro et Dipojono, « PWDFT.jl: A Julia package for electronic structure calculation using density functional theory and plane wave basis », Computer Physics Communications, vol. 256, , p. 107372 (DOI 10.1016/j.cpc.2020.107372, Bibcode 2020CoPhC.25607372F).
- ↑ Plane wave density functional theory (PWDFT) in Julia.
- ↑ Density-functional toolkit (DFTK). Plane wave Density Functional Theory in Julia.
- ↑ Rakhuba et Oseledets, « Calculating vibrational spectra of molecules using tensor train decomposition », J. Chem. Phys., vol. 145, no 12, , p. 124101 (PMID 27782616, DOI 10.1063/1.4962420, Bibcode 2016JChPh.145l4101R, arXiv 1605.08422).
- ↑ Takano, Nakata, Yonezawa et Nakamura, « Development of massive multilevel molecular dynamics simulation program, platypus (PLATform for dYnamic protein unified simulation), for the elucidation of protein functions », J. Comput. Chem., vol. 37, no 12, , p. 1125–1132 (PMID 26940542, PMCID 4825406, DOI 10.1002/jcc.24318).
- ↑ Shao et al., « Accelerating Nuclear Configuration Interaction Calculations through a Preconditioned Block Iterative Eigensolver », Computer Physics Communications, vol. 222, no 1, , p. 1–13 (DOI 10.1016/j.cpc.2017.09.004, Bibcode 2018CoPhC.222....1S, arXiv 1609.01689).
- ↑ Kang et al., « ACE-Molecule: An open-source real-space quantum chemistry package », The Journal of Chemical Physics, vol. 152, no 12, , p. 124110 (PMID 32241122, DOI 10.1063/5.0002959, Bibcode 2020JChPh.152l4110K).
- ↑ « A Survey of Eigen Solution Methods in LS-DYNA® » () (lire en ligne)
—15th International LS-DYNA Conference, Detroit. - ↑ A. Knyazev et A. Malyshev « Accelerated graph-based spectral polynomial filters » () (DOI 10.1109/MLSP.2015.7324315, arXiv 1509.02468)
—2015 IEEE 25th International Workshop on Machine Learning for Signal Processing (MLSP), Boston, MA. - ↑ Andrew V. Knyazev « Modern preconditioned eigensolvers for spectral image segmentation and graph bisection » () (lire en ligne)
—Clustering Large Data Sets; Third IEEE International Conference on Data Mining (ICDM 2003) Melbourne, Florida: IEEE Computer Society. - ↑ Andrew V. Knyazev « Multiscale Spectral Image Segmentation Multiscale preconditioning for computing eigenvalues of graph Laplacians in image segmentation » () (DOI 10.13140/RG.2.2.35280.02565, lire en ligne)
—Fast Manifold Learning Workshop, WM Williamburg, VA. - ↑ Andrew V. Knyazev « Multiscale Spectral Graph Partitioning and Image Segmentation » () (lire en ligne)
—Workshop on Algorithms for Modern Massive Datasets Stanford University and Yahoo! Research. - ↑ « Spectral Clustering — scikit-learn documentation ».
Liens externes
- LOBPCG dans MATLAB
- LOBPCG à Octave
- LOBPCG dans SciPy
- LOBPCG en Java chez Google Code
- LOBPCG dans Block Locally Optimal Eigenvalue Xolvers (BLOPEX) sur GitHub et archivé sur Google Code
Portail des mathématiques