¿Cuál es el costo computacional de

Uno de los principales problemas con los que tenemos que lidiar en las simulaciones moleculares es el cálculo de las fuerzas dependientes de la distancia. Si podemos restringir las funciones de fuerza y ​​distancia para tener potencias iguales de la distancia de separación , entonces podemos calcular el cuadrado de la distancia r 2 = r ⋅ r y no tener que preocuparnos por r . Sin embargo, si hay potencias impares, entonces debemos tratar con r = $r$$r^2 = {\bf r \cdot r}$$r$ .$r = \sqrt{r^2}$

Mi pregunta es: ¿qué tan caro es la computación como se implementa en las bibliotecas de lenguajes comunes (C / C ++, Fortran, Python), etc. ¿Realmente se pueden obtener muchas mejoras de rendimiento ajustando manualmente el código para arquitecturas específicas?$\sqrt{x}$

Como extensión de la respuesta de Moyner , el chip sqrtsuele ser una rsqrtraíz cuadrada recíproca que calcula . Entonces, si en su código solo va a usar1/r(si está haciendo dinámica molecular, sí), puede calculardirectamente y ahorrarse la división. La razón por la cualse calcula en lugar dees que su iteración de Newton no tiene divisiones, solo adiciones y multiplicaciones.$a \rightarrow 1/\sqrt{a}$$1/r$r = rsqrt(r2)rsqrtsqrt

Como nota al margen, las divisiones también se calculan de forma iterativa y son casi tan lentas como rsqrten el hardware. Si busca eficiencia, es mejor que intente eliminar divisiones superfluas.

Algunas arquitecturas más modernas, como las arquitecturas POWER de IBM, no proporcionan rsqrtper se, sino una estimación precisa de unos pocos bits, por ejemplo, FRSQRTE . Cuando un usuario llama rsqrt, esto genera una estimación y luego una o dos (tantas como sea necesario) iteraciones del algoritmo de Newton o Goldschmidt usando multiplicaciones y sumas regulares. La ventaja de este enfoque es que los pasos de iteración se pueden canalizar e intercalar con otras instrucciones sin bloquear la FPU (para obtener una visión general muy buena de este concepto, aunque en arquitecturas antiguas, consulte la Tesis doctoral de Rolf Strebel ).

Para potenciales de interacción, la sqrtoperación se puede evitar por completo utilizando un interpolante polinómico de la función potencial, pero mi propio trabajo (implementado en mdcore) en esta área muestra que, al menos en arquitecturas de tipo x86, la sqrtinstrucción es lo suficientemente rápida.

Actualizar

Dado que esta respuesta parece estar recibiendo bastante atención, también me gustaría abordar la segunda parte de su pregunta, es decir, ¿realmente vale la pena intentar mejorar / eliminar operaciones básicas como sqrt?

En el contexto de las simulaciones de Molecular Dynamics, o cualquier simulación basada en partículas con interacciones limitadas por corte, hay mucho que ganar con mejores algoritmos para encontrar vecinos. Si está utilizando listas de Celdas , o algo similar, para encontrar vecinos o crear una lista Verlet , calculará una gran cantidad de distancias espurias por pares. En el caso ingenuo, solo el 16% de los pares de partículas inspeccionados estarán dentro de la distancia de corte entre sí. Aunque no se calcula la interacción para tales pares, acceder a los datos de partículas y calcular la distancia espuria por pares conlleva un gran costo.

Mi propio trabajo en esta área ( aquí , aquí y aquí ), así como el de otros (por ejemplo, aquí ), muestra cómo se pueden evitar estos cálculos espurios. Estos algoritmos de búsqueda de vecinos incluso superan a las listas de Verlet, como se describe aquí .

El punto que quiero enfatizar es que, aunque puede haber algunas mejoras para obtener de un mejor conocimiento / explotación de la arquitectura de hardware subyacente, también hay ganancias potencialmente mayores para repensar los algoritmos de nivel superior.

La raíz cuadrada se implementa en el hardware en la mayoría de los procesadores, es decir, hay instrucciones de ensamblaje específicas y el rendimiento debe ser comparable en la mayoría de los idiomas porque es muy difícil arruinar la implementación. Probablemente nunca podrá superar la instrucción FSQRT, ya que fue diseñada por un diseñador de hardware inteligente.

La forma en que se implementa en el hardware puede variar, pero probablemente sea algún tipo de iteración de punto fijo, por ejemplo, el método de Newton-Raphson que realiza un número específico de iteraciones hasta que se calcula el número de dígitos necesarios. Los métodos iterativos en hardware son en general mucho más lentos que otras operaciones, ya que se deben completar varios ciclos antes de que el resultado esté listo.

También hay algunas instrucciones SIMD de transmisión que se pueden usar en los registros XMM para cálculos vectoriales rápidos que se encuentran aquí . Estos registros son bastante pequeños, pero si tiene un número conocido de coordenadas (por ejemplo, un sistema de coordenadas cartesianas tridimensionales) pueden ser bastante más rápidas.

Si su idioma es lo suficientemente bajo, siempre puede escribir a una precisión más baja o usar un número de precisión más bajo para sus coordenadas. La precisión simple a menudo es más que suficiente, y por lo que recuerdo será más rápido al calcular raíces cuadradas, ya que las iteraciones pueden terminarse antes.

Debería ser bastante fácil comparar diferentes idiomas: simplemente escriba una larga serie de números aleatorios en un archivo, cárguelo usando diferentes idiomas y luego cronometre las raíces cuadradas.

Puede haber mejoras en el rendimiento, pero primero debe perfilarse para saber que calcular el recíproco del sqrt es el cuello de botella (y no, por ejemplo, cargar las posiciones y salvar las fuerzas).

El proyecto GROMACS MD surgió de una idea para explotar los detalles del formato de punto flotante IEEE para sembrar un esquema de iteración Newton-Raphson para calcular una aproximación aceptable al recíproco de la raíz cuadrada (ver Apéndice B.3 de http: / /www.gromacs.org/Documentation/Manual ), pero no hay CPU HPC en uso donde GROMACS todavía usa esta idea.