Dlaczego naukowiec obliczeń musiałby wdrożyć własną wersję std :: complex?

Wiele bardziej znanych bibliotek C ++ w informatyce, takich jak Eigen , Trilinos i deal.II, używa standardowego obiektu biblioteki szablonów C ++ std::complex<>do reprezentowania złożonych liczb zmiennoprzecinkowych.

Gniazdo Poulson za odpowiedź na pytanie o domyślnych konstruktorów, wskazuje on, że ma własną implementację std::complexw Elemental „dla wielu powodów”. Jakie są te powody? Jakie są zalety i wady tego podejścia?

Uważam, że ta dyskusja wielokrotnie pojawiała się na liście PETSc. Moje główne powody to:

1. Standard C ++ stwierdza, że ​​std :: complex jest zdefiniowany tylko dla typów danych zmiennoprzecinkowych, podwójnych i długich podwójnych. Dlatego nie można go używać w przypadku innych typów danych, takich jak czteroprecyzja.

2. Norma nie gwarantuje stabilności złożonej arytmetyki.

3. Standard nie gwarantuje, że dane w std :: complex są przechowywane jako prawdziwy komponent, po którym następuje komponent urojony. Ma to kluczowe znaczenie dla interfejsów z bibliotekami zewnętrznymi, takimi jak BLAS i LAPACK. Dotyczy to wszystkich głównych wdrożeń, ale wolałbym móc to zapewnić.

4. Wolę móc bezpośrednio manipulować rzeczywistymi i wymyślonymi komponentami. std :: complex czyni to niepotrzebnie trudnym.

5. W końcu chciałbym mieć bardziej ogólną wersję, która wymaga, aby typ danych był pierścieniem zamiast pola. Obejmuje to liczby całkowite Gaussa.

Używam std::complex<>w swoich programach i muszę walczyć z flagami kompilatora i obejściem dla każdego nowego kompilatora lub aktualizacji kompilatora. Spróbuję opisać te walki w porządku chronologicznym:

1. Pomiary wydajności wykazały, że krok polegający jedynie na obliczeniu kwadratu wartości bezwzględnej pola liczb zespolonych zajął więcej czasu niż poprzednia FFT dla gcc-4.x. Wkopanie w wygenerowany kod asemblera pokazało, że std::norm( ) obliczyło wartość bezwzględną ( ) w sposób pozwalający uniknąć przepełnienia, a następnie podniosło wynik do kwadratu. Ten problem można rozwiązać za pomocą flagi kompilacji .| z $|z|^2$$|z|$-ffast-math
2. Kompilator intel icc na Linuksie (lub linkerze) skompilowany std::argdo non-opt w określonych konfiguracjach (zgodność linków z określoną wersją gcc). Problem pojawiał się zbyt często, więc std::argmusiał zostać zastąpiony przez atan2(imag(),real()). Ale zbyt łatwo było o tym zapomnieć podczas pisania nowego kodu.
3. Typ std::complexużywa innych konwencji wywołań (= ABI) niż wbudowany typ kompleksu C99 oraz wbudowany typ kompleksu Fortran dla nowszych wersji gcc.
4. Do -ffast-mathkompilacji współdziała flaga z obsługi wyjątków pływający punkt w nieoczekiwany sposób. Dzieje się tak, ponieważ kompilator wyciąga podziały z pętli, powodując division by zerowyjątki w czasie wykonywania. Te wyjątki nigdy nie miałyby miejsca w pętli, ponieważ odpowiedni podział nie miał miejsca z powodu otaczającej logiki. Ten był naprawdę zły, ponieważ była to biblioteka, która została skompilowana oddzielnie od programu, który używał obsługi wyjątków zmiennoprzecinkowych (używając różnych flag kompilacji) i napotkał te problemy (odpowiednie zespoły siedziały w przeciwnych częściach świata, więc ten problem naprawdę spowodował poważne problemy). Zostało to rozwiązane przez optymalizację wykonaną ręcznie przez kompilator z większą ostrożnością.
5. Biblioteka stała się częścią programu i nie używała już -ffast-mathflagi kompilacji. Po aktualizacji do nowszej wersji gcc wydajność spadła o ogromny czynnik. Nie zbadałem jeszcze szczegółowo tego problemu, ale obawiam się, że jest on związany z załącznikiem G C99 . Muszę przyznać, że jestem całkowicie zdezorientowany tą dziwną definicją mnożenia liczb zespolonych, a nawet wydaje się, że istnieją różne wersje tego z twierdzeniami, że inne wersje są błędne. Mam nadzieję, że -fcx-limited-rangeflaga kompilacji rozwiąże problem, ponieważ wydaje się, że istnieje inny problem związany z -ffast-mathtą nowszą wersją gcc.
6. -ffast-mathKompilacji flag sprawia zachowanie NaNcałkowicie nieprzewidywalny dla nowszych wersji gcc (nawet isnanma wpływ). Jedynym obejściem wydaje się być uniknięcie jakiegokolwiek wystąpienia NaNw programie, co przeczy celowi istnienia NaN.

Teraz możesz zapytać, czy planuję zrezygnować z wbudowanych typów złożonych i std::complexz tych powodów. Pozostanę przy typach wbudowanych tak długo, jak pozostanę przy C ++. W przypadku, gdy C ++ stanie się całkowicie bezużyteczny do obliczeń naukowych, wolę rozważyć zmianę języka na bardziej zajmujący się kwestiami związanymi z obliczeniami naukowymi.