diff options
| -rw-r--r-- | CMakeLists.txt | 9 | ||||
| -rw-r--r-- | src/lib/CMakeLists.txt | 5 | ||||
| -rw-r--r-- | src/lib/omp.h | 63 | ||||
| -rw-r--r-- | src/lib/openmp_stubs.c | 1402 | ||||
| -rw-r--r-- | src/lib/tng_io.c | 4 | ||||
| -rw-r--r-- | src/tests/CMakeLists.txt | 9 | ||||
| -rw-r--r-- | src/tests/md_openmp.c | 880 | ||||
| -rw-r--r-- | src/tests/tng_io_testing.c | 2 |
8 files changed, 2370 insertions, 4 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index c459d69..c36dcb8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,6 +6,15 @@ project(TRAJECTORY) set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) +option(USE_OPENMP "Try to use the OpenMP library (if available)" ON) +if(USE_OPENMP) + FIND_PACKAGE(OpenMP) + if(OPENMP_FOUND) + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}") + endif() +endif() + add_subdirectory(src) file(COPY example_files DESTINATION .) diff --git a/src/lib/CMakeLists.txt b/src/lib/CMakeLists.txt index b12bb42..dda277d 100644 --- a/src/lib/CMakeLists.txt +++ b/src/lib/CMakeLists.txt @@ -1,2 +1,7 @@ +if(NOT OPENMP_FOUND) + add_library(openmp SHARED openmp_stubs.c) + set(OpenMP_LIBS "openmp") +endif() + add_library(tng_io SHARED tng_io.c md5.c) diff --git a/src/lib/omp.h b/src/lib/omp.h new file mode 100644 index 0000000..179797f --- /dev/null +++ b/src/lib/omp.h @@ -0,0 +1,63 @@ +#ifndef _OMP_H_DEF +#define _OMP_H_DEF +/* + Define the lock data types +*/ +typedef void *omp_lock_t; +typedef void *omp_nest_lock_t; +/* + Define the schedule kinds. + Add vendor specific schedule constants at end. +*/ +typedef enum omp_sched_t +{ + omp_sched_static = 1, + omp_sched_dynamic = 2, + omp_sched_guided = 3, + omp_sched_auto = 4 +} omp_sched_t; +/* + Exported OpenMP functions +*/ +#ifdef __cplusplus +extern "C" +{ +#endif + +extern void omp_destroy_lock ( omp_lock_t *lock ); +extern void omp_destroy_nest_lock ( omp_nest_lock_t *lock ); +extern int omp_get_active_level ( void ); +extern int omp_get_ancestor_thread_num ( int level ); +extern int omp_get_dynamic ( void ); +extern int omp_get_level ( void ); +extern int omp_get_max_active_levels ( void ); +extern int omp_get_max_threads ( void ); +extern int omp_get_nested ( void ); +extern int omp_get_num_procs ( void ); +extern int omp_get_num_threads ( void ); +extern void omp_get_schedule ( omp_sched_t *kind, int *modifier ); +extern int omp_get_team_size ( int level ); +extern int omp_get_thread_limit ( void ); +extern int omp_get_thread_num ( void ); +extern double omp_get_wtick ( void ); +extern double omp_get_wtime ( void ); +extern int omp_in_final ( void ); +extern int omp_in_parallel ( void ); +extern void omp_init_lock ( omp_lock_t *lock ); +extern void omp_init_nest_lock ( omp_nest_lock_t *lock ); +extern void omp_set_dynamic ( int dynamic_threads ); +extern void omp_set_lock ( omp_lock_t *lock ); +extern void omp_set_max_active_levels ( int max_active_levels ); +extern void omp_set_nest_lock ( omp_nest_lock_t *lock ); +extern void omp_set_nested ( int nested ); +extern void omp_set_num_threads ( int num_threads ); +extern void omp_set_schedule ( omp_sched_t kind, int modifier ); +extern int omp_test_lock ( omp_lock_t *lock ); +extern int omp_test_nest_lock ( omp_nest_lock_t *lock ); +extern void omp_unset_lock ( omp_lock_t *lock ); +extern void omp_unset_nest_lock ( omp_nest_lock_t *lock ); + +#ifdef __cplusplus +} +#endif +#endif
\ No newline at end of file diff --git a/src/lib/openmp_stubs.c b/src/lib/openmp_stubs.c new file mode 100644 index 0000000..5ce31dc --- /dev/null +++ b/src/lib/openmp_stubs.c @@ -0,0 +1,1402 @@ +# include <stdio.h> +# include <stdlib.h> +# include <time.h> + +# include "omp.h" + +struct __omp_lock +{ + int lock; +}; + +enum { UNLOCKED = -1, INIT, LOCKED }; + +struct __omp_nest_lock +{ + short owner; + short count; +}; + +enum { NOOWNER = -1, MASTER = 0 }; + +/******************************************************************************/ + +void omp_destroy_lock ( omp_lock_t *arg ) + +/******************************************************************************/ +/* + Purpose: + + OMP_DESTROY_LOCK destroys a simple lock. + + Discussion: + + The routine is intended to return the state of the lock to the + uninitialized state. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, omp_lock_t *ARG, the simple lock. +*/ +{ + struct __omp_lock *lock = ( struct __omp_lock * ) arg; + lock->lock = INIT; + + return; +} +/******************************************************************************/ + +void omp_destroy_nest_lock ( omp_nest_lock_t *arg ) + +/******************************************************************************/ +/* + Purpose: + + OMP_DESTROY_NEST_LOCK destroys a nestable lock. + + Discussion: + + The routine is intended to return the state of the lock to the + uninitialized state. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, omp_nest_lock_t *ARG, the nestable lock. +*/ +{ + struct __omp_nest_lock *nlock = ( struct __omp_nest_lock * ) arg; + nlock->owner = NOOWNER; + nlock->count = UNLOCKED; + return; +} +/******************************************************************************/ + +int omp_get_active_level ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_ACTIVE_LEVEL returns the number of nested active parallel regions. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_GET_ACTIVE_LEVEL, the number of nested active parallel + regions enclosing the task that contains this call. +*/ +{ + return 0; +} +/******************************************************************************/ + +int omp_get_ancestor_thread_num ( int level ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_ANCESTOR_THREAD_NUM returns the thread number of the ancestor of this thread. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input, int LEVEL, the nested level. + + Output, int OMP_GET_ANCESTOR_THREAD_NUM, the thread number of the + ancestor of this thread. +*/ +{ + if ( level == 0 ) + { + return 0; + } + else + { + return -1; + } +} +/******************************************************************************/ + +int omp_get_dynamic ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_DYNAMIC reports if dynamic adjustment of thread number is allowed. + + Discussion: + + The user can request dynamic thread adjustment by calling OMP_SET_DYNAMIC. + + For this stub library, the value FALSE is always returned. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_GET_DYNAMIC, is TRUE (1) if dynamic adjustment of thread + number has been enable, by default or by a user call, and FALSE (0) + otherwise. +*/ +{ + return 0; +} +/******************************************************************************/ + +int omp_get_level ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_LEVEL returns the number of nested parallel regions enclosing this task. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_GET_LEVEL, the number of nested parallel regions + enclosing this task. +*/ +{ + return 0; +} +/******************************************************************************/ + +int omp_get_max_active_levels ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_MAX_ACTIVE_LEVELS gets the maximum number of nested active parallel regions. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_GET_MAX_ACTIVE_LEVELS gets the maximum number of + nested active parallel regions. +*/ +{ + return 0; +} +/******************************************************************************/ + +int omp_get_max_threads ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_MAX_THREADS returns the default number of threads. + + Discussion: + + If a parallel region is reached, and no number of threads has been + specified explicitly, there is a default number of threads that will + be used to form the new team. That value is returned by this function. + + For this stub library, the value 1 is always returned. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_GET_MAX_THREADS, the default number of threads. +*/ +{ + return 1; +} +/******************************************************************************/ + +int omp_get_nested ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_NESTED reports if nested parallelism has been enabled. + + Discussion: + + The user can request nested parallelism by calling OMP_SET_NESTED. + + For this stub library, the value FALSE is always returned. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_GET_NESTED, is TRUE (1) if nested parallelism has been + enable by default or by a user call, and FALSE (0) otherwise. +*/ +{ + return 0; +} +/******************************************************************************/ + +int omp_get_num_procs ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_NUM_PROCS returns the number of processors available to the program. + + Discussion: + + For this stub library, the value 1 is always returned. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_GET_NUM_PROCS, the number of processors available. +*/ +{ + return 1; +} +/******************************************************************************/ + +int omp_get_num_threads ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_NUM_THREADS returns the number of threads in the current team. + + Discussion: + + For this stub library, the value 1 is always returned. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_GET_NUM_THREADS, the number of threads in the + current team. +*/ +{ + return 1; +} +/******************************************************************************/ + +void omp_get_schedule ( omp_sched_t *kind, int *modifier ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_SCHEDULE returns information about the "runtime" schedule. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, omp_sched_t *KIND, may be + 1, omp_sched_static, + 2, omp_sched_dynamic, + 3, omp_sched_guided, + 4, omp_sched_auto. + + Output, int *MODIFIER; this contains the "chunk_size" information for + static, dynamic, or guided schedules, and is ignored for the auto schedule. + If the chunk_size is less than 1, then the default value is used instead. +*/ +{ + *kind = omp_sched_static; + *modifier = 0; + return; +} +/******************************************************************************/ + +int omp_get_team_size ( int level ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_TEAM_SIZE returns the size of the thread team for a given level. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input, int LEVEL, the nested level. + + Output, int OMP_GET_TEAM_SIZE, the size of the thread team for + this level. +*/ +{ + if ( level == 0 ) + { + return 1; + } + else + { + return -1; + } +} +/******************************************************************************/ + +int omp_get_thread_limit ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_THREAD_LIMIT returns the maximum number of OpenMP threads available. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_GET_THREAD_LIMIT, the maximum number of OpenMP + threads available. +*/ +{ + return 1; +} +/******************************************************************************/ + +int omp_get_thread_num ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_THREAD_NUM returns the thread number of a thread in a team. + + Discussion: + + Thread numbers start at 0. + + If this function is not called from a parallel region, then only one + thread is executing, so the value returned is 0. + + For this stub library, the value 0 is always returned. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_GET_THREAD_NUM, the thread number. +*/ +{ + return 0; +} +/******************************************************************************/ + +double omp_get_wtick ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_WTICK returns the precision of the timer used by OMP_GET_WTIME. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, double OMP_GET_WTICK, the number of seconds between + successive "ticks" of the wall clock timer. +*/ +{ + double value; + + value = 1.0 / ( double ) CLOCKS_PER_SEC; + + return value; +} +/******************************************************************************/ + +double omp_get_wtime ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_GET_WTIME returns elapsed wall clock time in seconds. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, double OMP_GET_WTIME, the current reading of the + wall clock timer. +*/ +{ + double value; + + value = ( double ) clock ( ) / ( double ) CLOCKS_PER_SEC; + + return value; +} +/******************************************************************************/ + +int omp_in_final ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_IN_FINAL is true if the routine is executed in a final task region. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_IN_FINAL, is true if the routine is executed in a + final task region. +*/ +{ + return 1; +} +/******************************************************************************/ + +int omp_in_parallel ( void ) + +/******************************************************************************/ +/* + Purpose: + + OMP_IN_PARALLEL returns TRUE if the call is made from a parallel region. + + Discussion: + + For this stub library, the value FALSE is always returned. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, int OMP_IN_PARALLEL, is "TRUE" (1) if the routine was called + from a parallel region and "FALSE" (0) otherwise. +*/ +{ + return 0; +} +/******************************************************************************/ + +void omp_init_lock ( omp_lock_t *arg ) + +/******************************************************************************/ +/* + Purpose: + + OMP_INIT_LOCK initializes a simple lock. + + Discussion: + + This routine is intended to initialize the lock to the unlocked state. + + For this stub library, the lock points to -1. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, omp_lock_t *ARG, the simple lock. +*/ +{ + struct __omp_lock *lock = ( struct __omp_lock * ) arg; + lock->lock = UNLOCKED; + return; +} +/******************************************************************************/ + +void omp_init_nest_lock ( omp_nest_lock_t *arg ) + +/******************************************************************************/ +/* + Purpose: + + OMP_INIT_NEST_LOCK initializes a nestable lock. + + Discussion: + + This routine is intended to initialize the lock to the unlocked state. + + For this stub library, the lock points to -1. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Output, omp_nest_lock_t *ARG, the lock. +*/ +{ + struct __omp_nest_lock *nlock = ( struct __omp_nest_lock * ) arg; + nlock->owner = NOOWNER; + nlock->count = 0; + return; +} +/******************************************************************************/ + +void omp_set_dynamic ( int dynamic_threads ) + +/******************************************************************************/ +/* + Purpose: + + OMP_SET_DYNAMIC enables dynamic adjustment of the number of threads. + + Discussion: + + If DYNAMIC_THREADS is TRUE, then the number of threads available for + execution in a parallel region may be dynamically adjusted. + + For this stub library, the input value is ignored. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input, int DYNAMIC_THREADS, is TRUE if the user wishes to allow + dynamic adjustment of the number of threads available for execution + in any parallel region. +*/ +{ + return; +} +/******************************************************************************/ + +void omp_set_lock ( omp_lock_t *arg ) + +/******************************************************************************/ +/* + Purpose: + + OMP_SET_LOCK sets a simple lock. + + Discussion: + + The lock must already have been initialized. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input/output, omp_lock_t *ARG, the simple lock. +*/ +{ + struct __omp_lock *lock = ( struct __omp_lock * ) arg; + if ( lock->lock == UNLOCKED ) + { + lock->lock = LOCKED; + } + else if ( lock->lock == LOCKED ) + { + fprintf ( stderr, "error: deadlock in using lock variable\n" ); + exit ( 1 ); + } + else + { + fprintf ( stderr, "error: lock not initialized\n" ); + exit ( 1 ); + } + return; +} +/******************************************************************************/ + +void omp_set_max_active_levels ( int max_active_levels ) + +/******************************************************************************/ +/* + Purpose: + + OMP_SET_MAX_ACTIVE_LEVELS limits the number of nested active parallel regions. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input, int MAX_LEVELS, the maximum number of nested active parallel + regions. +*/ +{ + return; +} +/******************************************************************************/ + +void omp_set_nest_lock ( omp_nest_lock_t *arg ) + +/******************************************************************************/ +/* + Purpose: + + OMP_SET_NEST_LOCK sets a nestable lock. + + Discussion: + + The lock must already have been initialized. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input/output, omp_nest_lock_t *ARG, the nestable lock. +*/ +{ + struct __omp_nest_lock *nlock = ( struct __omp_nest_lock * ) arg; + if ( nlock->owner == MASTER && nlock->count >= 1 ) + { + nlock->count++; + } + else if ( nlock->owner == NOOWNER && nlock->count == 0 ) + { + nlock->owner = MASTER; + nlock->count = 1; + } + else + { + fprintf ( stderr, "error: lock corrupted or not initialized\n" ); + exit ( 1 ); + } + return; +} +/******************************************************************************/ + +void omp_set_nested ( int nested ) + +/******************************************************************************/ +/* + Purpose: + + OMP_SET_NESTED controls the use of nested parallelism. + + Discussion: + + For this stub library, the input value is ignored. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input, int NESTED, is TRUE (1) if nested parallelism is to be enabled. +*/ +{ + return; +} +/******************************************************************************/ + +void omp_set_num_threads ( int num_threads ) + +/******************************************************************************/ +/* + Purpose: + + OMP_SET_NUM_THREADS sets the number of threads. + + Discussion: + + This routine sets the number of threads to be used in all subsequent + parallel regions for which an explicit number of threads is not + specified. + + For this stub library, the input value is ignored. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input, int NUM_THREADS, the number of threads to be used in all + subsequent parallel regions for which an explicit number of threads + is not specified. 0 < NUM_THREADS. +*/ +{ + return; +} +/******************************************************************************/ + +void omp_set_schedule ( omp_sched_t kind, int modifier ) + +/******************************************************************************/ +/* + Purpose: + + OMP_SET_SCHEDULE chooses the schedule when "runtime" is the schedule kind. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input, omp_sched_t KIND, may be + 1, omp_sched_static, + 2, omp_sched_dynamic, + 3, omp_sched_guided, + 4, omp_sched_auto. + + Input, int MODIFIER; this contains the "chunk_size" information for + static, dynamic, or guided schedules, and is ignored for the auto schedule. + If the chunk_size is less than 1, then the default value is used instead. +*/ +{ + return; +} +/******************************************************************************/ + +int omp_test_lock ( omp_lock_t *arg ) + +/******************************************************************************/ +/* + Purpose: + + OMP_TEST_LOCK tests a simple lock. + + Discussion: + + Calling this routine with an uninitialized lock causes an error. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input/output, omp_lock_t *ARG, the simple lock. + If the lock was initialized but not set, it is set by this call. + + Output, int OMP_TEST_LOCK: + TRUE (1), on input, the lock was initialized and not set; + FALSE (0), on input the lock was initialized, and set. +*/ +{ + struct __omp_lock *lock = ( struct __omp_lock * ) arg; + if ( lock->lock == UNLOCKED ) + { + lock->lock = LOCKED; + return 1; + } + else if ( lock->lock == LOCKED ) + { + return 0; + } + else + { + fprintf ( stderr, "error: lock not initialized\n" ); + exit ( 1 ); + } +} +/******************************************************************************/ + +int omp_test_nest_lock ( omp_nest_lock_t *arg ) + +/******************************************************************************/ +/* + Purpose: + + OMP_TEST_NEST_LOCK tests a nestable lock. + + Discussion: + + Calling this routine with an uninitialized lock causes an error. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input/output, omp_nest_lock_t *ARG, the nestable lock. + If the lock was initialized but not set, it is set by this call. + + Output, int OMP_TEST_NEST_LOCK, returns the new nesting count, + if the call was successful. Otherwise, the value 0 is returned. +*/ +{ + struct __omp_nest_lock *nlock = ( struct __omp_nest_lock * ) arg; + omp_set_nest_lock ( arg) ; + return nlock->count; +} +/******************************************************************************/ + +void omp_unset_lock ( omp_lock_t *arg ) + +/******************************************************************************/ +/* + Purpose: + + OMP_UNSET_LOCK unsets a simple lock. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input/output, omp_lock_t *ARG, the simple lock. +*/ +{ + struct __omp_lock *lock = ( struct __omp_lock * ) arg; + if ( lock->lock == LOCKED ) + { + lock->lock = UNLOCKED; + } + else if ( lock->lock == UNLOCKED ) + { + fprintf ( stderr, "error: lock not set\n" ); + exit ( 1 ); + } + else + { + fprintf ( stderr, "error: lock not initialized\n" ); + exit ( 1 ); + } + return; +} +/******************************************************************************/ + +void omp_unset_nest_lock ( omp_nest_lock_t *arg ) + +/******************************************************************************/ +/* + Purpose: + + OMP_UNSET_NEST_LOCK unsets a nestable lock. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 May 2012 + + Author: + + John Burkardt + + Reference: + + OpenMP Application Program Interface, + Version 3.1, + July 2011. + + Parameters: + + Input/output, omp_nest_lock_t *ARG, the nestable lock. +*/ +{ + struct __omp_nest_lock *nlock = ( struct __omp_nest_lock * ) arg; + if ( nlock->owner == MASTER && nlock->count >= 1 ) + { + nlock->count--; + if ( nlock->count == 0 ) + { + nlock->owner = NOOWNER; + } + } + else if ( nlock->owner == NOOWNER && nlock->count == 0 ) + { + fprintf ( stderr, "error: lock not set\n" ); + exit ( 1 ); + } + else + { + fprintf ( stderr, "error: lock corrupted or not initialized\n" ); + exit ( 1 ); + } + return; +} diff --git a/src/lib/tng_io.c b/src/lib/tng_io.c index c97a3fb..12b7a4e 100644 --- a/src/lib/tng_io.c +++ b/src/lib/tng_io.c @@ -8275,7 +8275,7 @@ tng_function_status tng_frame_particle_data_write(tng_trajectory_t tng_data, memcpy(copy, values, val_n_particles * n_values_per_frame * size); for(i = 0; i < val_n_particles * n_values_per_frame; i++) { - tng_swap_byte_order_64(tng_data, ©[i*size]); + tng_swap_byte_order_64(tng_data, (int64_t *)copy+i*size); } fwrite(copy, val_n_particles * n_values_per_frame, size, tng_data->output_file); @@ -8288,7 +8288,7 @@ tng_function_status tng_frame_particle_data_write(tng_trajectory_t tng_data, memcpy(copy, values, val_n_particles * n_values_per_frame * size); for(i = 0; i < val_n_particles * n_values_per_frame; i++) { - tng_swap_byte_order_32(tng_data, ©[i*size]); + tng_swap_byte_order_32(tng_data, (int32_t *)copy+i*size); } fwrite(copy, val_n_particles * n_values_per_frame, size, tng_data->output_file); diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 6747eb6..92c580d 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -1,3 +1,7 @@ +if(NOT OPENMP_FOUND) + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") +endif() + include_directories(${TRAJECTORY_SOURCE_DIR}/src/lib) link_directories(${TRAJECTORY_BINARY_DIR}/src/lib) @@ -7,4 +11,7 @@ configure_file(tng_io_testing.h.in ${CMAKE_BINARY_DIR}/generated/tng_io_testing. include_directories(${CMAKE_BINARY_DIR}/generated/) add_executable(tng_testing tng_io_testing.c) -target_link_libraries(tng_testing tng_io)
\ No newline at end of file +target_link_libraries(tng_testing tng_io) + +add_executable(md_openmp md_openmp.c) +target_link_libraries(md_openmp tng_io ${OpenMP_LIBS} m)
\ No newline at end of file diff --git a/src/tests/md_openmp.c b/src/tests/md_openmp.c new file mode 100644 index 0000000..0b5241a --- /dev/null +++ b/src/tests/md_openmp.c @@ -0,0 +1,880 @@ +# include <stdlib.h> +# include <stdio.h> +# include <time.h> +# include <math.h> +# include <omp.h> +# include <tng_io.h> +#include "tng_io_testing.h" + +int main ( int argc, char *argv[] ); +void compute ( int np, int nd, double pos[], double vel[], + double mass, double f[], double *pot, double *kin ); +double dist ( int nd, double r1[], double r2[], double dr[] ); +void initialize ( int np, int nd, double box[], int *seed, double pos[], + double vel[], double acc[] ); +double r8_uniform_01 ( int *seed ); +void timestamp ( void ); +void update ( int np, int nd, double pos[], double vel[], double f[], + double acc[], double mass, double dt ); + +/******************************************************************************/ + +int main ( int argc, char *argv[] ) + +/******************************************************************************/ +/* + Purpose: + + MAIN is the main program for MD_OPENMP. + + Discussion: + + MD implements a simple molecular dynamics simulation. + + The program uses Open MP directives to allow parallel computation. + + The velocity Verlet time integration scheme is used. + + The particles interact with a central pair potential. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 30 July 2009 + + Author: + + Original FORTRAN77 version by Bill Magro. + C version by John Burkardt. + TNG trajectory output by Magnus Lundborg. + + Parameters: + + None +*/ +{ + double *acc; + double *box; + double *box_shape; + double dt = 0.0001; + double e0; + double *force; + int i; + double kinetic; + double mass = 1.0; + int nd = 3; + int np = 100; + double *pos; + double potential; + int proc_num; + int seed = 123456789; + int step; + int step_num = 1000; + int step_print; + int step_print_index; + int step_print_num; + int step_save; + double *vel; + double wtime; + tng_trajectory_t traj; + tng_molecule_t molecule; + tng_chain_t chain; + tng_residue_t residue; + tng_atom_t atom; + int64_t n_frames_per_frame_set; + int frames_saved_cnt = 0; + int frame_set_cnt = 0; + double **data; + + timestamp ( ); + + proc_num = omp_get_num_procs ( ); + + acc = ( double * ) malloc ( nd * np * sizeof ( double ) ); + box = ( double * ) malloc ( nd * sizeof ( double ) ); + box_shape = (double *) malloc (9 * sizeof (double)); + force = ( double * ) malloc ( nd * np * sizeof ( double ) ); + pos = ( double * ) malloc ( nd * np * sizeof ( double ) ); + vel = ( double * ) malloc ( nd * np * sizeof ( double ) ); + + printf ( "\n" ); + printf ( "MD_OPENMP\n" ); + printf ( " C/OpenMP version\n" ); + printf ( "\n" ); + printf ( " A molecular dynamics program.\n" ); + + printf ( "\n" ); + printf ( " NP, the number of particles in the simulation is %d\n", np ); + printf ( " STEP_NUM, the number of time steps, is %d\n", step_num ); + printf ( " DT, the size of each time step, is %f\n", dt ); + + printf ( "\n" ); + printf ( " Number of processors available = %d\n", omp_get_num_procs ( ) ); + printf ( " Number of threads = %d\n", omp_get_max_threads ( ) ); + + + printf("\n"); + printf(" Initializing trajectory storage.\n"); + if(tng_trajectory_init(&traj) != TNG_SUCCESS) + { + tng_trajectory_destroy(&traj); + printf(" Cannot init trajectory.\n"); + exit(1); + } +#ifdef EXAMPLE_FILES_DIR + tng_output_file_set(traj, EXAMPLE_FILES_DIR "tng_md_out.tng"); +#else + tng_output_file_set(traj, "/tmp/tng_md_out.tng"); +#endif + + + + /* Set molecules data */ + printf(" Creating molecules in trajectory.\n"); + tng_molecule_add(traj, "water", &molecule); + tng_molecule_chain_add(traj, molecule, "W", &chain); + tng_chain_residue_add(traj, chain, "WAT", &residue); + if(tng_residue_atom_add(traj, residue, "O", "O", &atom) == TNG_CRITICAL) + { + tng_trajectory_destroy(&traj); + printf(" Cannot create molecules.\n"); + exit(1); + } + tng_molecule_cnt_set(traj, molecule, np); + + +/* + Set the dimensions of the box. +*/ + for(i = 0; i < 9; i++) + { + box_shape[i] = 0.0; + } + for ( i = 0; i < nd; i++ ) + { + box[i] = 10.0; + box_shape[i*nd + i] = box[i]; + } + + + if(tng_data_block_add(traj, TNG_TRAJ_BOX_SHAPE, "BOX SHAPE", TNG_DOUBLE_DATA, + TNG_NON_TRAJECTORY_BLOCK, 1, 9, 1, TNG_UNCOMPRESSED, + box_shape) == TNG_CRITICAL || + tng_file_headers_write(traj, TNG_USE_HASH) == TNG_CRITICAL) + { + tng_trajectory_destroy(&traj); + printf(" Cannot write trajectory headers and box shape.\n"); + exit(1); + } + + printf ( "\n" ); + printf ( " Initializing positions, velocities, and accelerations.\n" ); +/* + Set initial positions, velocities, and accelerations. +*/ + initialize ( np, nd, box, &seed, pos, vel, acc ); +/* + Compute the forces and energies. +*/ + printf ( "\n" ); + printf ( " Computing initial forces and energies.\n" ); + + compute ( np, nd, pos, vel, mass, force, &potential, &kinetic ); + + e0 = potential + kinetic; +/* + This is the main time stepping loop: + Compute forces and energies, + Update positions, velocities, accelerations. +*/ + printf ( "\n" ); + printf ( " At each step, we report the potential and kinetic energies.\n" ); + printf ( " The sum of these energies should be a constant.\n" ); + printf ( " As an accuracy check, we also print the relative error\n" ); + printf ( " in the total energy.\n" ); + printf ( "\n" ); + printf ( " Step Potential Kinetic (P+K-E0)/E0\n" ); + printf ( " Energy P Energy K Relative Energy Error\n" ); + printf ( "\n" ); + + step_print = 0; + step_print_index = 0; + step_print_num = 10; + + step = 0; + printf ( " %8d %14f %14f %14e\n", + step, potential, kinetic, ( potential + kinetic - e0 ) / e0 ); + step_print_index++; + step_print = ( step_print_index * step_num ) / step_print_num; + + /* Saving frequency */ + step_save = 5; + + wtime = omp_get_wtime ( ); + + /* Create a frame set for writing data */ + tng_num_frames_per_frame_set_get(traj, &n_frames_per_frame_set); + if(tng_frame_set_new(traj, 0, + n_frames_per_frame_set) != TNG_SUCCESS) + { + printf("Error creating frame set %d. %s: %d\n", + i, __FILE__, __LINE__); + exit(1); + } + frame_set_cnt++; + + /* Setup an empty data array - this will be written to the frame set + * (since the block size needs to be correct in the file. After that it + * is filled with data. */ + data = malloc(n_frames_per_frame_set * 3 * np * sizeof(double)); + + for(i = 0; i < n_frames_per_frame_set * 3 * np; i++) + { + data[i] = 0; + } + + + /* Add data blocks */ + if(tng_particle_data_block_add(traj, TNG_TRAJ_POSITIONS, + "POSITIONS", + TNG_DOUBLE_DATA, + TNG_TRAJECTORY_BLOCK, + n_frames_per_frame_set, 3, + 1, 0, np, + TNG_UNCOMPRESSED, + data) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + if(tng_particle_data_block_add(traj, TNG_TRAJ_VELOCITIES, + "VELOCITIES", + TNG_DOUBLE_DATA, + TNG_TRAJECTORY_BLOCK, + n_frames_per_frame_set, 3, + 1, 0, np, + TNG_UNCOMPRESSED, + data) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + if(tng_particle_data_block_add(traj, TNG_TRAJ_FORCES, + "FORCES", + TNG_DOUBLE_DATA, + TNG_TRAJECTORY_BLOCK, + n_frames_per_frame_set, 3, + 1, 0, np, + TNG_UNCOMPRESSED, + data) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + + /* Write the frame set to disk */ + if(tng_frame_set_write(traj, TNG_SKIP_HASH) != TNG_SUCCESS) + { + printf("Error writing frame set. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + + if(tng_frame_particle_data_write(traj, 0, + TNG_TRAJ_POSITIONS, 0, np, + pos, TNG_USE_HASH) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + if(tng_frame_particle_data_write(traj, 0, + TNG_TRAJ_VELOCITIES, 0, np, + vel, TNG_USE_HASH) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + if(tng_frame_particle_data_write(traj, 0, + TNG_TRAJ_FORCES, 0, np, + force, TNG_USE_HASH) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + frames_saved_cnt++; + + +// printf("pos: %f, vel: %f, force: %f\n", pos[0], vel[0], force[0]); + + for ( step = 1; step <= step_num; step++ ) + { + compute ( np, nd, pos, vel, mass, force, &potential, &kinetic ); + + if ( step == step_print ) + { + printf ( " %8d %14f %14f %14e\n", step, potential, kinetic, + ( potential + kinetic - e0 ) / e0 ); + step_print_index++; + step_print = ( step_print_index * step_num ) / step_print_num; + } + if(step % step_save == 0) + { + if(frames_saved_cnt % n_frames_per_frame_set == 0) + { + if(tng_frame_set_new(traj, frames_saved_cnt, + n_frames_per_frame_set) != TNG_SUCCESS) + { + printf("Error creating frame set %d. %s: %d\n", + i, __FILE__, __LINE__); + free(data); + exit(1); + } + /* Add data blocks */ + if(tng_particle_data_block_add(traj, TNG_TRAJ_POSITIONS, + "POSITIONS", + TNG_DOUBLE_DATA, + TNG_TRAJECTORY_BLOCK, + n_frames_per_frame_set, 3, + 1, 0, np, + TNG_UNCOMPRESSED, + data) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + if(tng_particle_data_block_add(traj, TNG_TRAJ_VELOCITIES, + "VELOCITIES", + TNG_DOUBLE_DATA, + TNG_TRAJECTORY_BLOCK, + n_frames_per_frame_set, 3, + 1, 0, np, + TNG_UNCOMPRESSED, + data) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + if(tng_particle_data_block_add(traj, TNG_TRAJ_FORCES, + "FORCES", + TNG_DOUBLE_DATA, + TNG_TRAJECTORY_BLOCK, + n_frames_per_frame_set, 3, + 1, 0, np, + TNG_UNCOMPRESSED, + data) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + + /* Write the frame set to disk */ + if(tng_frame_set_write(traj, TNG_SKIP_HASH) != TNG_SUCCESS) + { + printf("Error writing frame set. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + + } + if(tng_frame_particle_data_write(traj, frames_saved_cnt, + TNG_TRAJ_POSITIONS, 0, np, + pos, TNG_USE_HASH) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + if(tng_frame_particle_data_write(traj, frames_saved_cnt, + TNG_TRAJ_VELOCITIES, 0, np, + vel, TNG_USE_HASH) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + if(tng_frame_particle_data_write(traj, frames_saved_cnt, + TNG_TRAJ_FORCES, 0, np, + force, TNG_USE_HASH) != TNG_SUCCESS) + { + printf("Error adding data. %s: %d\n", __FILE__, __LINE__); + free(data); + exit(1); + } + frames_saved_cnt++; + + } + update ( np, nd, pos, vel, force, acc, mass, dt ); + } + wtime = omp_get_wtime ( ) - wtime; + + printf ( "\n" ); + printf ( " Elapsed time for main computation:\n" ); + printf ( " %f seconds.\n", wtime ); + + free ( acc ); + free ( box ); + free ( force ); + free ( pos ); + free ( vel ); +/* + Terminate. +*/ + printf ( "\n" ); + printf ( "MD_OPENMP\n" ); + printf ( " Normal end of execution.\n" ); + + printf ( "\n" ); + timestamp ( ); + + return 0; +} +/******************************************************************************/ + +void compute ( int np, int nd, double pos[], double vel[], + double mass, double f[], double *pot, double *kin ) + +/******************************************************************************/ +/* + Purpose: + + COMPUTE computes the forces and energies. + + Discussion: + + The computation of forces and energies is fully parallel. + + The potential function V(X) is a harmonic well which smoothly + saturates to a maximum value at PI/2: + + v(x) = ( sin ( min ( x, PI2 ) ) )**2 + + The derivative of the potential is: + + dv(x) = 2.0 * sin ( min ( x, PI2 ) ) * cos ( min ( x, PI2 ) ) + = sin ( 2.0 * min ( x, PI2 ) ) + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 21 November 2007 + + Author: + + Original FORTRAN77 version by Bill Magro. + C version by John Burkardt. + + Parameters: + + Input, int NP, the number of particles. + + Input, int ND, the number of spatial dimensions. + + Input, double POS[ND*NP], the position of each particle. + + Input, double VEL[ND*NP], the velocity of each particle. + + Input, double MASS, the mass of each particle. + + Output, double F[ND*NP], the forces. + + Output, double *POT, the total potential energy. + + Output, double *KIN, the total kinetic energy. +*/ +{ + double d; + double d2; + int i; + int j; + int k; + double ke; + double pe; + double PI2 = 3.141592653589793 / 2.0; + double rij[3]; + + pe = 0.0; + ke = 0.0; + +# pragma omp parallel \ + shared ( f, nd, np, pos, vel ) \ + private ( i, j, k, rij, d, d2 ) + + +# pragma omp for reduction ( + : pe, ke ) + for ( k = 0; k < np; k++ ) + { +/* + Compute the potential energy and forces. +*/ + for ( i = 0; i < nd; i++ ) + { + f[i+k*nd] = 0.0; + } + + for ( j = 0; j < np; j++ ) + { + if ( k != j ) + { + d = dist ( nd, pos+k*nd, pos+j*nd, rij ); +/* + Attribute half of the potential energy to particle J. +*/ + if ( d < PI2 ) + { + d2 = d; + } + else + { + d2 = PI2; + } + + pe = pe + 0.5 * pow ( sin ( d2 ), 2 ); + + for ( i = 0; i < nd; i++ ) + { + f[i+k*nd] = f[i+k*nd] - rij[i] * sin ( 2.0 * d2 ) / d; + } + } + } +/* + Compute the kinetic energy. +*/ + for ( i = 0; i < nd; i++ ) + { + ke = ke + vel[i+k*nd] * vel[i+k*nd]; + } + } + + ke = ke * 0.5 * mass; + + *pot = pe; + *kin = ke; + + return; +} +/******************************************************************************/ + +double dist ( int nd, double r1[], double r2[], double dr[] ) + +/******************************************************************************/ +/* + Purpose: + + DIST computes the displacement (and its norm) between two particles. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 21 November 2007 + + Author: + + Original FORTRAN77 version by Bill Magro. + C version by John Burkardt. + + Parameters: + + Input, int ND, the number of spatial dimensions. + + Input, double R1[ND], R2[ND], the positions of the particles. + + Output, double DR[ND], the displacement vector. + + Output, double D, the Euclidean norm of the displacement. +*/ +{ + double d; + int i; + + d = 0.0; + for ( i = 0; i < nd; i++ ) + { + dr[i] = r1[i] - r2[i]; + d = d + dr[i] * dr[i]; + } + d = sqrt ( d ); + + return d; +} +/******************************************************************************/ + +void initialize ( int np, int nd, double box[], int *seed, double pos[], + double vel[], double acc[] ) + +/******************************************************************************/ +/* + Purpose: + + INITIALIZE initializes the positions, velocities, and accelerations. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 21 November 2007 + + Author: + + Original FORTRAN77 version by Bill Magro. + C version by John Burkardt. + + Parameters: + + Input, int NP, the number of particles. + + Input, int ND, the number of spatial dimensions. + + Input, double BOX[ND], specifies the maximum position + of particles in each dimension. + + Input, int *SEED, a seed for the random number generator. + + Output, double POS[ND*NP], the position of each particle. + + Output, double VEL[ND*NP], the velocity of each particle. + + Output, double ACC[ND*NP], the acceleration of each particle. +*/ +{ + int i; + int j; +/* + Give the particles random positions within the box. +*/ + for ( i = 0; i < nd; i++ ) + { + for ( j = 0; j < np; j++ ) + { + pos[i+j*nd] = box[i] * r8_uniform_01 ( seed ); + } + } + + for ( j = 0; j < np; j++ ) + { + for ( i = 0; i < nd; i++ ) + { + vel[i+j*nd] = 0.0; + } + } + for ( j = 0; j < np; j++ ) + { + for ( i = 0; i < nd; i++ ) + { + acc[i+j*nd] = 0.0; + } + } + return; +} +/******************************************************************************/ + +double r8_uniform_01 ( int *seed ) + +/******************************************************************************/ +/* + Purpose: + + R8_UNIFORM_01 is a unit pseudorandom R8. + + Discussion: + + This routine implements the recursion + + seed = 16807 * seed mod ( 2**31 - 1 ) + unif = seed / ( 2**31 - 1 ) + + The integer arithmetic never requires more than 32 bits, + including a sign bit. + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 11 August 2004 + + Author: + + John Burkardt + + Reference: + + Paul Bratley, Bennett Fox, Linus Schrage, + A Guide to Simulation, + Springer Verlag, pages 201-202, 1983. + + Bennett Fox, + Algorithm 647: + Implementation and Relative Efficiency of Quasirandom + Sequence Generators, + ACM Transactions on Mathematical Software, + Volume 12, Number 4, pages 362-376, 1986. + + Parameters: + + Input/output, int *SEED, a seed for the random number generator. + + Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between + 0 and 1. +*/ +{ + int k; + double r; + + k = *seed / 127773; + + *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; + + if ( *seed < 0 ) + { + *seed = *seed + 2147483647; + } + + r = ( double ) ( *seed ) * 4.656612875E-10; + + return r; +} +/******************************************************************************/ + +void timestamp ( void ) + +/******************************************************************************/ +/* + Purpose: + + TIMESTAMP prints the current YMDHMS date as a time stamp. + + Example: + + 31 May 2001 09:45:54 AM + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 24 September 2003 + + Author: + + John Burkardt + + Parameters: + + None +*/ +{ +# define TIME_SIZE 40 + + static char time_buffer[TIME_SIZE]; + const struct tm *tm; + size_t len; + time_t now; + + now = time ( NULL ); + tm = localtime ( &now ); + + len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); + + printf ( "%s\n", time_buffer ); + + return; +# undef TIME_SIZE +} +/******************************************************************************/ + +void update ( int np, int nd, double pos[], double vel[], double f[], + double acc[], double mass, double dt ) + +/******************************************************************************/ +/* + Purpose: + + UPDATE updates positions, velocities and accelerations. + + Discussion: + + The time integration is fully parallel. + + A velocity Verlet algorithm is used for the updating. + + x(t+dt) = x(t) + v(t) * dt + 0.5 * a(t) * dt * dt + v(t+dt) = v(t) + 0.5 * ( a(t) + a(t+dt) ) * dt + a(t+dt) = f(t) / m + + Licensing: + + This code is distributed under the GNU LGPL license. + + Modified: + + 17 April 2009 + + Author: + + Original FORTRAN77 version by Bill Magro. + C version by John Burkardt. + + Parameters: + + Input, int NP, the number of particles. + + Input, int ND, the number of spatial dimensions. + + Input/output, double POS[ND*NP], the position of each particle. + + Input/output, double VEL[ND*NP], the velocity of each particle. + + Input, double F[ND*NP], the force on each particle. + + Input/output, double ACC[ND*NP], the acceleration of each particle. + + Input, double MASS, the mass of each particle. + + Input, double DT, the time step. +*/ +{ + int i; + int j; + double rmass; + + rmass = 1.0 / mass; + +# pragma omp parallel \ + shared ( acc, dt, f, nd, np, pos, rmass, vel ) \ + private ( i, j ) + +# pragma omp for + for ( j = 0; j < np; j++ ) + { + for ( i = 0; i < nd; i++ ) + { + pos[i+j*nd] = pos[i+j*nd] + vel[i+j*nd] * dt + 0.5 * acc[i+j*nd] * dt * dt; + vel[i+j*nd] = vel[i+j*nd] + 0.5 * dt * ( f[i+j*nd] * rmass + acc[i+j*nd] ); + acc[i+j*nd] = f[i+j*nd] * rmass; + } + } + + return; +} diff --git a/src/tests/tng_io_testing.c b/src/tests/tng_io_testing.c index 2e7f9ed..c45668f 100644 --- a/src/tests/tng_io_testing.c +++ b/src/tests/tng_io_testing.c @@ -1,7 +1,7 @@ #include <inttypes.h> #include <stdlib.h> #include <string.h> -#include "tng_io.h" +#include <tng_io.h> #include "tng_io_testing.h" static tng_function_status tng_test_setup_molecules(tng_trajectory_t traj) |