Implémentation efficace d'un algorithme à matrice tridiagonale

12

Je résous un problème physique en utilisant un schéma numérique implicite. Cela m'amène à résoudre une équation linéaire avec une matrice tridiagonale. J'ai codé cet algorithme à partir de Wikipedia. Je me demande s'il existe une bibliothèque efficace qui permet de résoudre ce type d'équation de manière optimisée. Une note importante est que la matrice elle-même ne change que lorsque les paramètres du système changent, j'ai donc eu l'occasion de précalculer certaines étapes de l'algorithme pour un bon bonus de performance. J'utilise C ++.

linear-algebra libraries c++ gmk
la source
Quelle est la taille d'un système, doit-il être parallèle?
aterrel
1
La taille dépend de la précision requise (de centaines à des dizaines de milliers de valeurs). Maintenant, je code sur un ordinateur monocœur, mais il est possible d'avoir accès à un superordinateur universitaire avec de nombreux processeurs disponibles, donc le support du parallélisme serait bien.
gmk

Réponses:

15

Vous devriez probablement commencer par l'implémentation de LAPACK,? Gtsv, par exemple, dgtsv . Si vous voulez une version à mémoire distribuée, vous voudrez peut-être commencer par p? Gtsv de ScaLAPACK.

EDIT: Étant donné que votre matrice ne change pas très souvent, vous pouvez éviter de factoriser de manière redondante la matrice tridiagonale en décomposant la routine LAPACK? Gtsv en l'étape de factorisation,? Gttrf, et l'étape de résolution,? Gttrs. Des routines portant le même nom existent dans ScaLAPACK et servent le même objectif.

Jack Poulson
la source
Merci, ça ressemble à ce dont j'ai besoin. Je vais essayer maintenant d'exécuter ces routines à partir de mon code.
gmk
1
Puisque vous l'appelez depuis C ++, assurez-vous de déclarer le prototype à l'intérieur d'un bloc externe "C" {}. Selon votre système, vous devrez peut-être ajouter un trait de soulignement au nom de la routine.
Jack Poulson
2

Pour les systèmes parallèles distribués : je n'ai pas essayé ScaLAPACK, qui a un solveur tridiagonal parallèle, pour lequel il existe des exemples disponibles en ligne . J'ai essayé avec un certain succès une méthode proposée par David Moulton dans une publication LANL . Coder cela peut être plus que ce que vous voulez faire, mais en utilisant LAPACK, c'est bien en avant.

Yann
la source
1

J'ai trouvé un algorithme récursif intéressant ici à la page 975. Il semble prometteur, je me demande ce que les gens plus expérimentés en disent.

tiam
la source
Recettes numériques contient des erreurs. En termes de source de codes à utiliser, ce n'est pas le meilleur, bien que certains le considèrent comme un classique. Je serais surpris si ScaLAPACK n'implémentait pas un algorithme au moins aussi efficace que la réduction cyclique récursive.
Geoff Oxberry