J'utilise des liaisons de bibliothèque numérique pour Boost Ublas pour résoudre un système linéaire simple. L'option suivante fonctionne bien, sauf qu'il s'agit de la manipulation de matrices A (m x m) pour relativement Petit 'M'.
En pratique, j'ai une matrice beaucoup plus grande avec dimension m = 10 ^ 6 (jusqu'à 10 ^ 7).
Existe-t-il une approche C ++ existante pour la résolution de l'AX = B qui utilise efficacement la mémoire. P>
#include<boost/numeric/ublas/matrix.hpp>
#include<boost/numeric/ublas/io.hpp>
#include<boost/numeric/bindings/traits/ublas_matrix.hpp>
#include<boost/numeric/bindings/lapack/gesv.hpp>
#include <boost/numeric/bindings/traits/ublas_vector2.hpp>
// compileable with this command
//g++ -I/home/foolb/.boost/include/boost-1_38 -I/home/foolb/.boostnumbind/include/boost-numeric-bindings solve_Axb_byhand.cc -o solve_Axb_byhand -llapack
namespace ublas = boost::numeric::ublas;
namespace lapack= boost::numeric::bindings::lapack;
int main()
{
ublas::matrix<float,ublas::column_major> A(3,3);
ublas::vector<float> b(3);
for(unsigned i=0;i < A.size1();i++)
for(unsigned j =0;j < A.size2();j++)
{
std::cout << "enter element "<<i << j << std::endl;
std::cin >> A(i,j);
}
std::cout << A << std::endl;
b(0) = 21; b(1) = 1; b(2) = 17;
lapack::gesv(A,b);
std::cout << b << std::endl;
return 0;
}
6 Réponses :
Assumer que vos énormes matrices sont clairsemées, ce que j'espère qu'ils sont à cette taille, consultez le Pardiso < / a> Projet qui est un solveur linéaire rare, c'est ce dont vous aurez besoin si vous souhaitez gérer des matrices aussi grand que vous avez dit. Permet un stockage efficace de seulement des valeurs non nulles et est beaucoup plus rapide que de résoudre le même système de matrices denses.
Sans parler de la complexité du temps O (m ^ 3) de la solution naïve! Même le knuth intelligent à propos de O (m ^ 2.7ish) ... Si ces matrices ne sont pas rares, vous avez besoin d'un cluster et d'un analyseur numérique de première classe ...
+1 pour une idée de matrice clairsemée. J'ai trouvé des bibliothèques numériques et des comparaisons dans Pardiso Papier sur la comparaison des bibliothèques varous de matrice clairsemées ftp.numerical.rl.ac.uk/pub/Reports/ghsral200505.pdf Ceci peut être utilisé pour trouver d'autres bibliothèques matricielles rapides reconnues.
Pas sûr des implémentations C ++, mais il y a plusieurs choses que vous pouvez faire si la mémoire est un problème en fonction du type de matrice que vous utilisez:
@duffymo: merci. J'ai examiné la mise en œuvre de l'approche itérative en C ++, mais ils ont toujours besoin de la stocker dans matrice. freenet-homepage.de/guwi17/ublas/examples Si je me trompe, faites Vous connaissez une implémentation mem efficace en C ++ pour itératif?
Correct, bulgaRatbra. J'aurais dû me souvenir de ça. J'éditionnerais des algorithmes parallèles, car le problème de la partition de l'œuvre sur les processeurs N et de la tricoter pour obtenir le résultat est germe au problème de la déplaçant temporairement sur le disque.
Je suppose que votre matrice est dense. Si cela est rare, vous pouvez trouver de nombreux algorithmes spécialisés comme indiqué par Deusaduro et Duffyme .
Si vous n'avez pas de grappe (assez large) à votre disposition, vous souhaitez regarder des algorithmes hors de noyau. Scalapack a quelques solveurs extrêmes dans le cadre de son Paquet prototype , voir la documentation ICI et Google pour plus de détails. Recherche sur le Web pour les solvateurs / forfaits LU / (Matrix) «Out-of-Core) vous donnera des liens vers une richesse d'algorithmes et d'outils supplémentaires. Je ne suis pas un expert sur ceux.
Pour ce problème, la plupart des gens utiliseraient une cluster. Le colis que vous trouverez sur presque n'importe quel groupe est Scalapack, encore une fois. De plus, il existe généralement de nombreux autres emballages sur le cluster typique, de sorte que vous pouvez choisir et choisir ce qui convient à votre problème (exemples ICI et ici ).
Avant de commencer à coder, vous voulez probablement vérifier rapidement combien de temps il faudra pour résoudre votre problème. Un solveur typique prend environ O (3 * N ^ 3) (N est la dimension de la matrice). Si N = 100000, vous regardez donc 3000000 gflops. En supposant que votre solveur en mémoire soit 10 gflops / s par noyau, vous envisagez 3 1/2 jours sur un seul noyau. Au fur et à mesure que les algorithmes d'algorithmes, augmentent le nombre de cœurs devraient réduire le temps proche de linéairement. Au-dessus de cela vient les E / S.
CAVEAT: Le dessus O (3 * N ^ 3) suppose que vous utilisez des nombres complexes. Pour les nombres réels, divisez tout par 6, c'est-à-dire quelque part autour de O (0,5 * N ^ 3).
Jetez un coup d'œil au Liste des logiciels librement disponibles pour la solution de problèmes d'algèbre linéaire , compilés par Jack Dongarra et Hatem Ltaief.
Je pense que pour la taille du problème que vous regardez, vous avez probablement besoin d'un algorithme itératif. Si vous ne voulez pas stocker la matrice A dans un format clairsemé, vous pouvez utiliser une implémentation sans matrice. Les algorithmes itératifs n'ont généralement pas besoin d'accéder à des entrées individuelles de la matrice A, ils n'ont besoin que de calculer des produits de matrice-vectoriels AV (et parfois un ^ t v, le produit de la matrice transposée avec le vecteur). Donc, si la bibliothèque est bien conçue, elle devrait suffire si vous le transmettez une classe qui sait faire des produits de matrice-vectoriels.
Réponse courte: N'utilisez pas Boost's Réponse longue: ci-dessous est l'exemple de code (basé sur la fonction Une factorisation symbolique est effectuée avec avec LAPACK Liaisons, celles-ci ont été conçues pour des matrices denses,
Pas de matrices rares, utilisez UMFPACK à la place. UMFPACK est l'une des meilleures bibliothèques pour la résolution de l'AX = B lorsque A est grand et clairsemé.
umfpack_simple.c ) qui génère un simple A et B
et résout AX = B . generate_sparse_matrix_problem crée la matrice A et le
côté droit b . La matrice est d'abord construite en forme de triplet. Les
Les vecteurs TI, TJ et TX décrivent complètement A. Triplet Forme est facile à créer mais
Méthodes de matrice rares efficaces nécessitent un format de colonne clairsemé à compression. Conversion
est effectué avec umfpack_di_triplet_to_col . UMFPack_DI_SYMBOLIC . Un peu plat
La décomposition LU de A est effectuée avec umfpack_di_numérique .
Les résolveurs triangulaires inférieurs et supérieurs sont effectués avec umfpack_di_solve . n à 500 000, sur ma machine, l'ensemble du programme prend environ une seconde à courir. .
Valgrind repère que 369 239 649 octets (un peu plus de 352 Mo) ont été alloués.
Comme la réponse acceptée suggère qu'il ya Umfpack. Mais si vous utilisez Boost, vous pouvez toujours utiliser les matrices compactes dans Boost et utiliser Umfpack pour résoudre le système. Il y a une liaison qui facilite la tâche:
http://mathema.tican.de/software/boost-numérique-bindings
C'est environ deux ans daté, mais c'est juste une liaison (avec quelques autres).
voir la question connexe: Umfpack and Boost's Ublas Sparse Matrix
Soulignez l'évidence, une matrice que la matrice de taille est de 4x10 ^ 12 à 4x10 ^ 14 octets, ou 4 à 400 téraoctets pour une seule matrice seule. (À moins que, comme indiqué ci-dessous, ses stars)