¿Por qué necesitaría un científico computacional implementar su propia versión de std :: complex?

14

Muchas de las bibliotecas de C ++ más conocidas en ciencia computacional como Eigen , Trilinos y deal.II usan el objeto de biblioteca de encabezado de plantilla C ++ estándar std::complex<>, para representar números complejos de punto flotante.

En la respuesta de Jack Poulson a una pregunta sobre constructores predeterminados, señala que tiene su propia implementación std::complexen Elemental "por varias razones". ¿Cuáles son esas razones? ¿Cuáles son las ventajas y desventajas de este enfoque?

Aron Ahmadia
fuente

Respuestas:

16

Creo que esta discusión ha aparecido varias veces en la lista PETSc. Mis razones principales son:

  1. El estándar C ++ establece que std :: complex solo se define para los tipos de datos float, double y long double. Por lo tanto, no se puede utilizar para otros tipos de datos, como la precisión cuádruple.

  2. La norma no garantiza la estabilidad de la aritmética compleja.

  3. El estándar no garantiza que los datos en un std :: complex se almacenen como el componente real seguido por el componente imaginario. Esto es crucial para las interfaces con bibliotecas externas, como BLAS y LAPACK. Es cierto para todas las implementaciones principales, pero preferiría poder asegurarlo.

  4. Prefiero poder manipular directamente los componentes reales e imaginarios. std :: complex hace que esto sea innecesariamente difícil.

  5. Finalmente, me gustaría tener una versión más general que solo requiera que el tipo de datos sea un anillo en lugar de requerir un campo. Esto incluiría los enteros gaussianos.

Jack Poulson
fuente
66
El punto 3 se ha abordado en C ++ 11. 26.4.4 establece que si zes una expresión de tipo lvalue cv std::complex<T> a continuación, reinterpret_cast<cv T(&)[2]>(z)y reinterpret_cast<cv T(&)[2]>(z)[0]designará la parte real de z, y reinterpret_cast<cv T(&)[2]>(z)[1]designará la parte imaginaria de z. También se abordan las matrices de números complejos.
James Custer
3
@JamesCuster: Estoy dispuesto a cambiar eventualmente a C ++ 11, pero los códigos científicos que desean permanecer portátiles a arquitecturas semi-exóticas probablemente tendrán que esperar al menos dos o tres años más para hacerlo. Además, C ++ 11 desafortunadamente solo aborda parte del problema.
Jack Poulson
Entiendo, lo estaba lanzando por si alguien mira esta pregunta en el futuro.
James Custer
2
Bueno, creo que es un error decir que tendrías que esperar hasta que los compiladores admitan C ++ 11. El requisito explícito se incluyó en el nuevo estándar porque todas las implementaciones existentes ya lo admiten. No puedo pensar en un caso en el que no sea seguro asumir este diseño en particular en los compiladores / bibliotecas existentes, ya que simplemente no habría tenido ningún sentido implementar std :: complex de ninguna otra manera.
Wolfgang Bangerth
1
@WolfgangBangerth: Fue más un comentario general sobre cambiar a C ++ 11. De cualquier manera, C ++ 11 no soluciona la mayoría de los problemas con std :: complex.
Jack Poulson
7

Lo uso std::complex<>en mis programas y tengo que luchar con los indicadores del compilador y la solución para cada nuevo compilador o actualización del compilador. Intentaré contar estas peleas en orden cronológico:

  1. Las mediciones de rendimiento mostraron que un paso que involucraba solo calcular el cuadrado del valor absoluto de un campo de números complejos tomó más tiempo que un FFT anterior para gcc-4.x. Excavar en el código ensamblador generado mostró que std::norm(El |zEl |2) calculó el valor absoluto (El |zEl |) evitando el desbordamiento, y luego cuadró el resultado. Este problema podría solucionarse mediante el indicador de compilación -ffast-math.
  2. El compilador intel icc en linux (o enlazador) compilado std::arga un no-opt bajo ciertas configuraciones (compatibilidad de enlace con una versión gcc específica). El problema resurgió con demasiada frecuencia, por lo que std::argtuvo que ser reemplazado por atan2(imag(),real()). Pero fue demasiado fácil olvidar esto al escribir código nuevo.
  3. El tipo std::complexutiliza diferentes convenciones de llamada (= ABI) que el tipo complejo C99 integrado y el tipo complejo Fortran incorporado para las nuevas versiones de gcc.
  4. El -ffast-mathindicador de compilación interactúa con el manejo de excepciones de punto flotante de maneras inesperadas. Lo que sucede es que el compilador extrae las divisiones de los bucles, lo que provoca division by zeroexcepciones en tiempo de ejecución. Estas excepciones nunca habrían ocurrido dentro del ciclo, porque la división correspondiente no tuvo lugar debido a la lógica circundante. Esa fue realmente mala, porque era una biblioteca que se compiló por separado del programa que usaba la entrega de excepciones de punto flotante (usando diferentes indicadores de compilación) y se topaba con estos problemas (los equipos correspondientes estaban sentados en partes opuestas del mundo, así que Este problema realmente causó problemas graves). Esto se resolvió haciendo la optimización utilizada por el compilador a mano con más cuidado.
  5. La biblioteca se convirtió en parte del programa y ya no utilizó el -ffast-mathindicador de compilación. Después de una actualización a una versión más nueva de gcc, el rendimiento se redujo en un factor enorme. Todavía no he investigado este problema en detalle, pero me temo que está relacionado con el Anexo G del C99 . Tengo que admitir que estoy completamente confundido por esta extraña definición de multiplicación para números complejos, e incluso parece que existen diferentes versiones de esto con afirmaciones de que las otras versiones están equivocadas. Espero que el -fcx-limited-rangeindicador de compilación resuelva el problema, porque parece haber otro problema relacionado con -ffast-mathesta nueva versión de gcc.
  6. El -ffast-mathindicador de compilación hace que el comportamiento sea NaNcompletamente impredecible para las versiones más nuevas de gcc (incluso isnanse ve afectado). La única solución parece ser evitar cualquier aparición NaNen el programa, lo que anula el propósito de la existencia de NaN.

Ahora puede preguntar si planeo abandonar los tipos complejos integrados y std::complexpor estas razones. Me quedaré con los tipos integrados, siempre y cuando me quede con C ++. En caso de que C ++ logre volverse completamente inutilizable para la informática científica, preferiría considerar cambiar a un lenguaje que se ocupe más de los problemas relevantes para la informática científica.

Thomas Klimpel
fuente
Parece que mis temores relacionados con el Anexo G de C99 se han hecho realidad, y ahora -fcx-limited-range se requiere para una velocidad de cálculo decente cuando se multiplican números complejos. Al menos eso es lo que obtengo de la siguiente historia de guerra reciente: medium.com/@smcallis_71148/…
Thomas Klimpel