Pourquoi un scientifique en informatique devrait-il implémenter sa propre version de std :: complex?

14

La plupart des bibliothèques C ++ les plus connues en science informatique telles que Eigen , Trilinos et deal.II utilisent l'objet de bibliothèque d'en-tête de modèle C ++ standard std::complex<>, pour représenter des nombres à virgule flottante complexes.

Dans la réponse de Jack Poulson à une question sur les constructeurs par défaut, il souligne qu'il a sa propre implémentation std::complexdans Elemental "pour un certain nombre de raisons". Quelles sont ces raisons? Quels sont les avantages et les inconvénients de cette approche?

Aron Ahmadia
la source

Réponses:

16

Je pense que cette discussion a été mentionnée à plusieurs reprises sur la liste PETSc. Mes principales raisons sont:

  1. La norme C ++ indique que std :: complex n'est défini que pour les types de données float, double et long double. Il ne peut donc pas être utilisé pour d'autres types de données, tels que la quadruple précision.

  2. La norme ne garantit pas la stabilité de l'arithmétique complexe.

  3. La norme ne garantit pas que les données d'un std :: complex soient stockées en tant que composant réel suivi du composant imaginaire. Ceci est crucial pour les interfaces avec des bibliothèques externes, telles que BLAS et LAPACK. C'est vrai pour toutes les implémentations majeures, mais je préférerais pouvoir le garantir.

  4. Je préfère pouvoir manipuler directement les composants réels et imaginaires. std :: complex rend cela inutilement difficile.

  5. Je voudrais éventuellement avoir une version plus générale qui ne nécessite que le type de données soit un anneau au lieu de nécessiter un champ. Cela comprendrait les entiers gaussiens.

Jack Poulson
la source
6
Le point 3 a été traité en C ++ 11. 26.4.4 indique que si zest une expression lvalue de type cv std::complex<T> alors reinterpret_cast<cv T(&)[2]>(z)et reinterpret_cast<cv T(&)[2]>(z)[0]doit désigner la partie réelle de z, et reinterpret_cast<cv T(&)[2]>(z)[1]doit désigner la partie imaginaire de z. Des tableaux de nombres complexes sont également traités.
James Custer
3
@JamesCuster: Je suis tout à fait favorable au passage au C ++ 11, mais les codes scientifiques qui veulent rester portables vers des architectures semi-exotiques devront probablement attendre au moins deux à trois ans de plus pour le faire. De plus, C ++ 11 ne résout malheureusement qu'une partie du problème.
Jack Poulson
Je comprends, je viens de le jeter là-bas au cas où quelqu'un se pencherait sur cette question à l'avenir.
James Custer
2
Eh bien, je pense que c'est une dérobade de dire qu'il faudrait attendre que les compilateurs prennent en charge C ++ 11. L'exigence explicite a été insérée dans la nouvelle norme car toutes les implémentations existantes la prennent déjà en charge. Je ne peux pas penser à un cas où il serait dangereux d'assumer déjà cette disposition particulière dans les compilateurs / bibliothèques existants car cela n'aurait tout simplement aucun sens d'implémenter std :: complex d'une autre manière.
Wolfgang Bangerth
1
@WolfgangBangerth: C'était plus un commentaire général sur le passage au C ++ 11. Quoi qu'il en soit, C ++ 11 ne résout pas la plupart des problèmes avec std :: complex.
Jack Poulson
7

J'utilise std::complex<>dans mes programmes, et je dois me battre avec des drapeaux de compilateur et une solution de contournement pour chaque nouveau compilateur ou mise à niveau du compilateur. Je vais essayer de raconter ces combats par ordre chronologique:

  1. std::norm|z|2|z|-ffast-math
  2. Le compilateur intel icc sur linux (ou éditeur de liens) compilé std::argen non-opt dans certaines configurations (compatibilité des liens avec une version spécifique de gcc). Le problème a refait surface trop souvent et a donc std::argdû être remplacé par atan2(imag(),real()). Mais il était trop facile d'oublier cela lors de l'écriture de nouveau code.
  3. Le type std::complexutilise des conventions d'appel différentes (= ABI) que le type complexe C99 intégré et le type complexe Fortran intégré pour les versions gcc plus récentes.
  4. L' -ffast-mathindicateur de compilation interagit avec la gestion des exceptions à virgule flottante de manière inattendue. Ce qui se passe, c'est que le compilateur extrait les divisions des boucles, provoquant ainsi des division by zeroexceptions lors de l'exécution. Ces exceptions ne se seraient jamais produites à l'intérieur de la boucle, car la division correspondante n'a pas eu lieu en raison de la logique environnante. Celui-là était vraiment mauvais, car c'était une bibliothèque qui était compilée séparément du programme qui utilisait la remise d'exception en virgule flottante (en utilisant différents drapeaux de compilation) et rencontrait ces problèmes (les équipes correspondantes étaient assis dans des parties opposées du monde, donc ce problème a vraiment causé de gros ennuis). Cela a été résolu en faisant l'optimisation utilisée par le compilateur à la main avec plus de soin.
  5. La bibliothèque est devenue une partie du programme et n'a plus utilisé l' -ffast-mathindicateur de compilation. Après une mise à niveau vers une nouvelle version de gcc, les performances ont chuté d'un facteur énorme. Je n'ai pas encore étudié cette question en détail, mais je crains qu'elle soit liée à l' annexe G de la C99 . Je dois admettre que je suis complètement confus par cette étrange définition de la multiplication pour les nombres complexes, et il semble même exister différentes versions de cela avec des affirmations selon lesquelles les autres versions sont erronées. J'espère que l' -fcx-limited-rangeindicateur de compilation résoudra le problème, car il semble y avoir un autre problème lié à -ffast-mathcette nouvelle version de gcc.
  6. L' -ffast-mathindicateur de compilation rend le comportement de NaNcomplètement imprévisible pour les nouvelles versions de gcc (même isnanest affecté). La seule solution de contournement semble être d'éviter toute occurrence de NaNdans le programme, ce qui va à l'encontre de l'objectif de l'existence de NaN.

Vous pouvez maintenant demander si je prévois d'abandonner les types complexes intégrés et std::complexpour ces raisons. Je vais rester avec les types intégrés, tant que je reste avec C ++. Dans le cas où le C ++ devrait réussir à devenir complètement inutilisable pour le calcul scientifique, je préférerais envisager de passer à un langage qui prend plus en charge les problèmes liés au calcul scientifique.

Thomas Klimpel
la source
On dirait que mes craintes liées à l'annexe G de C99 se sont réalisées, et -fcx-limited-range est maintenant en quelque sorte requis pour une vitesse de calcul décente lors de la multiplication de nombres complexes. C'est du moins ce que j'obtiens de l'histoire récente de la guerre: medium.com/@smcallis_71148/…
Thomas Klimpel